Two hundred and eleven first-episode drug-naïve SCZ patients and 145 demographically matched healthy controls (HCs) were enrolled at the Department of Psychiatry of the First Affiliated Hospital of Zhengzhou University, China. This study was approved by the Human Ethics Committee of the hospital (No. 2016-LW-17). All participants provided signed informed consent.

Patients were diagnosed by two trained psychiatrists at the department using the Diagnostic and Statistical Manual of Mental Disorders, fourth edition (DSM-IV) criteria obtained via the Structured Clinical Interview for DSM-IV (SCID) (44). Severity of psychosis was evaluated with the Positive and Negative Syndrome Scale (PANSS) (45), and a general clinical assessment was performed, including height and weight for calculation of body mass index (BMI). At the hospital visits, age, sex, history of smoking, and years of education were also collected for all participants.

On enrollment, participants were required to be free of previous treatments for any mental disorders, free of diagnoses with organic diseases (46), not pregnant or lactating, and free of any use of antibiotic or anti-inflammatory agents during the prior month; not to have have a BMI >28 kg/m²; and not to have a total PANSS score ≤ 60. Other than the PANSS score, the same exclusion criteria were applied to HCs.

Fresh fecal samples were collected from all participants between 8 a.m. and 9 a.m. and immediately stored in a freezer at −80°C for gut microbiota assay. Whole blood samples were also collected at the same time for genotyping. A fecal sample of 0.2 g was used to extract the microbiome DNA using the Cetyl Trimethyl Ammonium Bromide (CTAB)/SDS procedure (47). The extracted DNA sample was diluted to 1 ng/ul with sterile water before sequencing. The 16S rRNA gene of distinct regions (V3–V4) was amplified using a specific primer (i.e., 16S V4: 515F-806R) with barcodes. PCR was performed in reaction volume of 30 μL with 15 μL of Phusion^® High-Fidelity PCR Master Mix (New England Biolabs), 0.2 μM of forward and reverse primers, and 10 ng of template DNA. The obtained PCR products were purified using a GeneJETTM Gel Extraction Kit (Thermo Scientific^TM, USA) following the manufacturer’s protocol. Purified PCR products were sequenced using the Ion S5 TM XL platform, which generates DNA 400–600 bp long single-end reads. These raw DNA sequences were quality-controlled using the Cutadapt program (V1.9.1) (48). The outputs from Cutadapt were aligned to the SILVA reference database with default settings, and chimera sequences (49) were removed using UCHIME (50).

The UPARSE program (V7.0.1001) (51) was used to assign quality-controlled sequences to operational taxonomic units (OTUs) with a sequence similarity threshold of 97%. Based on the assigned OTUs, the five most widely used α-diversity indices were estimated by QIIME (52): Shannon (42), Simpson (42), the Abundance-based Coverage Estimator (ACE) (40), Chao1 (41), and Faith’s phylogenetic diversity metric (PD) (43). Since these indices are highly skewed, all indices (except Shannon) were transformed and standardized to a mean of 0 and standard deviation of 1. For the ACE, Chao1, and PD indices, the natural log transformation was further applied. For the Simpson index, which falls within the range of 0–1, the logit transformation was performed. For the Shannon index, no transformation was applied.

DNA was extracted from whole blood samples and was genotyped using the InfiniumOmniZhongHua_8v Array (Illumina, Inc.) on an iScan instrument at Bio Miao Biological Technology (Beijing). Genotype calling and quality control (QC) were performed using GenomeStudio V2011. Quality assessments resulted in the exclusion of 23 samples from subsequent analyses due to a low call rate, i.e., <95%. The resultant file was exported to PLINK (53) format from the GenomeStudio platform. The following PLINK parameters were used for further QC before imputation: −maf 0.01, −mind 0.1, −geno 0.05, and –hwe 0.0001. In addition, variant strands were checked with the command –flip, and the concordance between self-reported and genetically predicted sex was checked with the command –check-sex. These pre-imputation high-quality data were also used for checking relatedness between subjects after construction of a thinned dataset via PLINK with the command –indep-pairwise 1,500 150 0.2 –maf 0.05, followed by –Z-genome –min 0.06. After QC, 287 subjects remained. Among these subjects, gut microbiota data were available for 159 patients and 86 HCs. The program genipe (54) was used to call IMPUTE2 (55) and Shapeit2 (55) for imputation of our sample to the 1,000 Genomes samples (56). For post-imputation QC, single nucleotide polymorphisms (SNPs) with imputation r² < 0.8, MAF <0.05, or Hardy–Weinberg test <1×10⁻⁶ were removed. In total, 4,896,670 SNPs remained. These remaining SNP dosages were hard-coded to the best-guess genotype, i.e., 0, 1, or 2 copies of the reference allele.

