Approval was obtained from the ethics committee of Shenzhen Second People’s Hospital [No. 20210617001]. The procedures used in this study adhere to the tenets of the Declaration of Helsinki. Informed consent was obtained from all individual participants included in the study. This was a cross-sectional study, and all participants were recruited from Shenzhen Second People’s Hospital from April to November 2020. After excluding those with complications of the endocrine system (such as diabetes, hyperthyroidism, acromegaly, etc.) and a family history of endocrine neoplasia, 44 subjects who agreed to provide faecal samples were recruited for our study. Written informed consent was obtained from all participants. In total, 29 diagnosed patients were categorized into the PA group, including 13 IPA patients and 16 NIPA patients, while 15 healthy subjects were categorized into the normal control (HC) group. All samples were confirmed by pathological examination. The diagnostic criteria for IPA included Knosp classification, intraoperative findings of tumour invasion, and Ki-67 labelling index according to previous studies [2, 14]. The demographic and clinical characteristics were derived from the electronic medical records of the hospital information system, including age, gender and body mass index, etc. Stool samples were collected from the first defecation (almost in the morning) using sterile swabs about 5 g within 3 days. The samples were snap-frozen in dry ice right after collection, and transferred to −80°C freezers in 30 min and stored with sterile tubes until subsequent experiments. To avoid the patient identification, we gave each study subject and sample with a unique subject ID and sample ID, respectively.

Shotgun metagenome sequencing was performed based on the method of a previous study [15]. Briefly, total faecal genomic DNA was extracted using the Stool Genomic DNA Kit (CW2092S; CWBIO) according to the manufacturer’s instructions.

DNA concentrations were measured using a NanoDrop 2000 system (Thermo Fisher Scientific), and the DNA molecular size was estimated by agarose gel electrophoresis. The DNA library was constructed using the TruSeq DNA Sample Preparation Kit (Illumina, San Diego, CA, United States). Libraries were sequenced using the 150 base pair paired-end strategy on the Illumina X-ten platform at Novogene Bioinformatics Technology Co., Ltd. A total of 274.5 GB original sequencing data were obtained from 44 samples (6.2G per sample, ranged from 4.8 to 8.7 G) and deposited in the Sequence Read Archive (PRJNA799832), where 6.1 ± 0.41 G for IPA group, 6.5 ± 0.66 G for NIPA group, 6.0 ± 1.13 G for HC group.

Sequence adapter contamination and low-quality reads were discarded from the raw sequencing reads using Fastp software [16], and the remaining reads were filtered to eliminate human host DNA based on the human genome reference (hg19) using BWA software [17]. The relative abundance profiles of microbes were calculated using MetaPhlAn3 software with default parameters [18], and the abundance of microbial metabolic pathways against the MetaCyc database was calculated using HUMAnN3 software [18].

Alpha diversity was estimated on the basis of the microbial profile of each sample using Vegan’s diversity function in the R package. Partial least squares discriminant analysis (PLS-DA) was used to reveal taxonomic changes among the three groups, and variable importance in projection (VIP) scores were used to rank the ability of different taxa to discriminate between different groups. LEfSe (Linear Discriminate Analysis of Effect size) software was used to estimate the differentially abundant taxa between two groups with p < 0.05 [19]. The significantly different functional metabolic pathways were investigated using DESeq2 software [20]. The microbial ecological network was constructed using SpiecEasi software at the genus level and visualized with Gephi software [21].

Differences in characteristics were assessed with t tests (or analysis of variance, ANOVA) and chi-squared tests for continuous and categorical variables, respectively. All these analyses were performed using R software (v 3.6.3) at 0.05 significance levels.

A classification model was built to classify subjects with PA based on microbial profiles using the support vector machine (SVM) algorithm with 10-fold cross validation. The AUC (area under the curve) was calculated to evaluate the model performance.

A total number of 44 participants were enrolled in this study. We only collected the subjects without complications of the endocrine system and a family history of endocrine neoplasia. According to the severity, the PA patients were further divided into IPA (n = 13) and NIPA (n = 16) subgroups. On the other hand, 15 healthy subjects were enrolled as the control group (HC). The demographic and clinical characteristics of all the participants are summarized in Table 1. They are all Han Chinese and live in Guangdong Province. There was no significant difference among the three groups in terms of age, BMI (body mass index) or sex. In addition, adrenocorticotropic hormone (ACTH), prolactin (PRL), thyroid stimulating hormone (TSH), luteinizing hormone (LH), follicle stimulating hormone (FSH), and growth hormone (GH) presented without significant differences between IPA and NIPA groups.

To evaluate the overall characteristics of the gut microbiota in the HC, NIPA and IPA groups, we compared the α-diversity estimators (Shannon index). No differences existed in the mean values of the Shannon index between IPA (2.40 ± 0.82 versus 2.29 ± 0.73, p = 0.73) and NIPA (2.65 ± 0.44 versus 2.29 ± 0.73, p = 0.11) groups, compared to the HC group.

