Characteristics of Gut Microbiota in Patients with GH-Secreting Pituitary Adenoma

In total, GH patients (n = 28) and healthy controls (n = 67) were enrolled in this study, with age and body mass index (BMI) adjusted. In patients with GHPA, 5 of 28 were suffering from recurrent GHPA, while the others were newly diagnosed. Levels of nadir GH and IGF-1 were 14.40 (5.12 to 22.47) ng/mL and 758.93 ± 146.10 ug/L, respectively. Clinical and biochemical findings of GHPA and healthy controls are summarized in Table 1, and the details of both groups are provided in Table S1 and S2, respectively.

Clinical parameters, including BMI, fasting plasma glucose (FPG), fasting GH, nadir GH, IGF-1, and IGF-1 index, were compared preoperation and 6 months after trans-sphenoid endoscopic assisted surgery, the mainstay of GHPA treatment (Table 2). In our medical center, we use antibiotics after surgery rather than preoperatively for this surgery, and we collected fecal samples before surgery, therefore avoiding influence on the microbiota profile. In general, except for BMI (25.87 ± 2.95 kg/m² versus 25.51 ± 3.04 kg/m², P = 0.23), significant decrease was observed in clinical parameters, including FPG (5.26 ± 0.54 mmol/L versus 4.89 ± 0.49 mmol/L, P < 0.001), fasting GH (20.57 ± 15.55 ng/mL versus 2.17 ± 3.62 ng/mL, P < 0.001), nadir GH (16.23 ± 13.33 ng/mL versus 1.41 ± 2.33 ng/mL, P < 0.001), IGF-1 (758.93 ± 146.10 μg/L versus 369.42 ± 195.17 μg/lL, P < 0.001), and IGF-1 index (2.44 ± 0.74 versus 1.16 ± 0.59, P < 0.001). Nineteen out of 28 patients reached biochemical remission, defined as an age-normalized serum IGF-1 value and a random GH of <1.0 ng/mL (30), and 20 out of 28 patients achieved radiological remission, defined as no recurrence or residual tumor on follow-up magnetic resonance imaging.

To investigate the microbial diversity between patients with GHPA and healthy controls (age and BMI adjusted), we characterized the clustering effects and found that 516 of the total 742 operational taxonomy units (OTUs) overlapped between patients with GHPA as GH group and healthy controls as HC group, while 42 OTUs were found specifically in GH group (Fig. 1A). Then, we analyzed the compositions of bacterial community between the two groups and detected the prominently altered gastrointestinal microbial diversity. Compared with that in the control group, fecal microbial diversity, assessed by Shannon index (Fig. 1B) and Simpson index (Fig. 1C), was notably decreased in GH group (P = 0.036 and P = 0.013, respectively). To demonstrate microbiome space between the groups, beta diversity was evaluated by weighted UniFrac, and intergroup variance was further presented via nonmetric multidimensional scaling (NMDS) along with P value acquired by multiresponse permutation procedure (MRPP) (Fig. 1D) and principal coordinates analysis (PCoA) along with Adonis analysis (Fig. 1E). Notable change was observed between two groups (P = 0.002 for MRPP analysis and P = 0.002 for Adonis test). Collectively, compared to that in the healthy controls, microbial diversity was significantly decreased in GH group, suggesting the potential vulnerability to various pathological conditions (31, 32).

The signature sequence from OTUs was then aligned to the database of known species (Ribosomal Database Project database) and annotated with ranges of taxonomic levels. As shown in Fig. 2A, Bacteroides, Faecalibacterium, and Prevotella were the dominant genera in all fecal samples. Total composition of bacterial community at the genus level in all participants is presented in Fig. S2A in the supplemental material. In addition, the relative abundance of phyla Bacteroidetes and Firmicutes and the corresponding Firmicutes/Bacteroidetes ratio were illustrated in Fig. S2B, indicating an increased abundance of Bacteroidetes along with decreased Firmicutes/Bacteroidetes ratio, consistent with previous research.