Subtle population structures within our sample were estimated using the principal component analysis (PCA) methods implemented in PLINK. Before PCA, pre-imputation quality-controlled genotypes were further processed using the PLINK commands –maf 0.1 and –indep-pairwise 100 50 0.1. The remaining SNPs were included in the PCA estimation with the PLINK command –pca. Since the eigenvalue scree plot showed a very flat pattern, we included the top 10 principal components (PC1-10) in our analysis.

The imputed best-guess genotypes were used for construction of polygenic risk scores. SCZ GWAS summary statistics for both European and East Asian populations were downloaded from https://www.med.unc.edu/pgc/download-results/scz/ with permission. After filtering out SNPs that did not exist in our data from the summary statistics, the program PRC-CS (57) was used to compute the European-based polygenic risk score and, separately, the East-Asian-based polygenic risk score for SCZ. Following the recommendation from the program, only SNPs from the HapMap3 reference data were used for PRS computation. Since our sample size was small, the polygenic parameter phi was preset to 0.01, i.e., assuming a highly polygenic genetic architecture for SCZ. For the remaining parameters of PRC-CS, the default values were used. The estimated effect sizes from PRS-CS were used to compute PRS for our sample using the –score function in PLINK.

The relationships of demographic and lifestyle variables with SCZ were assessed individually using t-tests (numeric variables: age, BMI, and years of education) or Fisher’s exact test for independence (sex and smoking). BMI was assessed both for the whole sample and in three subgroups: less than 18.49 kg/m²; greater than or equal to 18.5 but less than 25.49 kg/m²; and greater than or equal to 25.5 but less than or equal to 28 kg/m².

The associations of PRS and α-diversity with SCZ were tested via logistic regression models. The SCZ diagnostic status of subjects was set as the binary dependent variable. PRS was set as the predictor, along with age, sex, years of education, whether the subject had ever smoked, and PC1-10 as covariates. For the association between α-diversity and SCZ, the covariates PC1-10 were replaced by BMI. To estimate the joint effects of PRS and α-diversity indices, both α-diversity indices and PRS were entered into the model, and both PC1-10 and BMI were included as covariates, along with age, sex, years of education, and whether the subject had ever smoked, following the recommendation given by Keller (58). In addition, a term for the interaction between PRS and α-diversity indices was included to investigate whether α-diversity moderates the effects of PRS on SCZ or whether the effects of α-diversity depend on PRS values. Correction for multiple tests was applied using the false discovery rate (59), and corrected p-values <0.05 were considered statistically significant. The variance at the observed scale (Nagelkerke-R2) explained by PRS was estimated by first estimating Nagelkerke-R2 explained by the full model and the reduced model, i.e., excluding the focal variable. Next, the difference between the two estimates was used to determine the variance explained by the focal variable. This procedure was performed separately for PRS-EAS and PRS-EUR and for each of the five α-diversity indices.

Linear regression models were used to evaluate the effect of α-diversity on the symptom severity of SCZ patients, as measured by PANSS subscale scores: positive (PANSS-P), negative (PANSS-N), general psychopathology (PANSS-G), and total (PANSS-T). These subscale scores were first scaled to achieve normal distributions and were then used as the dependent variables in the models. Each subscale was analyzed separately. Sex, BMI, years of education, smoking habits, and age were included as covariates in each model. Correction for multiple tests was applied using the false discovery rate (59), and corrected p-values <0.05 were considered statistically significant.

To assess the contributions of phylum- and genus-level microbial abundance to α-diversity, linear regression models were used. Here, the transformed diversity indices, except for the Shannon index, were used as dependent variables. The relative abundances of each taxon, separately, were used as predictors, along with the covariates of sex, BMI, years of education, smoking habits, and age. Correction for multiple tests was applied using the false discovery rate (59), and corrected p-values <0.05 were considered statistically significant.

As a coarse screen analysis, logistic regressions were used to identify taxa associated with SCZ risk. The SCZ diagnosis status was used as the dependent variable, and the relative abundance was used as a predictor, along with sex, BMI, education years, smoking habits, and age as the model covariates. Correction for multiple tests was applied using the false discovery rate (59), and corrected p-values <0.05 were considered statistically significant. Taxa that were associated with SCZ were evaluated for interaction effects with PRS by adding a PRS-by–taxon interaction term to these models. In these models, the covariates were additionally expanded to include the top 10 PCs. Correction for multiple tests was applied using the false discovery rate (59), and corrected p-values <0.05 were considered statistically significant.

