Bacterial microbiome profiles of the inflamed terminal ileum mucosa in active Crohn’s disease patients

Recruitment of participants

Twelve CD patients with ileal and colon inflammation in the active phase (CD activity index, CDAI > 150) who were not undergoing treatment were enrolled in this study. Samples of patients with CD were collected from January to November 2019 at the Affiliated Suzhou Hospital of Nanjing Medical University, Suzhou Municipal Hospital. The characteristics of enrolled patients are shown in Table 1. Patients who accepted therapeutic antibodies and antibiotic drug treatment for six months were excluded. Ethical approval for this study was obtained from the Ethics Committee of Nanjing Medical University. All experiments were performed in accordance with relevant guidelines and regulations by the Research Ethics Boards at the Ethics Committee of Suzhou Municipal Hospital, Nanjing Medical University (NO. K2021-064-K01). Written informed consent was obtained from all participants.

Collection of terminal ileum mucosal samples

The samples were collected as part of the routine diagnostic workflow. Only a subset of these samples was utilized for research purposes. Two inflamed mucosal biopsies were collected from the terminal ileum of CD patients. Two noninflamed mucosal biopsies located 10 cm from the inflamed mucosa of the terminal ileum were collected. Biopsy specimens were immediately immersed either in liquid nitrogen for rapid freezing (one inflamed and one noninflamed mucosal biopsy from each patient with CD) or in 4% formaldehyde (one inflamed and one noninflamed mucosal biopsy from each patient with CD) as they underwent the colonoscopy. The specimens that were frozen in liquid nitrogen were stored at −80 °C before DNA extraction and library construction. The biopsies fixed in 4% formaldehyde were subjected to HE staining. The degree of tissue inflammation was identified by HE staining.

DNA extraction, library construction, and 16S rRNA gene sequencing

A total of 24 mucosal biopsies were collected from CD patients. Genomic DNA was extracted with a QIAamp^R DNA mini kit (No. 51306; Qiagen, Hilden, Germany) according to the manufacturer’s instructions.

A total of 30 ng of qualified genomic DNA samples and their corresponding fusion primers were selected for PCR, and parameters were set for amplification. Furthermore, Agencourt AMPure XP magnetic beads were used to purify the PCR amplification products and subsequently dissolve them in an elution buffer, which was labeled to complete the database construction. An Agilent 2100 BioAnalyzer was used to detect the range and concentration of fragments in the library. The HiSeq platform was selected for sequencing of qualified libraries according to the size of the inserted fragments.

Data analysis of microbial load

OTU clustering and analysis

Tags were clustered into operational taxonomic units (OTUs) by USEARCH (v7 .0.1090) to compare with the gold database (v20110519) for chimeric removal. The Usearch global method was used to calculate the abundance of OTUs. The following analyses were conducted to characterize OTUs, including an OTU Venn diagram, principal components analysis (PCA), species accumulation curves, partial least squares discrimination analysis, PLS-DA, and OTU rank curve.

Species annotation, composition, and biodiversity analysis

Species annotation was carried out by comparing the OTU representative sequence with the database through RDP Classifier (V2.2) software. The reliability threshold was set to 0.8.

To obtain the species classification information corresponding to each OTU, the RDP classifier Bayesian algorithm was used to classify the OTU representative sequences. The community composition of each sample was counted at all levels including phylum, class, order, family, genus, and species.

The species diversity of a single sample was analyzed by alpha diversity. Alpha diversity analysis was performed using MOTHUR (version 1.31.2). The observed species index, chao index, ace index, Shannon index, simpson index, and good-coverage index are shown to describe diversity. The first four indices previously listed are directly correlated with species abundance, such that a greater value indicates a greater abundance of species in the sample. However, the fifth index listed (good-coverage index) has the opposite relationship, so a smaller value indicates more abundant species in the sample. In other words, as the value of the good-coverage index increases, the probability that the sequence in the sample will not be measured decreases. The difference in alpha diversity between groups is displayed in the alpha diversity box chart. The Wilcoxon test was used for comparisons between the two groups.

Beta diversity analysis was used to compare the degree of difference in species diversity between pairs of samples. Beta diversity was performed through QIIME (v1.80). The heatmap of beta diversity is shown to reflect the similarity of bacterial species between the samples. A box chart was generated by the ggplot package of R (v3.4.1) to display differences in beta diversity between groups.

Analysis of microbiome phylogeny

A phylogenetic branch tree was built at five levels including phyla, class, order, family, and genus (the top 50 most abundant were selected) using FastTree (version 2.1.3). The phylogenetic tree is ordered from the inside to the outside, and each circle represents a level, followed by the phylum, class, order, family, and genus levels. Individual phyla are differentiated by color, and node size is representative of species abundance. The outer ring is the abundance heatmap, each ring is a set of samples (a sample), each set of samples (each sample) corresponds to a color, and the color depth varies with species abundance. The length of the branch indicates the difference in evolutionary distance. The closer the evolutionary relationship is, the closer the evolutionary tree species is.

Differential bacterial composition analysis between groups