In order to identify specific bacterial taxa associated with GHPA, we compared fecal microbiota by means of Wilcoxon rank-sum test with false-discovery rate (FDR) correction and LEfSe analysis and found several differential taxa at genus level between GH group and HC group (age and BMI adjusted). Strikingly, compared to that in the control group, abundance of only two taxa, Oscillibacter (FDR-adjusted P = 0.003) and Enterobacter genus (FDR-adjusted P = 0.043), was substantially increased in GH group (Fig. 2B and C), while the remaining genera were all significantly decreased in GH group, and Blautia (FDR-adjusted P < 0.001) and Romboutsia genera (FDR-adjusted P < 0.001) ranked in the top 2 (Fig. 2B and D). The results of LEfSe analysis, sorting taxa according to results of Kruskal-Wallis test, confirmed the prominent roles of Oscillibacter (P < 0.001), Enterobacter (P = 0.007), Blautia (P < 0.001), and Romboutsia (P < 0.001), with linear discriminant analysis (LDA) scores of 3.79, 3.59, 4.44, and 3.91, respectively (Fig. 2E). According to existing reports, the increase of Oscillibacter is associated with inflammatory status as well as increased permeability of the gastrointestinal tract (33–37). Meanwhile, in addition to acting as an opportunistic pathogen (38), the expanded Enterobacter may reinforce the fermentation and utilization of various carbon sources to produce metabolites such as acetic acid (39–41). On the contrary, the general reduction of other genera, consistent with decreased biodiversity, might confront patients with susceptibility to opportunistic pathogens (42). The reduction of gastrointestinal microbial diversity has frequently been related to various diseases, including diabetes, dyslipidemia, obesity, and hypertension (32, 43, 44). Previous studies have proposed that decreased gastrointestinal microbial diversity is frequently associated with a bloom of commensals that may become detrimental (pathobionts) and contribute to systemic inflammation (42, 45). As growth hormone (GH) and insulin-like growth factor-1 (IGF-1) play important roles in regulating the body metabolism, patients with growth hormone-secreting pituitary adenoma (GHPA) present with multiple metabolic abnormalities, such as diabetes (20 to 56%) and hyperlipidemia (13 to 51%) (22–24). In the present study, our data showed that gut microbial diversity in patients with GHPA was significantly decreased, which was similar to that found in patients with diabetes and dyslipidemia, and we therefore suspect that decreased microbial diversity might intermediate the vulnerability of GHPA patients to metabolic disturbances.

To further illuminate the comprehensive landscape of taxonomic alteration at species level in GHPA, fecal samples of 9 patients with GHPA along with those of 10 participants from HC group were randomly selected and analyzed by metagenomics shotgun sequencing (age and BMI adjusted). The metagenome results data are presented in Table S3. Profiles of phylotype at species level were presented, and Bacteroides vulgatus, Alistipes putredinis, Bacteroides dorei, and Bacteroides ovatus were the top four species (Fig. 3A). LEfSe analysis, taxa stratified according to Kruskal-Wallis test, showed that Alistipes shahii (P = 0.041, LDA score = 3.59) and Odoribacter splanchnicus (P = 0.007, LDA score = 3.53), together with Prevotella stercorea (P = 0.002, LDA score = 3.41), were notably enriched in the GH group, while uncultured phage crAssphage (P = 0.002, LDA score = 3.80), Sutterella wadsworthensis (P = 0.022, LDA score = 3.78), and Sutterella sp. KLE1602 (P = 0.007, LDA score = 3.78) were strikingly enriched in the healthy control group (Fig. 3B). It is worth noting that species Enterobacter sp. DC1 (P = 0.041, LDA score = 2.44) and Enterobacter sp. 940 PEND (P = 0.014, LDA score = 2.28) that were derived from Enterobacter were both prominently enriched in the GH group (Fig. 3B). Results of rank-sum test demonstrated that certain species, including Odoribacter splanchnicus (P = 0.006) and Alistipes shahii (P = 0.043), were remarkably increased in the GH group, while species such as Carnobacterium sp. N15 MGS 207 (P = 0.028), Phascolarctobacterium sp. CAG 207 (P = 0.030), and Bacteroides sp. 3_1_19 (P = 0.013) were notably increased in the HC group (Fig. 3C). Since the numbers of participants enrolled in 16S and metagenomics sequencing were unparalleled, variations in taxonomic profiles were inevitable, and rank-sum test at genus level (Fig. S2D) with metagenomics has revealed deviation from 16S rRNA sequencing with respect to differential taxa between GH and HC groups, thus causing difficulties in alignment and interpretation of results. For example, genus Oscillibacter was enriched in the GH group in the 16S study while the abundance of species derived from this genus was decreased in the GH group with metagenomics, probably resulting from systemic bias due to unequal numbers of participants. Comparison of notably altered taxa at genus level as well as the abundance of Oscillibacter with metagenomics (Fig. S2C) confirmed the inconsistency in regard to relative abundance of genus Oscillibacter between 16S (Fig. 2C) and metagenomics sequencing (Fig. S2C). Yet, despite the potential shift in taxonomic construction, results of Enterobacter at both genus and species levels were compatible, further augmenting its significance in GH group.