The PRS-related analyses were performed separately for PRS constructed using GWAS based on East Asian and European samples. Correction for multiple tests was also performed separately for these two sets of analysis.

A total of 159 first-episode, drug-naïve SCZ patients and 86 demographically matched HCs were included in the analyses. SCZ patients and HCs did not differ on age, sex, or BMI. For BMI, there was also no effect in each subgroup. Patients had significantly fewer years of education than HCs (p = 1.34×10⁻⁴, Table 1). The rarefaction plot (Supplementary Figure S1) showed that our preprocessing procedure for the microbiome data was reliable. As shown in Supplementary Figures S2, S3, that the transformed α-diversity indices became less skewed and would not be problematic in our statistical analysis.

We confirmed significant associations between PRS and the diagnosis in first-episode, drug-naïve patients, with the East-Asian sample-based PRS (PRS-EAS) showing a stronger association (OR = 2.08, 95% CI = 1.51–2.92, p = 1.22×10⁻⁵) than the European sample-based PRS (PRS-EUR) (OR = 1.73, 95% CI = 1.30–2.35, p = 2.77×10⁻⁴). On the observed scale, PRS-EAS and PRS-EUR explained 8.9 and 5.8% of phenotypic variance, respectively. Stratifying participants into quintiles by PRS showed that the proportion of patients increased monotonically from Q1 to Q5 (Figure 1). In the highest stratum (Q5), 82.5 and 78.9% of participants were patients based on PRS-EAS and PRS-EUR, respectively. The corresponding numbers for the lowest stratum (Q1) were 51.7 and 56.9%, respectively. These results suggest that SCZ may have both shared and unique genetic architecture in East Asian and Caucasian populations confirming that matching the ancestry of GWAS samples used to generate a PRS with that of target prediction samples could improve risk prediction.

We discovered distinct patterns of association between the five α-diversity indices and SCZ (Figure 2). While all pairwise correlations among the five indices were significant (p < 0.05), these indices belong to three groups (Figure 2A). The ACE and Chao1 indices, both of which are richness-based measures, were almost identical (Pearson r = 0.99) but were less correlated to the Shannon and Simpson indices, both of which are evenness-based measures. The latter two showed a significant and strong correlation (Pearson r = 0.93). The phylogenic diversity-based index, PD, sits in between the two groups. The differences between these three groups of indices were also reflected in their associations with the risk of SCZ. An increased score on the Shannon and Simpson indices showed a negative association with SCZ (Shannon: OR = 0.29, 95% CI = 0.18–0.43, p = 1.15×10⁻⁸; Simpson: OR = 0.29, 95% CI = 0.19–0.44, p = 1.25×10⁻⁸). The proportion of phenotypic variance explained by both the Shannon and the Simpson indices was 33%. The PD index showed a weaker association with SCZ (OR = 1.43, 95% CI = 1.05–2.00, p = 0.03, Nagelkerke-R2 = 19%). Neither the ACE nor the Chao1 index was associated with the risk of SCZ (Figure 2B).

Moreover, α-diversity was also associated with symptom severity. Both the Shannon and Simpson indices were inversely associated with total, positive, and general psychopathology PANSS scores (p < 0.05; PANSS-T: Shannon, beta = −0.31, se = 0.1; Simpson, beta = −0.33, se = 0.09; PANSS-P: Shannon, beta = −0.25, se = 0.1; Simpson, beta = −0.28, se = 0.11; PANSS-G: Shannon, beta = −0.35, se = 0.1; Simpson, beta = −0.33, se = 0.11). Higher diversity as measured by these two indices was associated with lighter symptoms on these three subscales. There were also suggestive associations between PD and PANSS-G (beta = −0.29, se = 0.14) and PANSS-T (beta = −0.28, se = 0.13), but neither of these survived correction for multiple tests. There were no associations between α-diversity and PANSS-negative symptoms.