Linear discriminant analysis effect size (LEfSe) software was used to discover microbes that were significantly different (https://huttenhower.sph.harvard.edu/galaxy/) (Segata et al., 2011). LEfSe was applied to construct a cluster diagram to determine the differential bacterial composition between groups. Taxonomic analysis at the order level demonstrated that all significant associations represented entire orders, with no significant signals coming from unclassified members within these taxa. An LDA diagram is shown to display microbes with significant differences. The Wilcoxon rank-sum test was conducted to calculate the significant difference between groups.

Enterotype analysis

The Jensen–Shannon distance and PAM (partitioning around medoids) were calculated according to the relative abundance at the genus level for clustering. The optimal clustering K value was calculated by the Calinski–Harabasz (CH) index. Then, between-class analysis (BCA, K ≥ 3) was employed for visualization. The software R (v3.4.1) cluster and clusterSim packages were used in this analysis. Different samples with similar dominant flora structures can be grouped in this analysis, yet this is only suitable under certain conditions; these conditions include bacterial group typing of specific environmental samples, such as enterotypes, vaginal typing (cervicotype), and oral typing.

Function prediction and correlation analysis

Kyoto Encyclopedia of Genes and Genomes (KEGG) functional prediction of bacteria with upregulated abundance in inflamed mucosa was carried out by PICRUSt2. The Wilcoxon test was used to discover functional differences between the groups.

The different bacteria at the genus level were obtained by the rank sum test, and a Spearman correlation heatmap among dominant species was generated by R software. The map reveals important patterns and relationships among dominant species.

Characteristics of participants

The clinical characteristics of the participants in this study are shown in Table 1. Twelve CD patients ranging from 27 to 53 years old, including seven women (mean age 35.6) and five men (mean age 44) were enrolled.

OTU clustering and analysis

A total of 1,158 OTUs were shared between inflamed tissues and noninflammatory mucosa tissues (self-controls) of CD patients. Compared to self-control, 753 OTUs were specific to inflamed tissues (Fig. 1A). OTU Core-Pan displays common and unique OTUs for all samples in petal plots. Three OTUs were common in all tissues (Fig. 1B). Species accumulation curves showed an upward trend at the end of the curve tended to be flat, which indicated that the sampling amount was sufficient (Fig. 1C). Partial least squares discrimination analysis (PLS-DA) (Fig. 1D) showed the results could be divided into two distinct groups between inflamed tissues and noninflamed tissues of CD patients based on differences in the abundance of OTUs. The abscissa was sorted by sample OTU abundance (from high to low), and the ordinate is OTU abundance. The abundance of species in the sample is reflected by the length of the horizontal axis of the curve. As the curve widens, the abundance in species composition of the sample increases. The uniformity of the species in the sample is reflected by the shape of the vertical axis of the curve. If the curve is flatter, the uniformity of the species composition in the sample is greater (Fig. 1E).

Bacterial microbiota composition and correlations in inflamed mucosa

Histograms at the genus level showed that the main genus in both inflamed tissues and noninflamed tissues (control) of CD patients, included Escherichia (0.2771, 0.2904) (relative abundance of the genus in inflamed tissues and noninflamed tissues), Bacteroides (0.1212, 0.1571), Prevotella (0.0734, 0.0688), Fusobacterium (0.0663, 0.0847), Acinetobacter (0.0390, 0.0323), Pseudomonas (0.0294, 0.0315), Ruminococcus (0.0204, 0.0283), Megamonas (0.0184, 0.0241), Sphingomonas (0.0176, 0.0217) and Dorea (0.0162, 0.0350), which collectively composed the predominant bacterial microflora (Fig. 2A). Those species whose abundances were less than 0.5% of all samples were classified into others. Enterotype analysis showed that Enterobacteriaceae (Escherichia), Moraxellaceae (Acinetobacter), and Pseudomonadaceae (Pseudomonas) were the three principal components in the mucosa of CD patients (Fig. 2B). The abscissa represents principal component one, and the ordinate represents principal component two, which are the two principal components with the largest variance contribution rate.

According to the species composition graph, the intestinal mucosa of CD patients was rich in Bacteroidaceae (Bacteroides), Prevotellaceae (Prevotella), Fusobacteriaceae (Fusobacterium), and Enterobacteriaceae (Escherichia) (Fig. 2C). A phylogeny tree (Fig. 2D) shows that Actinobacteria, Bacteroidetes, Firmicutes, Proteobacteria, Verrucomicrobia, Nitrospirae, and Acidobacteria were the dominant phyla in the intestinal mucosa of CD patients. The biodiversity of bacterial microbiota in inflamed mucosa was not significantly different compared with noninflamed mucosa.

Alpha diversity analysis demonstrated that the biodiversity of the bacterial microbiota in the inflamed mucosa was increased, but it was not significantly different from that in the noninflamed mucosa (Fig. 3A). The mean, SD, and P-value of the observed species, as well as the Chao, ACE, Shannon’s diversity, Simpson’s diversity, and Good’s coverage, are shown in Table 2.

The beta diversity index box plot (Fig. 3B) and beta diversity matrix heatmap (Fig. 3C) showed that there were no significant differences between the noninflamed mucosa and inflamed mucosa of patients with CD.

Differential bacterial composition in inflamed mucosa

Heatmap of bacterial relative abundance in inflamed and noninflamed mucosa (control) at the genus level, which indicated the similarity and variation in the community composition of the samples (Fig. 4A). Genera clustered on the upper-left quadrant exhibited declining abundances from control to inflamed tissue. And genera clustered in the lower-left quadrant exhibited progressively increasing abundances from control to inflamed tissue. The LEfSe Clustergram showed Micrococcaceae, Bifidobacteriaceae, Bifidobacteriales, Flavobacteriaceae, and Methylobacteriaceae as potential microbiota significantly contributing to group differences between the noninflamed and inflamed mucosa of CD patients (Fig. 4B). A colored dot represents a significantly different microbe, and the name of the corresponding microbe is shown in the legend on the top right. The yellow nodes represented groups of microbes that did not play an important role in the different mucosa. Working from inside toward the outside, each circle represented the species at the phylum, class, order, family, and genus levels. The LEfSe-LDA showed Bifidobacteriales, Bifidobacterium, Bifidobacteriaceae, Micrococcaceae, Flavobacteriaceae, Gardnerella, Methylobacteriaceae, Rothia, Methyobacterium, Leucobacter, Actinomyces, Shinella, and Capnocytophaga as potential microbiota significantly contributing to group differences between the noninflamed and inflamed mucosa of CD patients. The default preset value was 2.0 (only the absolute value of LDA greater than two will be shown in the figure). The colors in the bar chart distinguished each group, while the length represented the LDA score. This score indicated the impact size of the species with significant differences between groups (Fig. 4C). The Wilcox test results that were used to assess differences in genera indicated that the relative abundance of Methylobacterium, Rothia, Shinella, Capnocytophaga, Actinomyces, Gardnerella, Leucobacter, and Bifidobacterium was significantly increased in the inflamed mucosa of CD patients compared with the noninflamed mucosa. The histogram of the relative abundance of each group was shown on the left. Listed in the center is the log2 value of the average relative abundance ratio of the same genus between the two mucosa. The figure on the right showed the P-value and FDR values obtained by the Wilcox test. If the P-value and FDR values were less than 0.05, the genus differed significantly between the two groups (Fig. 4D).

Functional prediction of bacteria with upregulated abundance in inflamed mucosa

Function prediction via KEGG indicated that bacteria in the mucosa of CD patients participated in functions such as amino acid transport and metabolism, translation, ribosomal structure, and biogenesis, general function prediction only, carbohydrate transport and metabolism, cell wall/membrane/envelope biogenesis, transcription, energy production and conversion (Fig. 5A). KEGG functional prediction shows the predicted function of the colony based on the KEGG database. The horizontal axis represents the sample and the vertical axis represents the relative abundance of the predicted function.

Functional prediction of bacteria with upregulated abundance in inflamed mucosa indicated that there were no significantly different pathways between these two groups. (Fig. 5B). The left side shows the relative abundance histogram of each group; the middle is the log2 value of the relative abundance mean ratio of the same channel between the two groups; the p-value and FDR value obtained by the Wilcox test are listed on the right. If the P-value and FDR value were less than 0.05, then the pathway was considered significantly different between the two groups.

Spearman correlation coefficient analysis of species

The genera were differentiated by the rank sum test, and a heatmap showing the Spearman correlation among dominant genera was generated by R software (Fig. 6A). Through this heatmap, important patterns and relationships among dominant genera can be found. The graph shows correlations between genera. A darker color indicates more closely correlated. The strongly positively correlated (correlation coefficient ≥0.6) bacterial genera were Dorea-Ruminococcus (0.65) (correlation coefficient), Haemophilus- Fusobacterium (0.62), Megamonas-Faecalibacterium (0.64), Blautia-Dorea (0.76), Blautia-Sutterella (0.64), Akkermansia-Ruminococcus (0.60), Coprococcus-Dorea (0.6), Coprococcus-Blautia (0.8), Neisseria-Fusobacterium (0.64), and Neisseria-Haemophilus (0.64). There were no strongly negatively correlated (absolute value of correlation coefficient ≥ 0.6) genera. The top 5 of the genera with moderate negative correlation (absolute correlation coefficient ≥0.3 and <0.6) were Akkermansia-Veillonella (−0.55), Coprococcus-Escherichia (−0.47), Rahnella,-Haemophilus (−0.46), Clostridium-Megamonas (−0.46), and Sphingomonas- Pseudomonas (−0.45). Genera that were negatively correlated include. Figure 6B showed the network map of genera. Each node in the graph represented an individual genus. The size of the node represented the magnitude of the average relative abundance of the genera. Straight lines were used to link genera, with pink indicating a positive correlation, blue indicating a negative correlation, and the thickness of the lines indicating the magnitude of the correlation. This figure only showed the relationship between genera with a correlation coefficient greater than 0.2.