Next, we set out to probe for the functional profiles and the differences in the metabolic potential of gut microbiota in GHPA and healthy controls. Together, a noteworthy proportion of the overall genes in both groups were involved in metabolism, genetic information processing, and environmental information processing (Fig. S3A), corroborated by depiction of relative abundance of comprehensive elements in each group at level 1 (Fig. S3B) and level 2 (Fig. S3C).

Kyoto Encyclopedia of Genes and Genomes (KEGG) orthologs (KOs) is a taxonomic system of proteins (enzymes), in which proteins with highly similar sequences and similar functions along the same pathway are grouped together (46). To depict the detailed landscape of differential functional modules between the GH group and the control group, we initiated further comparison of KOs between the two groups by means of rank-sum test as well as LEfSe analysis and further annotated the varied KOs to related pathways.

According to the KOs abundance table, the P value of KOs’ difference between the two groups was calculated by rank-sum test. KOs with significantly different expression between the experimental group and the control group were screened according to the P value of <0.05. As a result, putative ABC transport system permease protein (P = 0.002), alpha-glucosidase (P = 0.028), 23S rRNA pseudouridine1911/1915/1917 synthase (P = 0.028), and thiamine-phosphate pyrophosphorylase (P = 0.004) were proved to be the top 4 abundant KOs in HC group (Fig. 4A), while an uncharacterized protein (P = 0.017), DNA polymerase III subunit epsilon (P = 0.017), glucose-1-phosphate thymidylyltransferase (P = 0.035), and dTDP-4-dehydrorhamnose 3,5-epimerase (P = 0.043) were the 4 KOs with the highest abundance in the GH group (Fig. 4A).

As for LEfSe analysis, linear discriminant analysis (LDA) was used to estimate the effects on sample classification resulting from each KO. In the GH group, uncharacterized protein (P = 0.022, LDA score = 3.79), dTDP-4-dehydrorhamnose 3,5-epimerase (P = 0.027, LDA score = 3.15), uracil-DNA glycosylase (P = 0.011, LDA score = 3.03), and lysozyme (P = 0.022, LDA score = 2.87) were the most significantly enriched KOs, with RNA polymerase sigma-70 factor, the extracytoplasmic function (ECF) subfamily (P = 0.041, LDA score = 4.09), putative ABC transport system permease protein (P = 0.004, LDA score = 3.77), alpha-glucosidase (P = 0.014, LDA score = 3.38), and GDP-l-fucose synthase (P = 0.034, LDA score = 2.99) notably enriched in the counterpart (Fig. 4B). According to the relative abundance of each KO from the two groups, we further annotated the detected KOs to the existent pathways in KEGG database, acquiring the percentage of the overall KOs occupied by each pathway. Amino acid metabolism (17 in GH versus 8 in HC) together with global and overview maps (13 in GH versus 10 in HC) were the most remarkably enriched pathways in the GH group. In contrast, energy metabolism (4 in GH versus 14 in HC), membrane transport (2 in GH versus 12 in HC), and carbohydrate metabolism (13 in GH versus 18 in HC) were the most notably enriched pathways in the HC group (Fig. 4C).

Cooperatively, rank-sum test and LEfSe analysis illustrated the enrichment of KOs related to glycan metabolism (dTDP-4-dehydrorhamnose 3,5-epimerase and glucose-1-phosphate thymidylyltransferase), amino acid metabolism (putative endopeptidase, glutamate racemase, asparaginyl-tRNA synthetase, serine O-acetyltransferase, acetylornithine deacetylase, and phosphoribosylglycinamide formyltransferase 2), and utilization of certain carbohydrates (acetyl-CoA synthetase and galactokinase), in spite of the overall decline in KOs associated with carbohydrate metabolism. As reported by existed researches, amino acid metabolism and carbohydrate metabolism are involved in the production of short-chain fatty acids (SCFA) (47), while glycan metabolism is related to pathogen recognition and inflammation (48, 49). Thus, the functional profiles of gastrointestinal microbiota in GH patients may indicate the variance in SCFA metabolism and inflammatory status.

To capture the driving factors potentially causing the functional changes and respective contribution, FishTaco analysis was initiated (50). According to our results, Alistipes shahii was a major contributor with respect to amino acids metabolism, carbohydrate metabolism, and energy metabolism, while species of Enterobacter were the suppressive factors in amino acids metabolism, carbohydrate metabolism, and energy metabolism (Fig. S4). Though the FishTaco analysis outlined the prominence of GH-enriched taxa in respect to functional change, the results were not fully consonant with regard to the role of Enterobacter, defined as a suppressive factor in amino acids metabolism. Given the limitations of FishTaco analysis, further exploration into the exact role of this taxon is warranted.