To understand whether gut microbial diversity contributes additional information to PRS in predicting SCZ, we entered both diversity index and PRS as predictors in multiple logistic regression models (Materials and Methods section). In general, we observed slightly changed effects for both PRS and the Shannon index compared with the effects estimated via the univariate logistic models. However, all significant predictors from the corresponding univariate models were still significant in the multiple regression models. For example, in the model with the Shannon index and PRS-EAS, the effect of the Shannon index changed from on OR of 0.29 to 0.26, and that of PRS-EAS changed from an OR of 2.08 to 2.27. In the model including both PRS-EAS and the PD index, the effect of PRS-EAS was slightly reduced (OR = 1.97), and the weaker effect of the PD index was reduced to null (p = 0.14). In addition, there were no qualitative changes to the effects of PRS-EAS or the Chao1 index in these models compared to those described above. Replacing PRS-EAS with PRS-EUR produced the same pattern. Therefore, the Shannon and Simpson diversity indices may contribute additional information beyond PRS in explaining the risk of SCZ.

Next, we examined whether the effects of α-diversity indices on SCZ were dependent on the genetic risk of participants as measured by PRS. We added terms representing the interactions between α-diversity indices and PRS to the above-described multiple logistic regression models. As illustrated in Figure 3, we found significant interactions between PRS-EAS and the Shannon, Chao1, and PD indices (Shannon, p = 0.05; Chao1, p = 4.43×10⁻³; and PD, p = 6.31×10⁻³). These interaction terms explained 1.05, 3.15, and 2.91% of the phenotypic variance, respectively. There was no interaction effect detected between PRS-EAS and the Simpson index (p = 0.11). Interestingly, no significant interaction effects between PRS-EUR and the five α-diversity indices were found.

To illustrate the nature of these significant interaction terms, we stratified our samples by the median values of PRS and, separately, of the α-diversity indices (Figure 3 and Supplementary Figures S4–S7). The effects of PRS-EAS on SCZ risk were larger for participants having higher α-diversity (Figures 3A–C), as measured by the Shannon, Chao1, and PD indices. In other words, when participants were stratified by PRS-EAS, we observed that increasing the α-diversity index measure reduced SCZ risk for participants with lower PRS-EAS. For those with higher PRS-EAS, increasing α-diversity according to the Chao1 and PD indices also increased the SCZ risk (Figures 3D–F). Similar results were also obtained for the Simpson index (Supplementary Figures S4, S5). The ACE index showed similar results to the Chao1 index (Supplementary Figures S6, S7).

To test whether the interaction effect we observed could be due to correlations between PRS-EAS and α-diversity (i.e., the possibility that higher diversity may correlate with lower PRS-EAS), we performed a linear regression analysis of each α-diversity index against PRS-EAS based on HC samples (N = 86). In these models, age, sex, BMI, whether the participant had ever smoked, years of education, and PC1-10 were included as covariates. None of the five diversity indices were correlated with PRS-EAS in our sample (Supplementary Table S1). Thus, the effects of PRS-EAS on SCZ seem to interact with gut microbial α-diversity.

To gain a deeper understanding of the observed interaction effects, we tested whether there were interactions between taxon abundance and PRS. Before performing these tests, we first assessed which taxa contributed to α-diversity and SCZ risk (Materials and Methods section). There were four genera contributing to the Shannon index (three genera: unidentified_Lachnospiraceae, beta = 0.20, se = 0.05, p = 4.43×10⁻²; Megamonas, beta = −0.21, se = 0.05, p = 1.52×10⁻²; Agathobacter, beta = −0.3, se = 0.05, p = 1.71×10⁻⁶), Simpson index (one genus: Agathobacter, beta = −0.34, se = 0.04, p = 3.51×10⁻¹³), and PD (one genus: Enterobacter, beta = 0.26, se = 0.06, p = 6.35×10⁻³) (Figure 4A). No specific taxa contributed to the ACE or Chao1 indices. There were 10 taxa associated with SCZ risk, among which six showed a protective effect, i.e., higher abundance was associated with a lower risk, and four showed a negative effect (Figure 4A and Supplementary Table S2). While there was no interaction effect between the genus Agathobacter (which was associated with both α-diversity and SCZ) and PRS, three genera showed significant interactions with PRS in their effect on SCZ risk (Figures 4C,D): Romboutsia (interaction beta = 0.84, se = 0.3, p = 0.03), Streptococcus (interaction beta = 0.89, se = 0.34, p = 0.03), and Anaerostipes (interaction beta = 1.37, se = 0.51, p = 0.03). For these three genera, increased abundance was associated with greater protective effects against SCZ in individuals with a lower-than-median PRS than those having above-median PRS (Figures 4C,D). We obtained qualitatively the same interaction effects between these three genera and PRS-EUR (Supplementary Table S3).