We also performed partial least squares discriminant analysis (PLS-DA) to compare the gut microbial profiles. The results showed that there was a distinct clustering pattern between samples from subjects with pituitary adenoma and healthy controls (Figure 1A), thus reflecting their separations in microbial composition. The variable importance in projection (VIP) score for the gut microbiome showed that Oscillibacter sp. 57_20, and Clostridium innocuum was the top 2 microorganisms contributing to the group separation (Figure 1B).

The average relative abundance of the dominant phyla shown in Figure 2A included Actinobacteria, Bacteroidetes, Firmicutes and Proteobacteria, where Bacteroidetes and Firmicutes accounted for more than 80%. Notably, the ratio of average relative abundance of Firmicutes/Bacteroidetes increased from HC group (0.44) to IPA (0.75) and NIPA (0.71) groups. Meanwhile, the predominant species were Alistipes putredinis, Anaerostipes hadrus, Bacteroides caccae, Bacteroides dorei, Bacteroides ovatus, Bacteroides plebeius, Bacteroides stercoris, Bacteroides thetaiotaomicron, Bacteroides uniformis, Bacteroides vulgatus, Bifidobacterium pseudocatenulatum, Collinsella aerofaciens, Escherichia coli, Eubacterium rectale, Eubacterium sp. CAG_38, Faecalibacterium prausnitzii, Fusicatenibacter saccharivorans, Klebsiella pneumoniae, Lachnospira pectinoschiza, Parabacteroides distasonis, Prevotella copri, Roseburia faecis, Roseburia inulinivorans, Ruminococcus bromii and Ruminococcus gnavus (Figure 2B). Compared to the HC group, the mean relative abundances of Anaerostipes hadrus (1.85% versus 0.21%, 1.39% versus 0.21%), Bacteroides dorei (2.02% versus 0.68%, 4.41% versus 0.68%), Bacteroides stercoris (2.08% versus 1.45%, 5.65% versus 1.45%), Collinsella aerofaciens (1.48% versus 0.37%, 1.27% versus 0.37%), Eubacterium rectale (3.29% versus 2.19%, 3.81% versus 2.19%), and Fusicatenibacter saccharivorans (2.62% versus 0.36%, 1.22% versus 0.36%) increased, whereas Bacteroides vulgatus (6.62% versus 14.24%, 7.91% versus 14.24%), Faecalibacterium prausnitzii (4.24% versus 5.80%, 2.65% versus 5.80%), and Prevotella copri (10.43% versus 19.36%, 9.95% versus 19.36%) decreased in both the NIPA and IPA groups.

Furthermore, we performed LEfSe analysis and revealed significantly different bacteria among the three groups. At the species level, specifically, the relative abundance of Clostridium innocuum was higher, while Oscillibacter sp. 57_20 and Fusobacterium mortiferum were lower in the NIPA group than in the HC group (Figure 2C). Likewise, enrichment of Clostridium innocuum, along with decreases of Oscillibacter sp. 57_20 and Fusobacterium mortiferum were also observed in the IPA group, compared to the HC group (Figure 2D). This result was consistent with the PLS-DA results.

To explore the efficacy of pituitary adenoma classification based on the gut microbial profile, we used the SVM algorithm with 10-fold cross-validation to build a classification model at the species level. We identified 10 species that could be used to predict the occurrence of PA, including Oscillibacter sp. 57_20, Fusobacterium mortiferum, and Clostridium innocuum. The receiver operating characteristic curve of the model was shown in Figure 3A with an AUC of 0.88.

The metabolic pathways of the PA and HC groups were predicted using HUMAnN3 software against the MetaCyc database. The pathways that were significantly enriched in the PA group were “formaldehyde oxidation I,” “superpathway of haem b biosynthesis from glycine,” “Bifidobacterium shunt,” “fatty acid beta-oxidation I,” “ketogluconate metabolism,” “methanogenesis from acetate,” “superpathway of pyrimidine ribonucleoside degradation,” “ppGpp metabolism,” “trehalose degradation V,” “formaldehyde assimilation II,” “glycogen degradation I,” “superpathway of L-lysine, L-threonine and L-methionine biosynthesis I,” “lactose and galactose degradation I,” “superpathway of glycerol degradation to 1,3-propanediol,” “sucrose degradation IV,” “aspartate superpathway,” “C4 photosynthetic carbon assimilation cycle,” “L-methionine biosynthesis I,” “pyrimidine deoxyribonucleotides de novo biosynthesis IV,” “gamma glutamyl cycle,” “superpathway of L-homoserine and L-methionine biosynthesis,” “superpathway of S-adenosyl-L-methionine biosynthesis,” and “guanosine nucleotides degradation III.” However, the microbial functions related to “dDTP-beta-L-rhamnose biosynthesis,” “6-hydroxymethyl-dihydropterin diphosphate biosynthesis I and III,” and “nitrate reduction VI” were higher in the gut microbiome of the HC group (Figure 4).