In order to uncover the potential relationship between gastrointestinal microbiota and clinical characteristics of GH-secreting pituitary adenoma, we initiated spearman correlation test between differential taxa at the genus level and the species level and the clinical index, including Knosp grade, tumor size, granulation extent of adenoma, radiological remission, biochemical remission, preoperative (pre-op) nadir GH, change in nadir GH, preoperative (pre-op) fasting GH, change in fasting GH, preoperative (pre-op) IGF-1 index, preoperative (pre-op) IGF-1, preoperative (pre-op) FPG, change in IGF-1 index, ratio of change in IGF-1 index, ratio of change in fasting GH, and ratio of change in nadir GH.

At the genus level, notably, the Fusobacterium genus was positively correlated with both radiological remission (r = 0.46, P = 0.013) and biochemical remission (r = 0.39, P = 0.042), while Lachnospiraceae incertae sedis, another taxon enriched in the control group, was negatively correlated with radiological remission (r = −0.45, P = 0.016) and biochemical remission (r = −0.43, P = 0.022). Meanwhile, Oscillibacter, prominently enriched in the GH group as seen in Fig. 2B to E, was also negatively correlated with radiological remission (r = −0.42, P = 0.026) and biochemical remission (r = −0.40, P = 0.034) (Fig. 5A). Additionally, a significant positive correlation is detected between Lachnospiraceae incertae sedis and the densely granulated type of GH-secreting pituitary adenoma (r = 0.52, P = 0.005), which suggests the lower aggressiveness and higher response rate to pharmacological intervention compared to those of the sparsely granulated type (Fig. 5A). Furthermore, Enterobacter, which was enriched in GHPA as shown in Fig. 2B to E, was positively correlated with preoperation IGF-1 index (r = 0.39, P = 0.038), ratio of change in fasting GH (r = 0.41, P = 0.032), and ratio of change in nadir GH levels (r = 0.44, P = 0.018) (Fig. 5A).

At the species level, strikingly, Enterobacter sp. DC1 was positively correlated with ratio of change in fasting GH (r = 0.68, P = 0.043) and ratio of change in nadir GH levels (r = 0.72, P = 0.028), in line with the correlation analyses of Enterobacter. With regard to preoperation IGF-1 index, Enterobacter sp. DC1 was not significantly associated with this parameter in a positive manner (r = 0.58, P = 0.108), and this association was less remarkable than that observed between genus Enterobacter and preoperation IGF-1 index mentioned above (r = 0.39, P = 0.038). Collectively, correlation analyses of Enterobacter at genus and species levels revealed its potential fundamental role in GHPA. Meanwhile, species of Odoribacter, including Odoribacter splanchnicus (r = −0.82, P = 0.007), Odoribacter sp. UNK MGS 12 (r = −0.75, P = 0.020), and Odoribacter sp. N15 MGS 14 (r = −0.72, P = 0.030), were inversely correlated with ratio of change in nadir GH. However, Alistipes shahii, enriched in GHPA, was negatively related to the change in IGF-1 index (r = −0.72, P = 0.037) and ratio of change in IGF-1 index (r = −0.83, P = 0.008) (Fig. 5B).

Based on the results above, we then focused on the Enterobacter genus to investigate its possible influence on clinical parameters. According to the relative abundance of the Enterobacter genus, we divided patients into two groups: patients with low abundance of Enterobacter (low-abundance group, n = 14) and patients with high abundance of Enterobacter (high-abundance group, n = 14). As shown in Fig. 5C, compared with those in the low-abundance group, preoperative IGF-1 (P = 0.041) and IGF-1 index (P = 0.047) were significantly higher in patients with high abundance of Enterobacter. In addition, ratio of change in IGF-1 index (P = 0.035) and ratio of change in nadir GH (P = 0.022) were also higher in patients with high abundance of Enterobacter than in the low-abundance group.

