Genome-wide association study followed by trans-ancestry meta-analysis identify 17 new risk loci for schizophrenia

SCZ cases were from inpatient and outpatient services of collaborating mental health centers of China. Part of cases has been described in our previous studies for candidate gene analyses [32–34]. Diagnosis was based on Diagnostic and Statistical Manual of mental disorders (DSM-IV) criteria, with the use of Structured Clinical Interview for DSM-IV (SCID) Axis I Disorders. All relevant and detailed information (including onset of SCZ, first onset or remitted, symptoms and chief complaint, duration course, family history of psychiatric disorders, medication history) were carefully evaluated by at least two independent experienced psychiatrists to reach a consensus DSM-IV diagnosis. Detailed information about diagnosis had been described in previous studies [33, 34]. The average age of cases and controls were 35.67 ± 10.29 and 28.82 ± 6.83 years, respectively. 37.45% cases and 54.53% controls were males, respectively. All participants provided written informed consents. This study was approved by the Ethical Committee and internal review board of the Kunming Institute of Zoology (No: SMKX-20191215-07) and participating hospitals and universities (including the Second Affiliated Hospital of Xinxiang Medical University and Xi’an Jiaotong University). The samples were recruited from 2010 to 2018.

Genomic DNA was extracted from the peripheral blood with the use of QIAamp DNA blood mini kit (Cat. No: 51106). We used two types of genotyping platforms (arrays), including Illumina ASA (BeadChip Array Asian Screening Array-24+v1.0 HTS ASAMD-24v1-0) and GSA (BeadChip Array Global Screening Array-24+v2.0 HTS GSA v2.0+Multi-Disease). For ASA array (including 743,722 variants), 56.7% of the variants are common variants (with minor allele frequency > 0.05), 30.8% are low-frequency variant (with minor allele frequency between 0.01 and 0.05), and 12.5% are rare variants (minor allele frequency < 0.01). ASA array includes a broad spectrum of pharmacogenomics markers (N = 5588) obtained from CPIC guidelines (www.cpicpgx.org↗) and the PharmGKB database (www.pharmgkb.org↗). In addition, the ASA array contains about 50,000 SNPs selected from ClinVar database (www.ncbi.nlm.nih.gov/clinvar↗). For GSA array (N = 759,993 variants), 54.4% of the variants are common variants (with minor allele frequency > 0.05), 18.1% are low-frequency variant (with minor allele frequency between 0.01 and 0.05), and 27.4% are rare variants (minor allele frequency < 0.01). Similar with ASA array, the GSA array includes multiple expert curated variants obtained from ClinVar (www.ncbi.nlm.nih.gov/clinvar↗), PharmGKB (www.pharmgkb.org↗), NHGRI (www.genome.gov/↗), and other databases. More details about the ASA and GSA arrays can be found in the official Illumina website (https://www.illumina.com/products/by-type/microarray-kits↗). Genotyping assays were conducted at Guoke Biotechnology Co., LTD in Beijing (www.bioguoke.com↗).

We conducted quality control as previously described [23, 24, 35], with some minor revisions. The main consideration is to include more subjects and variants under the premise of ensuring data quality. The individual-level quality control (QC) was processed as follows: (1) Samples with genotype missing rate > 0.03 were excluded; (2) Samples with a heterozygosity (calculated by plink 1.9 –het command) deviated ± 3 stand deviation (s.d.) from the mean heterozygosity of all samples were excluded; (3) Related samples were inferred by KING software (http://people.virginia.edu/~wc9c/KING/↗) [36]. We used –related command implemented in KING to infer the potential kinship coefficients with 3rd degree. After inferring relatedness between samples, KING program exported related and unrelated samples. Related samples were then excluded. (4) We checked the sex of samples by using Plink (v1.09) (with --check-sex command). Samples with inconsistent sex information were excluded. We excluded samples if their sex estimated by genotype data and medical record information were inconsistent. In addition, we also excluded samples whose sex could not be accurately predicted based on the genotype data.

The SNP-level QC are as follows: (1) Excluding SNPs with a call rate < 97%; (2) As described in the study of Lam et al. [24], we excluded SNPs with a significant deviation from Hardy-Weinberg equilibrium (P < 1.0 × 10⁻⁶ in controls and P < 1.0 × 10⁻¹⁰ in cases); (3) Excluding SNPs with a minor allele frequency (MAF) < 0.01; (4) Only biallelic SNPs were retained for further analysis. The QCs were performed by Plink 1.9 software [37]. Our analysis started in March 2020 and finished in May 2021.

We performed principal component analysis (PCA) to assess population stratification and exclude the outliers. The subjects from the 1000 Genomes project (including CHB (Han Chinese in Beijing, China), CHS (Southern Han Chinese), JPT (Japanese in Tokyo, Japan), YRI (Yoruba in Ibadan, Nigeria), and CEU (Utah residents with northern and western European ancestry)) [38] were used as references and PCA was performed with GCTA software [39]. Samples that were not clustered with subjects of Han Chinese ancestry were excluded. We then used the Smartpca program to calculate principal components as covariates to correct potential population stratification [40, 41]. We excluded the MHC region (chr6: 25-34 MB) when performing PCA. Five PCA iterations were run and top 20 principal components (PCs) were calculated. Samples that were away from 6 standard deviations (s.d.) of the mean of each PC were excluded. After stringent quality control, 3493 cases and 4709 controls were retained for GWAS.

For each group (ASA subgroup 1 and subgroup 2, GSA), the top 20 PCs were calculated. The PCs were calculated as covariates to correct potential population stratification for each group (ASA subgroup 1 and subgroup 2, GSA), respectively. As suggested by Price et al. [41], we included a variable number of PCs as covariates to perform logistic regression analysis in each group (ASA subgroup1, ASA subgroup2, GSA). We used genomic control inflation factor (λ_GC) to estimate the effect of the number of PCs on population stratification adjustment [41, 42]. The selection of the final number of PCs was mainly based on the following criteria [43]: (1) the λ_GC should as close to 1 as possible (indicating population stratification is well controlled); (2) when including more PCs, the λ_GC did not change significantly (indicating that the number of PCs selected is enough to adjust population stratification). Based on these criteria, we finally selected 10 PCs for ASA subgroup 1, 4 PCs for ASA subgroup 2, and 12 PCs for GSA samples respectively.

The imputation of the genotypes was performed by Eagle [44] and minimac3 [45]. The Eagle was used for phasing the genotype data of each chromosome. The reference panel from the 1000 Genomes project (Phase 3) [38] was downloaded from the minimac3 website (https://genome.sph.umich.edu/wiki/Minimac3↗). The imputed data were further processed with the following QC steps: (1) SNPs with an imputation quality score < 0.8 were excluded; (2) SNPs with a MAF < 0.01 were excluded; (3) SNPs that significantly deviated from the Hardy-Weinberg equilibrium (P < 1 × 10⁻⁶ in controls, P < 1 × 10⁻¹⁰ in cases) were excluded; (4) Only SNPs with a genotype imputation rate > 97% were remained; (5) Only biallelic SNPs were included for further association analysis. After QC, the total number of SNPs included in GWAS analysis was 4,724,225 (for ASA arrays) and 5,107,135 (for GSA arrays). The number of overlapping SNPs between the two platforms is 3,937,527. To make our QC procedures more clear and easy to follow, we provided the detailed information about QC (to list all pre- and post-imputation QC steps with excluded and remaining SNPs and participants after each step) in Additional file 1: Figure S1.

As our samples were genotyped with Illumina ASA and GSA Chip arrays, to avoid potential effects of different array platforms on our analysis, we performed genome-wide association analysis separately. We firstly performed genetic association analysis in our samples (ASA subgroup 1, ASA subgroup 2, GSA, respectively) by using logistic regression. This association analysis was based on genotype data of each group (i.e., ASA subgroup 1, ASA subgroup 2, GSA), with PCs included as covariate. The results from these three genome-wide summary statistics were then meta-analyzed (summary statistics-based meta-analysis, fixed-effect model was used). GWAS were performed by Plink (v1.09) software [37]. We further performed a summary statistics-based meta-analysis in East Asian populations through combining the results of our study and GWAS summary statistics (only EAS were used) from a recent study by Lam et al. (including 22,778 case and 35,362 controls) [24]. Finally, we carried out a trans-ancestry meta-analysis (summary statistics-based, fixed effect model) through combining our results with the summary statistics from East Asians and Europeans (containing 56,418 cases and 78,818 controls) [24, 26].

To identify independent risk loci, we used FUMA [46] to clump the association results (with the default parameters). Different reference panels were used for clumping. For the meta-analysis in EAS, subjects of East Asian ancestry from the 1000 Genomes project [38] were used. For the meta-analysis of EAS and PGC2 [24], subjects of European ancestry from the 1000 Genomes project [38] were used. The clumping processes were as follows: The minimum r² threshold for independent significant SNPs was set to 0.6, which was used to define the boundaries of the genomic risk locus. When independent significant SNPs were defined in a locus, these SNPs were further processed to identify the lead SNP for each locus. The minimum r² for defining lead SNPs was set to r² = 0.1. The LD blocks of independent significant SNPs that < 250 kb were merged into a single genome locus. More details about LD clump can be found in FUMA website (https://fuma.ctglab.nl/tutorial↗). Visualization of association results of interest SNPs and its nearby variants were generated with locuszoom [47] (http://locuszoom.org/genform.php?type=yourdata↗).

We compared the genome-wide significant (GWS) loci identified in this study and published SCZ GWAS [21, 23, 24, 26], and we also included GWASs listed in GWAS catalog database in FUMA [46]. We used following approaches for risk loci comparison: (1) If the published SCZ GWAS provide GWS index SNPs and its corresponding genomic region (chromosomal coordinates), we overlapped genomic region of our GWS loci with these published risk loci. Non-overlapping loci were defined as new loci; (2) If the published SCZ GWAS only list GWS index SNPs, the newly identified risk loci (by us) should have no overlap with these GWS index SNPs, the index SNPs of the newly identified risk loci (by us) should also not be in linkage disequilibrium with the reported GWS index SNPs (R² < 0.1); (3) GWAS catalog database implemented in FUMA [46] were also used to identify new SCZ GWS loci.

Polygenetic risk score profiling was performed by PRSice2 [48]. The PRS scores derived from training datasets were used to predict the case-control status of our samples. And the samples included in this study were used as target sample. We used classic clumping and threshold method. We run PRS using default parameters in PRSice2 [48]. The detailed clumping parameters of PRS calculation were as follows: --clump-kb 250, --clump-p 1.0, --clump-r2 0.10. Summary statistics from EAS [24], PGC2 EUR [24], EAS + EUR [24], and CLOZUK+PGC2 [21] were used as training datasets in PRS analysis. And we set 10 P value thresholds (5 × 10⁻⁸, 5 × 10⁻⁵, 0.0001, 0.001, 0.01, 0.05, 0.1, 0.2, 0.5, 1) to analyze the phenotypic variance explained at different P value cutoffs.

To explore if schizophrenia associations were enriched in specific human tissues, we performed tissue enrichment analysis by using MAGMA (Version 1.08) [49], which is implemented in FUMA software (Version 1.3.6a) [46]. Briefly, gene expression data of different human tissues (RNA sequencing data from the Genotype-Tissue Expression (GTEx) [50] consortium) were used to identify the genes that differentially expressed in a specific tissue. Based on the GWAS P values, MAGMA quantifies the degree of association between a gene and schizophrenia (i.e., obtain a gene-level P value) by using a multiple linear principal component regression model. MAGMA then tests if schizophrenia associations were enriched in the specifically expressed genes in a specific tissue. More detailed information about tissue enrichment analysis can be found in FUMA website (https://fuma.ctglab.nl/↗).

Cell-type enrichment analysis was also performed using MAGMA (Version 1.08) [49]. Single-cell RNA sequencing data from the mouse central nervous system (CNS) [51] were downloaded and processed as described in Bryois et al. [52]. Briefly, a total of 160,769 high-quality single-cell RNA-seq data were analyzed. Human genes were mapped to orthologous mouse genes (one to one) based on MGI annotations (http://www.informatics.jax.org/homology.shtml↗), and genes that were not expressed in the mouse central nervous system (CNS) were excluded. The expression specificity of a specific gene was calculated as described in Bryois et al. [52] and the top 10% that were most specifically expressed genes in each cell type were used for enrichment analysis. The P values of MAGMA enrichment analysis were corrected by false discovery rate (FDR).

Gene set analysis was performed with MAGMA software (Version 1.08) [49]. MAGMA gene set analysis includes two main procedures. Firstly, by using the GWAS summary statistics and gene annotation files (NCBI b37, downloaded from official MAGMA website (https://ctg.cncr.nl/software/magma↗)), MAGMA maps SNPs to genes (with gene boundaries setting parameters “--annotate window = 35,10”). MAGMA then calculates the association strength between a specific gene and schizophrenia. A gene and phenotype association matrix was generated in this step with default “snp-wise = mean” model. The MHC region (from 25 to 34 Mb) on chromosome 6 was excluded in analysis. Secondly, gene set analysis which was also known as competitive gene set analysis was performed. MAGMA tests whether the association of a target gene set with phenotype is greater than other genes that are not included in the gene set. The P value of competitive gene set was used to determine the significance level. We downloaded gene sets from MSigDB database (v7.1) (http://software.broadinstitute.org/gsea/msigdb/↗) [53] and all GO terms (including cellular component, biological processes, and molecular functions). KEGG pathway gene sets were also included and a total of 10,378 gene set terms were compiled. We further retained 6194 gene sets with a gene number range from 10 to 200 for the final analysis. The final P values (i.e., the association strength between gene sets and schizophrenia) of MAGMA competitive gene sets were corrected by false discovery rate (FDR).

We conducted stratified linkage disequilibrium score regression (LDSC) [54] to partition SCZ SNP heritability and test enrichment of SCZ heritability in different functional annotations [54]. In addition, we also performed genetic correlation analysis between our Chinses cohort and the reported EAS samples (based on summary statistics) using LDSC [54].

We performed TWAS using the summary statistics from all combined samples (including our Han Chinese cohort, EAS and PGC2 [24], a total of 59,911 cases and 83,527 controls) and brain eQTL (gene expression SNP weights) from the CommonMind Consortium (CMC, N = 452) [55]. CMC measured gene expression in the dorsolateral prefrontal cortex (DLPFC) of human brain with the use of RNA sequencing. The CMC gene expression SNP weights were derived from the FUSION pipeline (http://gusevlab.org/projects/fusion/↗). Detailed information about sample collection, RNA extraction and sequencing, genotyping and statistical analysis of CMC dataset can be found in the original studies [55]. TWAS was performed on autosomal chromosomes using the FUSION.assoc.test.R script across all predictive models, such as LASSO, GBLUP, Elastic Net, and BSLMM. To determine which is the best model for TWAS, FUSION performed a five-fold cross-validating of each model. TWAS associations (i.e., genes) were considered “transcriptome-wide significant” if they passed a strict Bonferroni-corrected threshold for all genes tested in the dataset (corrected significance P value: 0.05/3551 = 1.41 × 10^{− 5}). Detailed description of the principle of FUSION and statistical model can be found in the original paper [56].

We performed eQTL analysis using 4 public available brain eQTL resources, including the CommonMind Consortium (CMC) [55], LIBD BrainSeq Phase II RNA-seq Project (LIBD2) [57], xQTL [58], and GTEx [50] project. Please refer to the original papers [50, 55, 57, 58] about the sample collection, RNA-seq data processing, and eQTL calculation.

We examined the expression level of genome-wide significant genes in brains of SCZ cases (n = 559) and controls (n = 936) by using the PsychEncode data [59]. Briefly, brain gene expression data of 559 schizophrenia cases and 936 controls were quantified and analyzed by the PsychEncode. Detailed information on tissue collection, RNA sequencing, gene expression quantification, and differential expression analysis were provided in PsychEncode papers [59] and website (http: http://resource.psychencode.org/↗).

Our samples had no overlap with the previously published SCZ GWAS of Han Chinese population [17, 19, 23, 60]. We first performed a PCA analysis using the samples genotyped with Illumina Asian Screening Arrays (ASA) and found population stratification of our samples (Fig. 1a); we thus divided our samples into two genetically matched subgroups. After stringent QC and excluding outliers, 2055 cases and 1823 controls were included in subgroup 1, and 607 cases and 1186 controls were included in subgroup 2 (Fig. 1b, c). For samples genotyped with GSA arrays, no obvious population stratification was observed. After strict QC, the final samples genotyped with GSA arrays included in this study were 831 cases and 1700 controls (Additional file 1: Figure S2). PCA analysis using genotype data of our samples and the 1000 Genomes project [38] showed that all of our samples clustered with samples of Han Chinese ancestry (Additional file 1: Figures S3, S4). After strict QC, imputation using the 1000 Genomes project phase 3 panel [45] and post-imputation QC, a total of 3,937,527 biallelic SNPs from 3493 SCZ cases and 4709 healthy controls were retained for GWAS. Principal components (PCs) were included as covariates [40, 41] (10 PCs for ASA subgroup1, 4 PCs for ASA subgroup2, 12 PCs for GSA samples) when performing GWAS. We firstly meta-analyzed the samples genotyped with ASA arrays (subgroups 1 and 2). We then conducted a meta-analysis through combining the results from ASA and GSA arrays with the fixed effect model. The genomic inflation (λ_GC) of our combined meta-analysis (including ASA and GSA samples) was 1.10, and the λ₁₀₀₀ (scaled to a sample size of 1000 cases and 1000 controls) was 1.02 (which was very close to the reported values in previous Chinese GWAS [17, 23], 1.02 in Li et al. [23] and 1.019 in Yu et al. [17]) (Additional file 1: Figure S5), indicating that population stratification unlikely confounds our genome-wide association results.

Two genome-wide significant risk loci were identified in our sample (Fig. 1d). The lead risk SNP for the first locus is rs7192086 (P = 4.92 × 10^-8, OR = 1.22), which is located in the intron 2 of the SHISA9 gene (Fig. 1e). Of note, a previous study (7308 cases and 12,834 controls) showed nominal association between rs7192086 and SCZ (P = 1.34 × 10^{− 05}, OR = 1.12) in European population [22]. The lead SNP for the second risk locus is rs57016637 (P = 2.33 × 10⁻¹¹, OR = 1.34) located in the intron 2 of PES1 (Fig. 1f). This genomic loci contain three independent significant SNPs, and the other two SNPs are rs117961127 (P = 2.73 × 10⁻¹¹, OR = 1.35) (located in the intron 2 of OSBP2) and rs116976860 (1.42 × 10⁻⁰⁹, OR = 1.33) (located in the intergenic region of SEC14L6 and GAL3ST1). Of note, rs57016637 is an East Asian-specific polymorphism (Additional file 1: Figure S6).

We then carried out a meta-analysis through meta-analyzing GWAS data obtained from our sample and a recent GWAS conducted in EAS [24]. Our meta-analysis using combined EAS samples (26,271 cases, 40,071 controls) [24] identified 24 risk loci (Additional file 1: Figure S7). However, all of these loci have been reported in a previous study [24]. Intriguingly, we noticed that rs3845188 (P = 6.50 × 10⁻⁰⁸, OR = 0.91) (which is located in the intron 1 of the NEBL gene) reached the suggestive significance level (P = 5.00 × 10⁻⁰⁶), suggesting this locus may be associated with SCZ (Additional file 1: Figure S8).

We further performed a meta-analysis through combining results of our study and GWAS results from EAS and PGC2 European samples (59,911 cases and 83,527 controls) [24]. We identified 153 genome-wide independent risk loci in the combined samples (Fig. 2a) (The definition of independent risk loci was described in methods section). Among these 153 risk loci, 15 were novel (Table 1, Additional file 1: Table S1). In total, we identified 17 new risk loci for SCZ.

To identify the potential target genes of the newly identified risk SNPs, we performed eQTL analysis using data from the human brain tissues. Among the 17 lead SNPs, 9 showed associations (uncorrected P < 0.05) with expression of 33 genes in the human brain (Additional file 1: Table S2). Of note, rs2106747 (P = 3.36 × 10⁻⁰⁸, OR = 1.05, Fig. 2b) showed a strong association with FAM221A expression in all four eQTL datasets. Another interesting SNP is rs7301566 (P = 2.91 × 10⁻⁰⁸, OR = 1.06, Fig. 2c), which was associated with expression of several genes, including COX14, DIP2B, CERS5, RP4-605O3.4, SPATS2, ASIC1, and ATF1 (Additional file 1: Table S2). These results suggest that the newly identified risk variants may contribute to SCZ risk through regulating expression of these eQTL genes.

We explored the expression level of potential target genes of the newly identified risk variants in SCZ cases and controls using expression data from the PsychEncode (including 559 cases and 936 controls) [59]. Among the 33 potential target genes, 14 showed nominal difference in expression in SCZ cases compared with controls (Additional file 1: Table S3), supporting that the newly identified risk variants may confer SCZ risk by regulating the expression of these target genes. In addition, we also examined the expression of genes located near the two newly identified loci in our Chinese samples (Table 1, Fig. 1e, f). We found that SHISA9 showed a trend of upregulation in brains of SCZ cases compared with controls (P = 0.053). Interestingly, OSBP2 was significantly downregulated in brains of SCZ cases compared to controls (P = 1.07 × 10⁻⁰⁷). Taken together, these expression data provide further evidence that support the newly identified risk variants may confer risk of SCZ through modulating expression level of these genes.

We conducted PRS analysis to predict the case-control status of our samples (ASA subgroup 1) and estimate the phenotypic variance (of our samples) that can be explained by the published GWAS summary statistics data [24]. When using the summary statistics from EAS as training set [24], the explained variance (estimated by Nagelkerke R²) ranged from 0.4 to 5.9% and the training data has the largest variance explanation at P value = 0.2 (P = 2.19 × 10⁻³⁷) (Fig. 3). The training set from EAS + EUR has an overall better prediction performance and variance explanation at each P value threshold than EAS [24], and the explained variance ranged from 2.0% to 6.5%. At P value = 0.01 threshold, the training dataset has the largest variance explanation (P = 4.69 × 10⁻⁴¹). The EUR and CLOZUK+PGC2 training sets had relatively poor prediction performance than the EAS and EAS + EUR training set. The explained variance (estimated by Nagelkerke R²) ranged from 0.6 to 2.7% and 0.7 to 4.4% for EUR and CLOZUK+PGC2 training datasets, respectively. The PRS analysis indicated that EAS and EAS + EUR training sets had relative good power to predict the SCZ and healthy controls status of our sample, and the EAS + EUR training set had better prediction performance.

We conducted tissue enrichment analysis by using GWAS associations from the combined samples (i.e., our cohort, EAS and PGC2 EUR [24]) and gene expression from GTEx [50], with the use of FUMA [46]. SCZ heritability was mainly enriched in brain tissues (Additional file 1: Figure S9a). Of note, two brain cerebellum tissues showed the strongest associations (brain cerebellar hemisphere, P = 2.54 × 10⁻²³, and brain cerebellum, P = 8.51 × 10⁻²³). The frontal cortex (BA9) (P = 1.28 × 10⁻¹⁹) and brain cortex (P = 2.20 × 10⁻¹⁸) also showed significant enrichment.

We further conducted cell-type enrichment analysis (Additional file 1: Figure S9c). Similar with previous findings [52], two telencephalon neuronal cell types, including telencephalon projecting excitatory interneurons (P = 1.12 × 10⁻⁰⁹, FDR < 0.05) telencephalon projecting inhibitory interneurons (P = 1.15 × 10⁻⁰⁶, FDR < 0.05), showed significant enrichment for SCZ associations. We also identified other novel associations, including oligodendrocytes (P = 8.03 × 10⁻⁰⁷, FDR < 0.05), cholinergic and monoaminergic neurons (P = 6.06 × 10⁻⁰⁵, FDR < 0.05), and cerebellum neurons (P = 4.99 × 10⁻⁰⁵, FDR < 0.05). Collectively, our results suggest that SCZ risk genes are actively expressed in these identified tissues and cell types, implying the pivotal roles of these tissues and cell types in SCZ.

We carried out gene set enrichment analysis and identified 6 significant enriched terms, including neuron spine (P = 7.06 × 10⁻⁰⁷, FDR = 0.0044), membrane depolarization during action potential (P = 4.72 × 10⁻⁰⁶, FDR = 0.015), cytosolic calcium ion transport (P = 4.03 × 10⁻⁰⁵, FDR = 0.045), voltage gated sodium channel complex (P = 4.09 × 10⁻⁰⁵, FDR = 0.045), t tubule (P = 4.32 × 10⁻⁰⁵, FDR = 0.045), and voltage gated sodium channel activity (P = 2.57 × 10⁻⁰⁵, FDR = 0.045). In addition, a few gene sets also showed a trend of enrichment, including regulation of synaptic plasticity (P = 1.5 × 10⁻⁰⁴, FDR = 0.074) and neurotransmitter receptor complex (P = 2.00 × 10⁻⁰⁴, FDR = 0.074) (Additional file 1: Table S4).

LDSC analysis showed that SCZ associations were mainly enriched in various conserved genomic regions, promoters, and enhancers (SuperEnhancer_Hnisz, H3K27ac_Hnisz) (Additional file 1: Figure S9b), which were consistent with previous findings [26, 54].

We conducted a TWAS to identify genes whose cis-regulated expression were associated with SCZ. We identified 76 transcriptome-wide significant genes (corrected by multiple comparison testing) (Fig. 4 and Additional file 1: Table S5). Among these identified genes, 24 have been reported in a previous study [59]. Of note, GIGYF1 (P = 3.03 × 10⁻⁰⁵) (located near the newly identified risk loci) showed a trend of TWAS significance (Additional file 1: Table S5). Taken together, our TWAS identified 52 new risk genes whose cis-regulated expression change may have a role in SCZ.

All participants provided written informed consents and this study was approved by the Ethical Committee and internal review board of the Kunming Institute of Zoology, the Second Affiliated Hospital of Xinxiang Medical University and Xi’an Jiaotong University. Ethics approval No: SMKX-20191215-07.

Not applicable.

The authors declare no competing interests and financial relationships with commercial interests.

The datasets generated and/or analyzed during the current study are not publicly available due to the regulations of the Ministry of Science and Technology of the People’s Republic of China (on human genetic resource usage). We need to comply with relevant regulations regarding the release of original genetic data (including the summary statistics) of the Chinese patients. This will take some time to finish according to the related policy. Data depositing in a third-party website (without application) is not allowed before approval. We understand that sharing of the genome-wide summary statistics publicly will contribute to the progress of the relevant field and we are planning to promote the progress of sharing summary statistics after this study is finalized. In addition, these GWAS samples were recruited from multiple hospitals across Mainland China, and this study could not have been completed without the great efforts from all the collaborators, we will need to get approval from all the authors (and their affiliated agencies) before sharing the summary statistics. However, these data are available from the corresponding author on reasonable request.

URLs:

PLINK, http://www.cog-genomics.org/plink2↗; KING, http://people.virginia.edu/~wc9c/KING/↗; EIGENSOFT, https://www.hsph.harvard.edu/alkes-price/software/↗; GCTA, http://cnsgenomics.com/software/gcta/;↗ The 1000 Genome project website, https://www.internationalgenome.org/↗; Minimac3, https://genome.sph.umich.edu/wiki/Minimac3↗; LDSC, https://github.com/bulik/ldsc↗; MAGMA, https://ctg.cncr.nl/software/magma↗; MSigDB, http://software.broadinstitute.org/gsea/msigdb/↗; PsychEncode, http://www.psychencode.org/↗;; Locuszoom, http://locuszoom.org↗; FUMA, https://fuma.ctglab.nl/↗. Single cell RNA-seq analysis, https://github.com/jbryois/scRNA_disease↗; PRSice2, http://www.prsice.info/↗; FUSION, http://gusevlab.org/projects/fusion/↗; SZDB, http://www.szdb.org↗