To investigate the latent role of gut microbiome as the biomarker of GHPA, we constructed a random forest classifier model which could discriminate the patients from healthy controls with better specificity and sensitivity. We first conducted a 10-fold cross-validation on a random forest model between the patients (n = 15) and healthy controls (n = 51). As a result, 6 genera were selected as the optimal markers and were ranked according to the relative degree of contribution as Blautia, Dorea, Romboutsia, Anaerostipes, Saccharibacteria genera incertae sedis, and Oscillibacter in regard to mean decrease accuracy (Fig. 6A), while genera Blautia, Dorea, Anaerostipes, Romboutsia, Oscillibacter, and Fusicatenibacter ranked in the top 6 with respect to mean decrease Gini (Fig. 6B). Finally, we utilized the remaining samples (samples of 13 GHPA and 16 healthy controls) to validate the diagnostic efficacy of this model. Area under the curve (AUC) value was 0.981, indicating the potential of gut microbiome as a biomarker for GHPA (Fig. 6C). This machine-learning algorithm calculated mean decrease accuracy along with mean decrease Gini to determine the order in which taxa were selected and evaluate the relative importance of individual taxa in respect to classification (GH group or control group). Nevertheless, this algorithm reflected relative importance of individual taxa with respect to classification, and it did not generate the result that these taxa were the most significantly varied taxa between the two groups. Given the fact that genus Enterobacter was not ranked top among the most notably altered taxa in 16S rank-sum test and LEfSe analysis, it may explain the reason why the specific taxon was not included in the parameters for robust differentiation. Collectively, our results demonstrate a substantial link between gastrointestinal microbiome and vital parameters of GH patients, hinting at the prospective role of microbiome in GHPA.

Clinical indices mentioned in this study included demographic features (age, gender, and body mass index), comorbidities (hypertension [HBP] and diabetes mellitus [DM]; diagnoses were based on 2017 American Heart Association [73] and 1999 World Health Organization criteria [74], respectively), tumor features (including the recurrence of tumor, tumor size, Knosp grading of tumor, a grading system to evaluate invasiveness of pituitary adenoma as well as granulation extent of tumor, a pathological index associated with aggressiveness and resistance to pharmacological therapy [75, 76]), laboratory tests (fasting plasma glucose [FPG], fasting GH, nadir GH, IGF-1, and IGF-1 index), and endpoints of management (biochemical remission was defined as an age-normalized serum IGF-1 value and a random GH of <1.0 ng/mL according to existing consensus (30), and we defined radiological remission as no recurrence or residual tumor on follow-up magnetic resonance imaging).

GH group and HC group are abbreviations of growth hormone-secreting pituitary adenoma group and healthy controls group, respectively. The recurrence of tumor denoted the status that the patient carried the recurrent tumor despite past treatment. Body mass index (BMI) was calculated by weight/(height)², where weight is in kilograms and height is in meters. Tumor size was evaluated via the formula 1/2 × length × width × height, where all measurements are in millimeters. IGF-1 index was the ratio of the measured IGF-1 value to the upper limit of normal. Fasting GH refers to the GH level at the 0 point during oral glucose tolerance test (OGTT), while nadir GH is the nadir GH level during OGTT.

Clinical parameters were calculated as follows. Change in fasting GH = preoperation fasting GH – postoperation fasting GH. Change in nadir GH = preoperation nadir GH – postoperation nadir GH. Change in IGF-1 index = preoperation IGF-1 index – postoperation IGF-1 index. Ratio of change in fasting GH = (preoperation fasting GH – postoperation fasting GH)/preoperation fasting GH. Ratio of change in nadir GH = (preoperation nadir GH – postoperation nadir GH)/preoperation nadir GH. Ratio of change in IGF-1 index = (preoperation IGF-1 index – postoperation IGF-1 index)/preoperation IGF-1 index. Postoperative values were obtained 6 months after surgery.

This study was approved by the local ethics committee of the Huashan Hospital (KY2015-256), Fudan University, and all participants gave informed consent. Patients with newly diagnosed GH-secreting pituitary adenoma (GHPA) were recruited from the Department of Neurosurgery from March 2018 to April 2019. Diagnostic criteria of GHPA were based on the endocrine society clinical practice guideline (30). To exclude the potential effect on microbiome analysis from potential confounding factors, exclusion criteria were (i) antibiotics and/or probiotics usage within 6 months, (ii) a history of gastrointestinal tumors and/or inflammatory diseases, and (iii) immunosuppressive drug usage within 6 months.

Clinical parameters were collected during the perioperative period and follow-up hospitalization on the sixth month after surgery. Clinical information included demographic features (age, gender, and body mass index), comorbidities (hypertension and diabetes mellitus), tumor features (recurrence, tumor size, Knosp grade, and dense/sparse granulation), laboratory tests (fasting plasma glucose, fasting and nadir growth hormone level during oral glucose tolerance test, insulin-like growth factor-1, and IGF-1 index), and endpoints of management (biochemical remission as well as radiological remission). Further, preoperation BMI, FPG, fasting and nadir GH, IGF-1, IGF-1 index, and their counterparts were compared 6 months after surgery.

Healthy controls with general demographics matched were selected from a large cohort established by the Department of General Surgery, Shanghai Tenth People’s Hospital, School of Medicine, Tongji University (77). The exclusion criteria mentioned above were also applied to the control group.

A total of 28 patients with GH-secreting adenoma and 67 healthy controls were included in this study, and their stool samples were submitted to 16S rRNA amplicon sequencing. Meanwhile, stool samples of 9 patients and 10 healthy controls were randomly selected and further analyzed by metagenomics sequencing. Detailed composition, biodiversity, taxonomic difference, and potential functional implication of the gut microbiota were delineated. Furthermore, the correlation between gut microbiome and clinical characteristics was analyzed. A classifier to discriminate GH-secreting adenoma from the general population was also constructed and validated.

GH and IGF-1 were measured as described previously (78). Briefly, GH was measured by a two-site chemiluminescent immunometric assay (AutoDELFIA hGH, PerkinElmer Life and Analytical Sciences, Wallac Oy), with intraassay coefficient of variation (CV) of 5.3 to 6.5%, interassay CV of 5.7 to 6.2%, and sensitivity up to 0.01 μg/L (0.026 mU/L).

IGF-1 was measured with the Immulite 2000 solid-phase, enzyme-labeled chemiluminescent immunometric assay (Siemens Healthcare Diagnostic Products Limited, UK); normal age-appropriate ranges are as follows:1 to 6 years, 49 to 327 μg/L; 7 to 11 years, 57 to 551 μg/L; 12 to 13 years, 143 to 850 μg/L; 14 to 16 years: 220 to 996 μg/L; 17 to 18 years, 163 to 731 μg/L; 19 to 20 years, 127 to 483 μg/L; 21 to 35 years, 115 to 358 μg/L; 36 to 50 years, 94 to 284 μg/L; >50 years, 55 to 238 μg/L; intraassay CV, 2.3 to 3.5%; interassay CV, 7.0 to 7.1%; sensitivity, 20 μg/L. IGF-1 index = IGF-1/upper limit of normal range (ULN) (79).

Plasma glucose was measured by HITACHI 7600 Biochemical Analyzer (Tokyo, Japan).

Fecal samples of patients were collected preoperatively via Sarstedt feces collection containers (SARSTEDT, Nümbrecht, Germany) following the manufacturer’s instructions at the hospital, immediately frozen, transferred to our lab on ice, and immediately stored at −80°C until sample processing. Fecal samples were collected by the same staff. As for the healthy controls group, fecal samples were processed and completed during community gastrointestinal disorders screening according to the identical protocol as mentioned above.

DNA was extracted from each fecal sample using improved protocol based on the QIAamp Fast DNA stool minikit (Qiagen, Germany). In detail, 1 mL of InhibitEX buffer and a proper amount of glass beads (0.5 mm diameter, Qiagen) were added to each 200-mg sample of feces. The mixture was homogenized and beat with 60 Hz for 1 min twice with a homogeneous instrument (FASTPREP-24, Aosheng Biotech, China). Afterwards, the DNA purification was performed according to the manufacturer’s instructions.

The V3-V4 region of the bacteria 16S rRNA genes was amplified by PCR (95°C for 3 min, followed by 30 cycles at 98°C for 20s, 58°C for 15s, and 72°C for 20s and a final extension at 72°C for 5 min) using barcoded primers 341F 5′-CCTACGGGRSGCAGCAG-3′ and 806R 5′-GGACTACVVGGGTATCTAATC-3′. PCRs were performed in a 30 μL mixture containing 15 μL of 2× KAPA Library Amplification ReadyMix, 1 μL of each primer (10 μM), 50 ng of template DNA, and ddH₂O. Negative controls consisting of empty sterile storage tubes were processed for DNA extraction and amplification using the same procedures and reagents used for the 95 samples. There was no detectable amplification in the negative controls. Amplicons were extracted from 2% agarose gels and purified using the AxyPrep DNA gel extraction kit (Axygen Biosciences, Union City, CA, USA) according to the manufacturer’s instructions and quantified using Qubit 2.0 (Invitrogen, USA). All quantified amplicons were pooled to equalize concentrations for sequencing using Illumina MiSeq/HiSeq (Illumina, Inc., CA, USA). The paired-end reads of 250 bp were overlapped on their 3 ends for concatenation into original longer tags by using PANDAseq (https://github.com/neufeld/pandaseq↗, version 2.9) (80). DNA extraction, library construction, and sequencing were conducted at Realbio Genomics Institute (Shanghai, China).

Assembled tags, trimmed of barcodes and primers, were further checked on their rest lengths and average base quality. 16S tags were restricted between 220 bp and 500 bp such that the average Phred score of bases was no worse than 20 (Q20) and no more than 3 ambiguous bases (N). The copy number of tags was enumerated and redundancy of repeated tags was removed. Only the tags with frequency more than 1, which tend to be more reliable, were clustered into OTUs, each of which had a representative tag. Operational taxonomic units (OTUs) were clustered with 97% similarity using UPARSE (http://drive5.com/uparse/↗) (81), and chimeric sequences were identified and removed using Usearch (version 7.0.1090) (82). Each representative tag was assigned to a taxon by RDP Classifier against the RDP database (http://rdp.cme.msu.edu/↗) (83) using confidence threshold of 0.8. OTU profiling table, and alpha diversity analyses were also achieved by python scripts of QIIME (version 1.9.1) (84).

Weighted UniFrac uses phylogenetic information to compare microbial community differences between samples, taking the abundance of sequences into account. The calculated results, taking into account the evolutionary distance between taxa, are used as an index to measure beta diversity. The larger the index, the greater the difference between samples. By clustering the results of UniFrac, samples with smaller UniFrac distance are clustered together, reflecting the similarity among samples. In the current study, we exploited QIIME to perform weighted UniFrac analysis, with detailed codes provided in github repository (https://github.com/xxhym/bioinformatic_codes↗), and original script can be found on the website http://qiime.org/scripts/beta_diversity.html↗. The results were presented via principal coordinates analysis (PCoA) and nonmetric multidimensional scaling (NMDS) analysis. To depict the variance between different groups, Adonis test together with principal coordinates analysis (PCoA) were exploited. In addition, multiresponse permutation procedure (MRPP), a mathematic form to present difference in weighted UniFrac (85, 86), was utilized to collaborate with NMDS analysis to better illustrate the intergroup difference.

To acquire differential taxa, LEfSe analysis and intergroup rank-sum test analysis were conducted. LEfSe was conducted by means of nonparametric Kruskal-Wallis test, and it used linear discriminant analysis (LDA) to estimate the influence of the abundance of each component (taxon) on the difference effect, thus illustrating the communities or taxa with significant difference influence on sample division (87). The detailed code is provided in our github repository (https://github.com/xxhym/bioinformatic_codes↗), and the analysis can also be achieved either online (http://huttenhower.sph.harvard.edu/galaxy/root/index↗) or via local software (https://toolshed.g2.bx.psu.edu/↗). On the other hand, the rank-sum test was used to analyze the significant differences between different groups to identify the species that had significant effects on the division between groups. In this analysis, the Wilcox test function of the stats package in R was used for the difference analysis between two groups, and the Kruskal test function of the stats package in R was used for the difference analysis among more than two groups.

Thereafter, using phylotype at genus level and clinical indices of GH patients, the correlations between the two components are calculated to unfold important internal relationships, graphically shown through a heatmap generated through corrplot package in R.

DNA was extracted from each fecal sample using improved protocol based on the manual of QIAamp Fast DNA stool minikit (Qiagen, Germany). In detail, 1 mL of InhibitEX buffer and a proper amount of glass beads (0.5 mm diameter, Qiagen) were added to each 200-mg sample of feces. The mixture was homogenized and beat with 60 Hz for 1 min twice with a homogeneous instrument (FASTPREP-24, Aosheng Biotech, China). Afterwards, the DNA purification was performed according to the manufacturer’s instructions, the DNA concentration was measured with a NanoDrop (Thermo Scientific) and Qubit 2.0 (Invitrogen, USA), and the molecular size was estimated by agarose gel electrophoresis.

Following the Illumina TruSeq DNA Sample Prep v2 guide (Illumina, Inc., San Diego, CA, USA), we constructed the DNA libraries with approximately 500 bp insert sizes for each sample. The quality of all libraries was evaluated using an Agilent 2100 bioanalyzer (Agilent Technologies, Wokingham, UK) and the Agilent 2100 DNA 1000 kit. All samples were subject to 150 bp paired-end sequencing on an Hiseq X-10 platform (Illumina, Inc., San Diego, CA, USA).

Illumina raw reads were screened according to the following criteria: (i) adaptor contamination reads were removed, (ii) reads containing more than three ambiguous N bases were removed, (iii) reads containing low-quality (Q < 20) bases were trimmed, and (iv) reads containing less than 60% of high-quality bases (Phred score ≥ 20) were deleted. Then, clean reads were subjected to bacterial genomes from the National Center for Biotechnology Information GenBank with SOAPaligner (version 2.21) (88) and reads that mapped to the host genome were abandoned. The subsequent reads were selected for further analysis.

Clean reads were aligned to the NCBI GenBank database (https://ncbi.nlm.nih.gov/genbank/↗) for the detection of known bacteria, fungi, viruses, and archaea by SOAPaligner 2.21. Then, the aligned reads were classified as kingdom, phylum, class, order, family, genus, and species to count classification and abundance, and taxonomic relative abundance profile at different levels was generated.

Preprocessed reads were assembled by SOAPdenovo (version 1.05) (89). During the assembly process, the assembled scaffolds were split into contigs at ambiguous Ns and contigs no less than 500 bp were kept for further analysis. The k-mers with maximum N₅₀ were selected as the final assembly results. Software MetaGeneMark (http://exon.gatech.edu/GeneMark/meta_gmhmmp.cgi↗) (90) was used to predict open reading frames (ORFs) in the assembled scaffolds, and predicted ORFs less than 100 bp were removed.

The nonredundant gene catalog set was constructed by pairwise comparison of predicted ORFs (gene length > 100 bp) using CD-HIT (version 4.5.7) (91). Two sequences with coverage of ≥90% and identity of ≥95% were considered redundant, and the longer one was regarded as the representative. The final nonredundant gene catalog contained 885,596 ORFs with an average length of 818.90 bp.

Clean reads were aligned to the genes in the nonredundant catalogue by using SOAPaligner. The calculation formula of gene abundance used was from the study by Qin et al. (92).

Using BLAST (version 2.2.28+) (93), the nonredundant genes were annotated against the KEGG (Kyoto Encyclopedia of Genes and Genomes) database. When assembled protein sequence was similar (score of ≥ 60 and E value of <1e–5) to a protein sequence in the database, the assembled protein was considered to play the same role as the protein in the database. We accumulated the relative abundance of all orthologous genes to generate the relative abundance of each KO.

According to the abundance of samples, the Bray-Curtis distance between samples was calculated and the hierarchical clustering was carried out to analyze the distance relationship between samples and whether there were outliers.

Afterwards, based on the abundance of species, the P value of each taxon at species level was calculated by Wilcox rank-sum test. Then, the species with significant difference between the two groups were screened out according to the P value and exhibited in the form of boxplots, with the general threshold of a P value of ≤0.05. Meanwhile, LEfSe analysis was initiated to probe microbial biomarkers and reveal genomic characteristics. Similar to LEfSe analysis with 16S rRNA, we used the nonparametric factorial Kruskal-Wallis (KW) sum-rank test to detect taxa with significantly different abundance between the two groups, and thereafter linear discriminant analysis (LDA) was used to estimate relative influence resulting from abundance of each component (species) on the total intragroup difference (87), thus locating species fundamental in respect to sample division.

Then, we combined the abundance of detected genes and the corresponding annotation from KEGG database to describe relative abundance of comprehensive functional elements at level 1 and level 2 of each sample in the form of barplots.

Thereafter, integrating taxa at species level and clinical parameters of GH patients, the spearman correlations between the two components were calculated to unveil vital internal relationships, presented as a heatmap generated through corrplot package in R.

Last, but not least, the FishTaco analysis (50), based on permutation-based approach to link species to gene function/metabolic pathways, was carried out in participants with metagenomics. Therefore, driving factors responsible for the functional changes and their corresponding contributions to the alterations were attained, probing into the possible mechanism of changes in gene function/metabolism pathway caused by species.

Statistical analysis of the clinical metadata was conducted in SPSS. Continuous data were described as mean ± standard deviation (interquartile range for nonnormalized data), and the count data were described as count and proportion. Paired t tests were used to assess the differences between normally distributed continuous variables, and Mann-Whitney U tests were used with variables that failed the normality test. The false-discovery rate (FDR) algorithm was used for multiple comparisons to control the chance of generating unwanted false positives.

All bioinformatic codes are available on Github repository (https://github.com/xxhym/bioinformatic_codes↗). 16S rRNA amplicon sequencing resources are available on BioProject ID PRJNA743650, and metagenomics sequencing resources are available on BioProject ID PRJNA759070.

All bioinformatic codes are available on Github repository (https://github.com/xxhym/bioinformatic_codes↗). 16S rRNA amplicon sequencing resources are available on BioProject ID PRJNA743650, and metagenomics sequencing resources are available on BioProject ID PRJNA759070.