Front. Microbiol.Frontiers in MicrobiologyFront. Microbiol.1664-302XFrontiers Media S.A.10.3389/fmicb.2022.875238MicrobiologyOriginal ResearchMucosa-Colonizing Microbiota Correlate With Host Autophagy Signaling in Patients With Inflammatory Bowel DiseaseWangWenxue12*†LiuZhongjian3†YueWei1†ZhuLing1ZhongHuijie12YangChao4HeTian4WanPing4*GengJiawei125*1Department of Infectious Disease and Hepatic Disease, First People’s Hospital of Yunnan Province, Affiliated Hospital of Kunming University of Science and Technology, Kunming, China2School of Medicine, Kunming University of Science and Technology, Kunming, China3Institute of Basic and Clinical Medicine, First People’s Hospital of Yunnan Province, Affiliated Hospital of Kunming University of Science and Technology, Kunming, China4Department of Gastroenterology, First People’s Hospital of Yunnan Province, Affiliated Hospital of Kunming University of Science and Technology, Kunming, China5Faculty of Life Science and Technology, Kunming University of Science and Technology, Kunming, China
Edited by: Junjun Wang, China Agricultural University, China
Reviewed by: Zhaolai Dai, China Agricultural University, China; Junling Shi, Northwestern Polytechnical University, China
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Both bacteria and autophagy are implicated in inflammatory bowel disease (IBD) pathogenesis. However, how bacteria crosstalk with autophagy signaling remains largely known, especially in intestinal mucosa. This study aimed to profile the internal complex autophagy signaling cascade and their external correlation with these bacteria, and consequently provide a systematic and precise target for future IBD diagnosis and therapy. We found the Ulcerative colitis (UC) patients exhibited more severe dysbiosis than the Crohn’s disease (CD) patients, as represented by alpha diversity, community phenotypes, and functional annotation compared with the control population. Meanwhile, CD patients showed greater transcriptional signaling activities of autophagy, endoplasmic reticulum (ER) stress, and bile acid production. Dominant bacteria (e.g., Rhodococcus, Escherichia, Shigella, and Enterococcus) were positively correlated and low-abundance bacteria (e.g., Bacillus, Acidovorax, Acinetobacter, and Stenotrophomonas) were negatively correlated with the autophagy signaling cascade (184 autophagy genes, 52 ER stress genes, and 22 bile acid production genes). Our observations suggested UC patients showed temporary and widespread microbiota turbulence and CD patients showed processive and local autophagy activity during IBD progression. Intestinal mucosa-colonizing bacteria were correlated with the bile/ER stress/autophagy signaling axis in IBD pathogenesis.
autophagyER stressbilemicrobiometranscriptomeinflammatory bowel disease8186043782160514202005AE160010YNWR-MY-2019-0722019LCZXKF-XH052020LCZXKF-XH012021LCZXXF-XH012021LCZXXF-XH15202101AT070275202101AY070001-2362018FE001-130National Natural Science Foundation of China10.13039/501100001809Yunnan Province Innovation Team of Intestinal Microecology-Related Disease Research and Technological TransformationEminent Doctors Program of Yunnan ProvinceYunnan Digestive Endoscopy Clinical Medical Center FoundationFundamental Research Projects of Yunnan ProvinceIntroduction
Intestinal microbes and their metabolic products play key roles in human diseases, especially inflammatory bowel disease (IBD). Meanwhile, both adherent-invasive and diffusely adherent bacteria stick to the intestinal mucosa and induce host immune cell activity and are linked to Crohn’s disease (CD) and ulcerative colitis (UC), respectively (Mirsepasi-Lauridsen et al., 2019). Notably, mucosa samples can directly reflect bacteria-host interactions during IBD progression. However, few multi-omics studies have investigated these interactions using mucosal biopsies. Furthermore, both the classification and mechanistic dissection of pathways involved in IBD remain challenging due to the complex correlative and interactive networks between host genetics and microbes (Plichta et al., 2019). Therefore, mucosal biopsy specimens are required to further explore these networks between host genetic factors and microbiota.
Genome-wide association studies have indicated that autophagy is an important mechanistic dissection of IBD pathogenesis, especially that of CD (Nguyen et al., 2013). Autophagy not only decreases intestinal epithelial permeability by inducing lysosomal degradation of tight junction proteins (Nighot et al., 2015), but also modulates programmed cell death in the intestinal epithelium (Matsuzawa-Ishimoto et al., 2017). In addition, autophagy is responsible for the elimination of intercellular bacteria from endoplasmic reticulum (ER) function defect-induced intestinal barrier leakage. For instance, autophagy can be induced by Autophagy related 16 like 1 (ATG16L1)-Transmembrane protein 59 (TMEM59) interactions in response to bacterial infection (Boada-Romero et al., 2016). Furthermore, activation of the bile acid receptor can strongly suppress the induction of autophagy (Lee et al., 2014). Thus, the profiles of mucosal bacteria-correlated autophagy gene networks and their interactions with ER stress and bile production signaling should be clarified in IBD patients.
Here, we collected intestinal mucosal biopsies from control population and IBD patients for dual-omics analysis of mucosal transcriptome and mucosa-colonizing bacterial diversity. Analysis focused on the autophagy signaling cascade and their complex correlations with intestinal mucosa-colonizing bacteria in IBD patients and control population. Combined with previous studies, our analysis suggested that mucosal bacteria regulated both ER stress and bile acid production and consequent autophagic activity, finally acting on host IBD progression. We systematically profiled the mucosal bacteria-autophagy correlation network in IBD patients, especially in those with CD, and provided a reliable analysis model of microbe-correlated intestinal diseases.
Materials and MethodsStudy Design
The control populations and IBD patients were recruited from the First People’s Hospital of Yunnan Province, China, from October 2017 to December 2019. Each participant provided information on age, occupation, and smoking and alcohol drinking history. The control population had no history of digestive tract disease or serious medical illnesses. Each IBD patient was diagnosed with UC or CD and had received no IBD therapy or antibiotics (within 3 months). All participants provided signed informed consents and completed a questionnaire regarding their age, sex, occupation, and antibiotic use, with family assistance if necessary.
Diagnosis of IBD
We used clinical feature inquiry, laboratory examination, and endoscopy and biopsy histopathological analysis to diagnose UC or CD according to clinical practice guidelines (Gomollón et al., 2017; Magro et al., 2017; Inflammatory Bowel Disease Group, Chinese Society of Gastroenterology, Chinese Medical Association, 2021). Firstly, outpatients with enteric symptoms of abdominal pain, diarrhea, blood in stools, and loose stools for more than 6 weeks, and/or parenteral manifestations, such as fatty liver and cholecystitis, were suggested for laboratory examination, which included a full blood count, electrolyte, liver function, and inflammatory marker [C-reactive protein (CRP)] tests. Secondly, when clinical features and laboratory tests suggested further investigation, an endoscopy (colon and/or intestinal endoscopy) was performed to obtain evidence and collect mucosal biopsies for histopathological analysis. Thirdly, histopathological investigations were performed to provide reliable information on mucosal architecture, lamina propria cellularity, neutrophil granulocyte infiltration, and epithelial abnormality. Finally, UC or CD was diagnosed when comprehensive analysis of clinical features, laboratory examinations, endoscopy, and histopathology gave positive responses.
Exclusion Criteria and Control Population Selection
To eliminate the possible effects of other intestinal diseases on UC and CD diagnosis, we carefully discriminated infectious enteritis, Clostridium difficile infection, intestinal tuberculosis, Behçet syndrome, amebiasis, schistosomiasis, and other intestinal diseases from IBD. Those participants without IBD or other digestive tract diseases and without serious diseases of other tissues and organs were chosen as the control population group. A total of 138 participants were initially recruited, 38 of which were excluded according to exclusion criteria (i.e., UC and/or CD and antibiotic use within the last 3 months). Finally, 25 control population, 26 CD patients, and 51 UC patients were enrolled in the study.
Sample Collection
The IBD patients and control population were first given an intravenous injection of propofol and etomidate, and then an enteroscopy was performed for mucosal sampling. All samples were collected from October 2017 to December 2019. Upon collection, the mucosal samples were immediately placed on ice and frozen at −80°C within 1 h for microbiome and transcriptome analyses.
Bacterial DNA Isolation and 16S rRNA Gene Sequencing
Microbial community genomic DNA was extracted from intestinal mucosa using a QIAamp DNA Mini Kit (Qiagen, Germany) according to the manufacturer’s instructions. The DNA extract was checked on 1% agarose gel, and DNA concentration and purity were determined using a NanoDrop 2000 UV–vis spectrophotometer (Thermo Scientific, Wilmington, United States). The hypervariable V3–V4 region of the bacterial 16S rRNA gene was amplified with primer pairs 338F (5′-ACTCCTACGGGAGGCAGCAG-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) using an ABI GeneAmp® 9700 PCR thermocycler (ABI, CA, United States). PCR amplification of the 16S rRNA gene was performed as follows: initial denaturation at 95°C for 3 min, followed by 27 cycles of denaturing at 95°C for 30 s, annealing at 55°C for 30 s, extension at 72°C for 45 s, single extension at 72°C for 10 min, and end at 4°C. The PCR mixture contained 5× TransStart FastPfu buffer 4 μl, 2.5 mM dNTPs 2 μl, forward primer (5 μM) 0.8 μl, reverse primer (5 μM) 0.8 μl, TransStart FastPfu DNA Polymerase 0.4 μl, template DNA 10 ng, and finally ddH2O up to 20 μl. PCR was performed in triplicate. The PCR products were extracted from 2% agarose gel and purified using an AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, United States) according to the manufacturer’s instructions and quantified using a Quantus™ Fluorometer (Promega, United States).
Purified amplicons were pooled in equimolar concentrations and paired-end sequenced on the Illumina MiSeq PE300/NovaSeq PE250 platforms (Illumina, San Diego, United States) using standard protocols.
The raw 16S rRNA gene sequencing reads were demultiplexed, quality-filtered using fastp v0.20.0 (Chen et al., 2018), and merged using FLASH v1.2.7 (Magoč and Salzberg, 2011) with the following criteria: (i) 300-bp reads were truncated at any site receiving an average quality score of <20 over a 50-bp sliding window, with truncated reads shorter than 50 bp and reads containing ambiguous characters discarded; (ii) only overlapping sequences longer than 10 bp were assembled according to their overlapping sequence. The maximum mismatch ratio of the overlapping region was 0.2. Reads that could not be assembled were discarded; and (iii) Samples were distinguished according to the barcode and primers, and the sequence direction was adjusted, exact barcode matching, two nucleotide mismatches in primer matching.
Operational taxonomic units (OTUs) with 97% similarity cutoff (Stackebrandt and Goebel, 1994; Edgar, 2013) were clustered using UPARSE v7.1, and chimeric sequences were identified and removed. The taxonomy of each OTU representative sequence was analyzed using RDP Classifier v2.2 (Wang et al., 2007) against the 16S rRNA database (Silva v132) with a confidence threshold of 0.7. To remove possible contamination, we sequenced three samples of pure water as a negative control with the same procedures, including DNA extraction, PCR amplification, cDNA library construction, and final sequencing. Based on the negative control sequencing results, the same DNA sequence count detected in the negative control was extracted from all control and IBD samples. To confirm the reliability of 16S rRNA sequencing data, we also performed real-time PCR and found the data are replicable.
Bioinformatics Analysis of 16S rRNA Gene Sequencing DataAlpha Diversity Analysis
Community richness (Sobs) and diversity (Shannon) indices were used to estimate α diversity and examined using the Welch’s t-test. Relative abundance of the top 47 genera was displayed using a histogram.
Community Hierarchical Clustering
Community clustering heatmap was conducted to show community variation at the genus level (top 47). Multiple sample comparisons were performed using Tukey–Kramer one-way ANOVA with false-discovery rate (FDR) correction and 95% confidence. Significance was assumed for adjusted values of p ≤ 0.05 (R software, vegan package).
Analysis of Abundance Differences
Significant differences between control population and IBD patients were determined using the Kruskal-Wallis H test at genus level. The FDR and Tukey–Kramer methods (CI = 0.95) were used for multiple testing correction and post hoc tests, respectively.
Functional Microbial Composition Analyses
The Functional Annotation of Prokaryotic Taxa (FAPROTAX) database extrapolates functions of cultured prokaryotes to estimate metabolic and other ecologically relevant functions (Louca et al., 2016). We used FAPROTAX to assess how the metabolic activities of intestinal mucosa-colonizing bacteria affect the host signaling network using the Tukey–Kramer method (CI = 0.95).
Phenotypic Prediction of Intestinal Mucosa-Colonizing Bacteria
Intestinal mucosa-colonizing bacteria phenotypes were predicted and compared using BugBase (Ward et al., 2017) and the Kruskal-Wallis H test. Briefly, BugBase uses an OTU table as an input file, which is normalized by the predicted 16S copy number. The preprocessed database and BugBase tool then automatically select thresholds to predict bacterial phenotypes.
RNA Extraction and Sequencing Data ProcessingRNA Extraction
Total RNA was extracted from tissue using TRIzol® Reagent according to the manufacturer’s instructions (Invitrogen), and genomic DNA was removed using DNase I (Takara). RNA quality was then determined using a 2100 Bioanalyzer (Agilent) and quantified using a ND-2000 spectrophotometer (NanoDrop Technologies). Only high-quality RNA samples (OD260/280 = 1.8–2.2, OD260/230 ≥ 2.0, RIN ≥ 6.5, and 28S:18S ≥ 1.0, >1 μg) were used to construct the sequencing library.
Library Preparation and Illumina Hiseq X Ten/NovaSeq 6000 Sequencing
The RNA-seq transcriptome library was prepared using a TruSeq™ RNA Sample Preparation Kit (Illumina, San Diego, CA, United States) with 1 μg of total RNA. In brief, messenger RNA (mRNA) was isolated according to the polyA selection method by oligo(dT) beads and then fragmented using fragmentation buffer. The double-stranded cDNA was synthesized using a SuperScript Double-Stranded cDNA Synthesis Kit (Invitrogen, CA, United States) with random hexamer primers (Illumina). The synthesized cDNA was then subjected to end repair, phosphorylation, and “A” base addition according to the library construction protocols of Illumina. Libraries were size-selected for 300-bp cDNA target fragments on 2% low-range ultra-agarose gel, followed by PCR amplification using Phusion DNA polymerase (NEB) for 15 PCR cycles. After quantification by TBS380, the paired-end RNA-seq library was sequenced using an Illumina HiSeq X Ten/NovaSeq 6000 sequencer (2 × 150-bp read lengths).
Read Mapping
The raw paired-end reads were trimmed and quality controlled using SeqPrep1 and Sickle2 with default parameters. The clean reads were then separately aligned to the reference genome in orientation mode using HISAT23 software (Kim et al., 2015). The mapped reads of each sample were assembled using StringTie4 with a reference-based approach (Pertea et al., 2015).
Differential Expression Analysis and Functional Enrichment
To identify differentially expressed genes (DEGs) between two samples, the expression level of each transcript was calculated according to the transcripts per million reads (TPM) method. RSEM5 was used to quantify gene abundances (Li and Dewey, 2011). Essentially, differential expression analysis was performed using DESeq2 (Love et al., 2014)/EdgeR (Robinson et al., 2010) with Q value ≤ 0.05 and DEGs with |log2FC| > 1 and Q value ≤ 0.05 (DESeq2 or EdgeR)/Q value ≤ 0.001 (DEGseq) deemed significant. In addition, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses were performed to identify DEGs significantly enriched in GO terms and metabolic pathways (Bonferroni-corrected value of p ≤ 0.05) compared to the whole transcriptome background. GO functional enrichment and KEGG pathway analysis were carried out using GOATOOLS6 and KOBAS (http://kobas.cbi.pku.edu.cn/home.do; Xie et al., 2011).
Alternative Splicing Event Identification
All alternative splicing events that occurred in the samples were identified using the recently released program rMATS (http://rnaseq-mats.sourceforge.net/index.html; Shen et al., 2014). Only isoforms similar with the reference or comprising novel splice junctions were considered; and splicing differences were detected as exon inclusion and exclusion, alternative 5′ and 3′, and intron retention events.
Transcriptome AnalysisTime Series Expression Trend Analysis
STEM is used for clustering, comparing, and visualizing time series gene expression data (Ernst and Bar-Joseph, 2006). The STEM clustering algorithm is a supervised algorithm, that is, clustering is classified into artificially set trends. First, the software simulates n of the most representative possible trends according to the preset, and then calculates the correlation coefficient between each gene and the preset trends. Finally, each gene is classified into the trend to which it is most similar. We used STEM (v1.3.11) to explore the gene expression pattern of the healthy controls and CD and UC patients. Analysis explored functional enrichment of genes with a certain expression pattern and predicted the genetic regulatory network of intestinal mucosa.
Gene Expression Correlation Analysis
Genes obtained from expression trend analysis and/or related to autophagy, ER stress, and bile acid production were used to calculate Spearman correlation coefficients based on the correlation of gene expression. In the correlation networks, the larger the node, the greater the number of correlations between the expression of the gene and other genes.
Protein–Protein Interaction Analysis
We used the STRING database7 to perform Protein–Protein Interaction (PPI) network analysis of genes of interest. The interactions corresponding to genes of interest were directly extracted from the database to construct the network. NetworkX (Python) was used to visualize the network of the genes of interest.
KEGG Functional Enrichment Analysis
We used R script and Fisher’s precision probability test to perform KEGG pathway enrichment analysis. To control the false positives of enrichment analysis, multiple testing corrections (Benjamini-Hochberg, BH) were carried out. KEGG pathways reaching a corrected p value of 0.05 were defined as significantly enriched.
Ethics Approval
All study protocols and procedures were approved by the Medical Ethics Board of the First People’s Hospital of Yunnan Province (GXBSC-2021001, 2021 updated), China, and were carried out in accordance with all relevant provincial, national, and international guidelines, including the Declaration of Helsinki. Written informed consent was obtained from all participants prior to their inclusion in the study.
Statistical Analysis
GraphPad Prism v9.0.0 (San Diego, CA, United States) and the R stats package were used to analyze all data (R Core Team, 2013). The Kruskal-Wallis H test with FDR correction (CI = 0.95) for multiple comparisons was conducted to compare abundant bacterial taxa. Spearman rank-order correlation was used to evaluate associations between genus-level bacterial relative abundances and autophagy-related gene expression levels. In all analyses, p < 0.05 was considered statistically significant.
ResultsCharacteristics of IBD Patients and Controls
We initially recruited 138 participants, 38 of which were excluded due to antibiotic use within the last 3 months and/or suffering from other serious illnesses. In total, 23 healthy control population, 26 CD patients, and 51 UC patients were enrolled (Table 1). Most IBD patients were male (68 out of 77 subjects) and had no history of smoking. The control population was recruited from outpatients suffering abdominal discomfort, but with no gastrointestinal or other serious illness. For transcriptome analysis, an additional biopsy sample was collected from six control population, five CD patients, and 12 UC patients (Supplementary Table S1) with informed consent.
Basic characteristics of all IBD patients.
Index
Control(n = 23)
CD(n = 26)
UC(n = 51)
Gender (male: female)
11:12
15:11
28:23
Age (mean ± SD)
47.48 ± 11.14
38.77 ± 9.40
44.86 ± 11.13
Smoking (yes: no)
0:23
9:17
11:40
Alcoholic drinking (yes: no)
11:12
7:19
34:17
Antibiotic use (within 1 month)
No
No
No
Occupation
Farmer
7
6
8
Factory worker
8
7
11
Office staff
3
3
21
Teacher
0
4
3
Others
5
6
8
Chief complaint
Abdominal pain
5
15
13
Diarrhea
0
2
16
Blood in stools
18
3
19
Others
0
6
3
IBD Patients, Especially UC Patients, Show Intestinal Mucosa-Colonizing Microbial Dysbiosis and Corresponding Dysfunction
The intestinal microbiota play a key role in IBD progression. IBD patients exhibit clear dysbiosis, mainly represented by a decrease in bacterial diversity (Ott et al., 2004; Manichanh et al., 2006). Mucosal biopsies are advantageous as they can be used to explore the real interactions between dysbiosis and host responses. Here, we collected intestinal mucosal biopsies from control population and IBD patients to investigate the compositional movement and functional variation potential of intestinal mucosa-colonizing bacteria. Consistent with previous study (Ott et al., 2004), we found that alpha diversity (Sobs and Shannon indices) of the intestinal mucosal bacterial community was significantly decreased in IBD patients compared with the control population. Furthermore, UC patients showed poorer alpha diversity than CD patients (Figures 1A,B). The decrease in the Sobs index mainly originated from the reduced abundances of Proteobacteria and Bacteroidetes (Supplementary Figure S1). Correspondingly, FAPROTAX analysis indicated that the metabolic activities of nitrite, nitrate, and fumarate were significantly decreased in IBD patients, especially UC patients (Figure 1D; Supplementary Figure S3). Notably, these metabolic activities had brought focuses because of their key roles in IBD incidence, prevention, and therapy (Dykhuizen et al., 1996; Saijo et al., 2010; Jädert et al., 2014; Casili et al., 2016). Redox of these nitrogenous salts is also associated with IBD pathology (Bourgonje et al., 2020). Based on community phenotype analysis, we also found that IBD patients showed a decrease in oxidative stress tolerance, primarily due to the lower abundances of Escherichia, Shigella, and Pseudomonas (Figure 1E; Supplementary Figures S4A,B). Unexpectedly, pathogen causing diarrhea and gastroenteritis potentials also decreased significantly in IBD patients. Potential pathogens, such as Rhodococcus, Streptococcus, Enterococcus, Veillonella, Ruminococcus torques, and Ruminococcus gnavus, were enriched in the intestinal mucosa of IBD patients (Figure 1C; Supplementary Figure S2), which contributed to high Gram-positive performance, another community phenotype (Supplementary Figures S4C,D).
Intestinal mucosa-colonizing bacteria and their potential functional profiles and community phenotypes in control population and inflammatory bowel disease (IBD) patients. Alpha diversity of intestinal mucosa-colonizing microbiota based on Sobs index (A) and Shannon index (B). (C) Kruskal-Wallis H test compared average composition of intestinal mucosa-colonizing bacteria (top 8) at genus level. Multiple testing correction: false discovery rate (FDR); Post hoc test: Tukey–Kramer (CI = 0.95). All bacteria were named to genus level unless otherwise noted in brackets. (D) Functional Annotation of Prokaryotic Taxa (FAPROTAX) approach predicted potential functional profiles of intestinal mucosa-colonizing bacteria using Kruskal-Wallis H test. Multiple testing correction: FDR; Post hoc test: Tukey–Kramer (CI = 0.95). (E) BugBase predicted community phenotypes of intestinal mucosa-colonizing bacteria using Kruskal-Wallis H test. Multiple testing correction: FDR; Post hoc test: Tukey–Kramer (CI = 0.95). * 0.01 < p ≤ 0.05, **0.001 < p ≤ 0.01, and *** p ≤ 0.001. Control population: n = 23; CD patient: n = 26; and UC patient: n = 51.
Active Intestinal Mucosa Signaling in CD Patients Exceeds That in UC Patients
Functional genomics network analysis is a powerful tool for identifying the key regulatory networks involved in IBD progression (Peters et al., 2017). Here, gene expression clustering analysis revealed that the intestinal mucosa signaling network was comprised of eight clusters based on differential expression in control population and UC and CD patients (Figure 2A). Three clusters, i.e., cluster 4 (3,339 genes, p = 0.0000012), cluster 6 (6,355 genes, p = 1.2e–156), and cluster 7 (9,773 genes, p = 0), were significant according to the differential transcription levels in the control population and UC and CD patients. Cluster 6 genes were highly expressed in the intestinal mucosa of both UC and CD patients, whereas cluster 7 and 4 genes showed a gradual decrease in expression in UC patients (Figures 2B–D). Thus, these findings indicate that intestinal mucosa signaling activity was higher in CD patients than in UC patients.
Clustering profile of intestinal mucosa signaling in IBD patients. (A) Trend analysis chart showing gene expression profile in intestinal mucosa in IBD patients using a time course analysis. Chart shows all active genes were divided into eight clusters, three of which [i.e., cluster 4 (p = 0.0000012), cluster 6 (p = 1.2e–156), and cluster 7 (p = 0)] were significant. Time course clustering algorithm: STEM, p < 0.05. Gene expression heatmaps of cluster 4 (B), cluster 6 (C), and cluster 7 (D), which contained 3,339, 6,355, and 9,773 genes, respectively. Con: control population, n = 6; UC: ulcerative colitis patients, n = 12; and CD: Crohn’s disease patients, n = 5.
Environmental Factor Correlation Analysis Reveals Bacterial Community-Matched Autophagy Signaling in Intestinal Mucosa of IBD Patients
Increasing studies emphasize that crosstalk between the mucosal microbiome and host signaling network greatly affects IBD progression and clinical outcome (Morgan et al., 2015; Ryan et al., 2020). Recent review work provides new insights of the interplay between autophagy and intestinal bacteria and suggested that IBD-associated autophagy alleles and their interactions with environmental triggers, such as resident microbiota, are crucial for developing new therapeutic strategies for IBD treatment (Larabi et al., 2020). Thus, we performed correlation analysis between intestinal mucosa-colonizing bacteria and autophagy gene expression levels in clusters 6 and 7, which covered most KEGG autophagy-related genes (Supplementary Figure S7). In this analysis, 144 autophagy genes (58 in cluster 6 and 86 in cluster 7) were significantly correlated with intestinal mucosa-colonizing bacteria (Figures 3A,D). Gene expression heatmap analysis showed that these bacteria-correlated autophagy genes were more highly expressed in CD patients than in UC patients (Supplementary Figure S8). Furthermore, these autophagy genes not only formed a correlation network, but also a PPI network. For example, bacteria-correlated autophagy genes Autophagy related 4A cysteine peptidase (ATG4A), Autophagy related 5 (ATG5), ATG16L1, Autophagy related 9A (ATG9A), Signal transducer and activator of transcription 3 (STAT3), NRAS proto-oncogene, GTPase (NRAS), and Heat shock protein 90 alpha family class A member 1 (HSP90AA1) in cluster 6 (Figures 3B,C) and Beclin 1 (BECN1), Mitogen-activated protein kinase 1 (MAPK1), Phosphatidylinositol 3-kinase catalytic subunit type 3 (PIK3C3), Phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit alpha (PIK3CA), Charged multivesicular body protein 2A (CHMP2A), Charged multivesicular body protein 3 (CHMP3), Charged multivesicular body protein 4A (CHMP4A), Vacuolar protein sorting 4 homolog A (VPS4A), Vacuolar protein sorting 25 homolog (VPS25), Translocase Of outer mitochondrial membrane 20 (TOMM20), and AKT serine/threonine kinase 2 (AKT2) in cluster 7 (Figures 3E,F) were core members of both the correlation and interaction networks.
Correlation pattern between intestinal mucosa-colonizing microbiota and host autophagy signaling in IBD patients. Spearman correlation heatmap shows that 58 autophagy-related genes in cluster 6 (A) and 86 autophagy-related genes in cluster 7 (D) were significantly correlated with intestinal mucosa-colonizing bacteria (top 47). All bacteria were named to genus level unless otherwise noted in brackets. Expression correlation network of bacteria-patterned autophagy genes in cluster 6 (B) and cluster 7 (E). Circle area is positively correlated with number of connected genes. Correlation network was constructed based on Spearman rank correlation coefficients (∣Spearman Coef∣ ≥ 0.8, p < 0.05). Multiple testing correction: Benjamini and Hochberg (BH). Protein–protein interaction network of bacteria-patterned autophagy genes in cluster 6 (C) and cluster 7 (F). Circle or equilateral triangle areas are positively correlated with number of connected genes. Interaction between circle-labeled genes with others has been reported. Interaction between equilateral triangle-labeled genes with others was predicted based on primary structure of gene-coding proteins. p ≤ 0.05. Control population: n = 6; UC patient: n = 12; and CD patient: n = 5.
The above bacteria could also be clustered into two groups, i.e., those showing positive correlations with autophagy genes in both cluster 6 and cluster 7 (e.g., Haemophilus, Ruminococcus torques, and Rhodococcus) and those showing negative correlations (e.g., Bacillus, Acidovorax, Acinetobacter, and Stenotrophomonas; Figures 3A,D). It is worth noting that most of the intestinal mucosa-colonizing bacteria showing positive correlation with the above autophagy genes have been reported in previous IBD studies (Goyal et al., 2018; Soltys et al., 2020). For those bacteria negatively correlated with the autophagy genes, Bacillus is recognized as a probiotic bacterium in IBD prevention and therapy (Ghavami et al., 2020; Shinde et al., 2020; Zhang et al., 2020; Liu et al., 2021), but the remaining bacteria (e.g., Acidovorax, Acinetobacter, and Stenotrophomonas) were first reported to involve in IBD progression. Interestingly, the co-occurrence network indicated that low-abundance intestinal mucosa-colonizing bacteria were negatively and densely associated with autophagy genes, whereas high-abundance bacteria were positively and sparsely associated with autophagy genes (Supplementary Figure S6).
The autophagy genes in cluster 4, which were highly expressed in CD patients but not in UC patients compared with control population, showed more positive correlations and fewer negative correlations with the intestinal mucosa-colonizing bacteria than found in clusters 6 and 7. The positive and negative correlations were represented by Haemophilus, Ruminococcus torques, and Rhodococcus, and Bacillus, Acidovorax, Acinetobacter, and Stenotrophomonas, respectively. These cluster 4 bacteria-correlated autophagy genes also showed correlation and interaction networks, both of which dominated by Tumor protein P53 (TP53), Stimulator of interferon response CGAMP interactor 1 (STING), and interferon gamma inducible protein 16 (IFI16). These results are highly consistent with previous study suggesting that STING senses bacterial viability to orchestrate autophagy (Moretti et al., 2017). The cluster 4 co-occurrence network also showed more positive correlations between autophagy genes and intestinal mucosa-colonizing bacteria (e.g., Haemophilus, Ruminococcus torques, Enterococcus, Veillonella, Escherichia, and Shigella), but strong negative correlations between Bacillus and Streptococcus and autophagy genes [e.g., Ring finger protein 41 (RNF41), BNIP3P5, and Interleukin 10 (IL10); Supplementary Figure S5]. However, this correlation only existed in CD patients. Although similar correlations between mucosa-colonizing bacteria and autophagy activity have been reported (Nakagawa et al., 2004; Sudhakar et al., 2019; Wu et al., 2019), we used clinical biopsies to explore real crosstalk between mucosal bacteria and autophagy genes during IBD progression and to clarify the complex signaling pathways connecting intestinal microbes and host disease pathogenesis.
Autophagy-Induced ER Stress Is Correlated With Intestinal Mucosa-Colonizing Bacteria in IBD Patients
Moretti et al. (2017) found that STING initiates autophagy after bacterial infection and that this autophagy is mediated by ER stress. ER stress also regulates autophagy processes during IBD progression (Kaser and Blumberg, 2009; Kaser and Blumberg, 2011; Grootjans et al., 2019). According to our analysis, 15 ER stress-related genes in cluster 6 were significantly correlated with intestinal mucosa-colonizing bacteria (Figure 4A; Supplementary Figure S12A), and the main correlations were presented with a co-occurrence network [e.g., negative correlation between Acinetobacter and Heat shock protein family A (Hsp70) member 1A (HSPA1A); Supplementary Figure S10A]. ER stress-related genes in cluster 6, such as HSPA1A, Heat shock protein family A (Hsp70) member 5 (HSPA5), Heat shock protein 90 beta family member 2 (HSP90B2I), Protein disulfide isomerase family A member 3 pseudogene 1 (PDIA3P1), and Phorbol-12-myristate-13-acetate-induced protein 1 (PMAIP1), formed correlation and interaction networks (Supplementary Figures S10B,C). In addition, 31 ER stress-related genes in cluster 7 were significantly correlated with mucosal bacteria (Figure 4D; Supplementary Figure S12B), and formed correlation and interaction networks (Supplementary Figures S10E,F). The main correlations were also presented with a co-occurrence network (e.g., negative correlation of Acidovorax with BCL2 like 11 (BCL2L11) and Wolframin ER transmembrane glycoprotein (WFS1) in cluster 7; Supplementary Figure S10D). The correlation and interaction networks of ER stress-related genes in cluster 7 mainly consisted of Activating transcription factor 6 (ATF6), Homocysteine inducible ER protein with ubiquitin like domain 1 (HERPUD1), BAG cochaperone 6 (BAG6), BCL2 binding component 3 (BBC3), Mesencephalic astrocyte derived neurotrophic factor (MANF), DnaJ heat shock protein family (Hsp40) member B9 (DNAJB9), BCL2L11, Membrane bound transcription factor peptidase, site 1 (MBTPS1), Eukaryotic translation initiation factor 2 subunit alpha (EIF2S1), Valosin containing protein (VCP), Mitogen-activated protein kinase kinase kinase 5 (MAP3K5), and Thioredoxin domain containing 12 (TXNDC12; Supplementary Figures S10E,F). HSPA5-induced ATF6 activity is a typical representative of ER stress signaling in the intestinal epithelium (Stengel et al., 2020) that bridging clusters 6 and 7 profile in this study. Furthermore, combinative analysis of both correlation and interaction network, also exhibited the complexity of ER stress-autophagy signaling cascade in IBD (Figures 4B,C,E,F). In addition, we also found similar ER stress signaling based on KEGG functional pathway analysis (Supplementary Figure S11).
Correlation pattern of endoplasmic reticulum (ER) stress with intestinal mucosa-colonizing bacterial community and autophagy signaling in IBD patients. Spearman correlation heatmap shows 15 ER stress-related genes in cluster 6 (A) and 31 ER stress-related genes in cluster 7 (D) were significantly correlated with intestinal mucosa-colonizing bacteria (top 47). All bacteria were named to genus level unless otherwise noted in brackets. * 0.01 < p ≤ 0.05, ** 0.001 < p ≤ 0.01, and *** p ≤ 0.001. Correlation network of both autophagy signaling and ER stress-related genes in cluster 6 (B) and cluster 7 (E). Autophagy signaling genes and ER stress genes are labeled in blue and red, respectively. Circle area is positively correlated with number of connected genes. Correlation network was constructed based on Spearman rank correlation coefficients (∣Spearman Coef∣ ≥ 0.8, p < 0.05). Multiple testing correction: BH. Interaction network of autophagy signaling and ER stress genes in cluster 6 (C) and cluster 7 (F). Circle or equilateral triangle areas are positively correlated with number of connected genes. Interaction between circle-labeled genes with others has been reported. Interaction between equilateral triangle-labeled genes with others was predicted based on primary structure of gene-coding proteins. Control population: n = 6; UC patient: n = 12; and CD patient: n = 5.
In cluster 4, six ER stress-related genes [i.e., endoplasmic reticulum protein 27 (ERP27), NHL repeat containing E3 ubiquitin protein ligase 1 (NHLRC1), Fc gamma receptor IIb (FCGR2B), Stanniocalcin 2 (STC2), Caspase 17, pseudogene (CASP17P), and ChaC glutathione specific gamma-glutamylcyclotransferase 1 (CHAC1)] were significantly correlated with intestinal mucosa-colonizing bacteria and were highly expressed in CD patients but not in UC patients compared with the control population. The co-occurrence network showed that both STC2 and FCGR2B were negatively and positively correlated with Proteobacteria (Acinetobacter and Klebsiella), respectively. Interestingly, both STC2 and FCGR2B are involved in ER stress activities during IBD progression (Franke et al., 2016; Coope et al., 2019). In addition, co-occurrence network analyses also linked autophagy signaling of STING and IL10 with the ER stress activities of CASP17P and FCGR2B (Supplementary Figure S9).
ER Stress-Related Bile Acid Production Signaling Is Correlated With Mucosa-Colonizing Bacteria in IBD Patients
Increasing evidence indicates that intestinal dysbiosis induces bile acid dysmetabolism and consequent IBD progression (Duboc et al., 2013; Lloyd-Price et al., 2019; Quinn et al., 2020; Sinha et al., 2020). Bacterial-driven bile acid metabolites have also been investigated to describe the pathophysiological basis of bacteria-bile acid associations (Long et al., 2017; Lavelle and Sokol, 2020). Here, we expanded the association between bile acid production-related genes and intestinal mucosa-colonizing bacteria. Eleven genes [e.g., Leptin (LEP), cytochrome P450 family 46 subfamily A member 1 (CYP46A1), KIT proto-oncogene, receptor tyrosine kinase (KIT), Phospholipase A2 group IB (PLA2G1B), and Retinoid X receptor alpha (RXRA)] in cluster 6 and eight genes [e.g., Sulfotransferase family 2A member 1 (SULT2A1), Fibroblast growth factor receptor 4 (FGFR4), and fatty acid binding protein 1 (FABP1)] in cluster 7 were significantly correlated with the intestinal mucosa-colonizing bacteria (Figures 5A,D). Similar to autophagy and ER stress signaling, genes in both cluster 6 and cluster 7 formed correlation and interaction networks, dominated by RXRA (Supplementary Figures S14B,C) and FABP1 (Supplementary Figures S14D,E), respectively. Notably, KIT and PLA2G1B in cluster 6 (Supplementary Figure S14A) and SULT2A1 in cluster 7 (Supplementary Figure S14D) showed strong correlations with the intestinal mucosa-colonizing bacteria. Bile-regulated lipid metabolism can induce ER stress in intestinal epithelial cells and intestinal bacteria modulate this metabolic process (Ko et al., 2020). Here, we identified bile acid-related genes in cluster 6 [RXRA, Aldo-keto reductase family 1 member C2 (AKR1C2), and Oxysterol binding protein like 7 (OSBPL7)] and cluster 7 [FGFR4, Nuclear receptor coactivator 2 (NCOA2), and FABP1], which were all significantly correlated with intestinal mucosa-colonizing bacteria and formed both correlation and interaction networks with ER stress signaling (Figures 5B,C,E,F). KEGG functional analysis also revealed the upregulation of genes related to bile biosynthesis (Supplementary Figure S15) and secretion (Supplementary Figure S16). Furthermore, the expression differentiation of bile acid-producing gene between UC and CD patients increased in cluster 7 than in cluster 6 (Supplementary Figure S17).
Correlation pattern of bile acid production with intestinal mucosa-colonizing bacteria community and ER stress in IBD patients. Spearman correlation heatmap shows 11 bile-related genes in cluster 6 (A) and eight bile-related genes in cluster 7 (D) were significantly correlated with intestinal mucosa-colonizing bacteria (top 47). All bacteria were named to genus level unless otherwise noted in brackets. * 0.01 < p ≤ 0.05, ** 0.001 < p ≤ 0.01, and *** p ≤ 0.001. Correlation network of both ER stress- and bile acid production-related genes in cluster 6 (B) and cluster 7 (E). ER stress- and bile acid production-related genes are labeled in red and green, respectively. Circle area is positively correlated with number of connected genes. Correlation network was constructed based on Spearman rank correlation coefficients (∣Spearman Coef∣ ≥ 0.8, p < 0.05). Multiple testing correction: BH. Interaction network of ER stress- and bile acid production-related genes in cluster 6 (C) and cluster 7 (F). Circle or equilateral triangle areas are positively correlated with number of connected genes. Interaction between circle-labeled genes with others has been reported. Interaction between equilateral triangle-labeled genes with others was predicted based on primary structure of gene-coding proteins. Control population: n = 6; UC patient: n = 12; and CD patient: n = 5.
In the other hand, the bacteria that correlated with intestinal mucosa bile acid signaling, mainly comprised of Proteobacteria. The correlations between Stenotrophomonas, Acinetobacter, Bacillus, and Corynebacterium and cluster 6 genes (Figure 5A; Supplementary Figure S14A) and Rhodococcus and Acidovorax and cluster 7 genes (Figure 5B; Supplementary Figure S14B), exhibited the most significant correlations with host bile acid signaling. In addition, FGF19 in cluster 4, which is reported to modulate the connection between intestinal microbiota and host inflammation (Gadaleta et al., 2020), was correlated with multiple intestinal bacteria, including Staphylococcus, Pelomonas, Haemophilus, Ruminococcus torques, and Bacteroides (Supplementary Figure S13).
Discussion
Increasing studies focus on the mucosal microbiome and local host immune activities in IBD patients (Howell et al., 2018; Lloyd-Price et al., 2019; Ryan et al., 2020) and the role of autophagy in IBD progression, especially in CD patients (Khor et al., 2011; Larabi et al., 2020). To date, however, no studies have reported on how mucosal bacteria interact with host autophagy activity based on biopsy samples even that obviously hold a great interest of precise diagnosis and therapy of IBD. Here, we used intestinal mucosal biopsies to elucidate the interaction network between mucosal bacteria and host autophagy signaling, as this network can reveal the bacterial infection-related pathogenesis of IBD and highlight potential and precise therapeutic targets for IBD.
Our results showed that active mucosal signaling, including autophagy, could be divided into cluster 4, 6, and 7 based on signaling extensity in UC patients. In general, all genes in the three clusters were highly activated in CD patients, but only genes in cluster 6 were highly expressed in UC patients. In addition, genes in cluster 7 showed increased expression in UC patients compared with that in the control population. The genes in cluster 4 also showed the same activities in UC patients and control population. These results are consistent with previous study showing that autophagy genes are pivotal for intestinal homeostasis and antimicrobial function in IBD, especially CD (Larabi et al., 2020).
In our study, UC patients exhibited greater disturbance in mucosal microbial diversity, community phenotype, and underlying functional spectrum, whereas CD patients exhibited greater autophagy and associated signaling cascades. These differences may have originated from the different pathological durations of UC and CD in patients. The long disease course in CD patients can strongly induce autophagy signaling to fight pathogenic microbial invasion. In another hand, pathogenic microbial invasion can induce CD progression, whereas diffusion behavior tends to increase UC progression (Mirsepasi-Lauridsen et al., 2019). Certainly, mucosal biopsies are ideal samples to determine bacterial diffusion behavior. Previous studies using fecal samples found that CD patients show more extensive microbial community disturbance than UC patients (Halfvarson et al., 2017; Pascal et al., 2017). Based on our findings, we suggest that CD patients continuously repair their gut defenses and establish new gut immune hemostasis during the long-term pathological progression of the disease, which is different from normal states and UC patients.
We also found that most dominant bacteria were positively correlated with autophagy genes and related signaling cascades, whereas less abundant bacteria showed negative correlations. Dominant bacteria may form stronger community structures that can withstand serious environmental shifts or may change into opportunistic pathogens for survival. In contrast, low-abundance bacteria communities may be more sensitive to environmental shift, and therefore tightly related to host metabolism and healthy conditions and finally process a high sensitivity and great value for intestinal disease diagnosis and therapy, including IBD. However, the low-abundance bacteria also brought some uncertainty. For example, Acidovorax (<0.2% of total abundance) was negatively correlated with most autophagy genes and related signaling cascades, even these genes or signaling play an adversarial role in autophagy. In addition, we observed several bacterial species that were not significantly correlated with any host genes (data not shown); it was difficult to explain how they survive under severe immune stress from host and fierce competition from other microbes if they do not clearly show cooperative (or confrontational) properties.
Our study has several limitations. Both the recruited patients (average age > 45) and control population (average age > 35) were older than participants in other studies, which is likely due to the consciousness deficiency of IBD risk in the local young population. In this study, the microbiota-correlated autophagy genes concurrently showed a strong correlation with ER stress. Integrating previous study (Kökten et al., 2018), the autophagy in our study could be either macroautophagy or chaperone-mediated autophagy, but the specific subtype was not clear and needs to be further explored. The relationship between mucosal bacteria and host autophagy activity needs more mucosa biopsies to conform their correlation profiles, but our matched analysis of microbiome and transcriptome from same subjects and location still drew an interactive map of bacteria-autophagy-IBD progression and could provide more valuable diagnostic and therapeutic targets based on autophagy mechanism.
Conclusion
In conclusion, we found that UC patients exhibited more severe dysbiosis and the functional phenotype of intestinal mucosa-colonizing bacteria, whereas CD patients exhibited more active autophagic signaling compared to the controls. Dominant bacteria and low-abundance bacteria showed positive and negative correlations with host autophagy genes, respectively. In addition, correlation and interaction networks were found between bile acid production and ER stress, with the latter showing an interesting interaction with autophagy activity. Thus, the present study elucidated how the intestinal mucosa-colonizing bacterial community interacts with host bile/ER stress/autophagy signaling cascades during IBD progression, which could aid in disease diagnosis and autophagy-targeted therapy.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://www.ncbi.nlm.nih.gov/, PRJNA797846.
Ethics Statement
The studies involving human participants were reviewed and approved by the Medical Ethics Board of the First People’s Hospital of Yunnan Province (GXBSC-2021001, 2021 updated). The patients/participants provided their written informed consent to participate in this study.
Author Contributions
WW and JG designed the project and reviewed and revised the final version of the manuscript. WW and ZL analyzed the 16S rRNA sequencing data. LZ, WY, HZ, CY, TH, PW, and JG collected mucosal biopsies and clinical documents. WW, ZL, and JG finished the draft. WW, LZ, PW, and JG supervised the study and rendered foundation supports. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by the National Natural Science Foundation of China (81860437 and 82160514), the Yunnan Province Innovation Team of Intestinal Microecology-Related Disease Research and Technological Transformation (202005AE160010), Eminent Doctors Program of Yunnan Province (YNWR-MY-2019-072), Yunnan Digestive Endoscopy Clinical Medical Center Foundation [2X2019-01-02]-(2019LCZXKF-XH05, 2020LCZXKF-XH01, 2021LCZXXF-XH01, and 2021LCZXXF-XH15), and Fundamental Research Projects of Yunnan Province (202101AT070275, 202101AY070001-236, and 2018FE001-130).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
We thank WY for help with identification of the UC and CD patients.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2022.875238/full#supplementary-material
ReferencesBoada-RomeroE.Serramito-GómezI.SacristánM. P.BooneD. L.XavierR. J.Pimentel-MuiñosF. X.. (2016). The T300A Crohn’s disease risk polymorphism impairs function of the WD40 domain of ATG16L1. Nat. Commun.7:11821. doi: 10.1038/ncomms11821, PMID: 27273576BourgonjeA. R.FeelischM.FaberK. N.PaschA.DijkstraG.van GoorH. (2020). Oxidative stress and redox-modulating therapeutics in inflammatory bowel disease. Trends Mol. Med.26, 1034–1046. doi: 10.1016/j.molmed.2020.06.006, PMID: 32620502CasiliG.CordaroM.ImpellizzeriD.BruschettaG.PaternitiI.CuzzocreaS.. (2016). Dimethyl fumarate reduces inflammatory responses in experimental colitis. J. Crohns Colitis10, 472–483. doi: 10.1093/ecco-jcc/jjv231, PMID: 26690241ChenS.ZhouY.ChenY.GuJ. (2018). fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics34, i884–i890. doi: 10.1093/bioinformatics/bty560, PMID: 30423086CoopeA.PascoalL. B.BotezelliJ. D.da SilvaF. A. R.AyrizonoM. D. L. S.RodriguesB. L.. (2019). ER stress activation in the intestinal mucosa but not in mesenteric adipose tissue is associated with inflammation in Crohn’s disease patients. PLoS One14:e0223105. doi: 10.1371/journal.pone.0223105, PMID: 31557250DubocH.RajcaS.RainteauD.BenarousD.MaubertM.-A.QuervainE.. (2013). Connecting dysbiosis, bile-acid dysmetabolism and gut inflammation in inflammatory bowel diseases. Gut62, 531–539. doi: 10.1136/gutjnl-2012-302578, PMID: 22993202DykhuizenR. S.MassonJ.McKnightG.MowatA. N.SmithC. C.SmithL. M.. (1996). Plasma nitrate concentration in infective gastroenteritis and inflammatory bowel disease. Gut39, 393–395. doi: 10.1136/gut.39.3.393, PMID: 8949643EdgarR. C. (2013). UPARSE: highly accurate OTU sequences from microbial amplicon reads. Nat. Methods10, 996–998. doi: 10.1038/nmeth.2604, PMID: 23955772ErnstJ.Bar-JosephZ. (2006). STEM: a tool for the analysis of short time series gene expression data. BMC Bioinformatics7:191. doi: 10.1186/1471-2105-7-191, PMID: 16597342FrankeL.el BannoudiH.JansenD. T. S. L.KokK.TrynkaG.DiogoD.. (2016). Association analysis of copy numbers of FC-gamma receptor genes for rheumatoid arthritis and other immune-mediated phenotypes. Eur. J. Hum. Genet.24, 263–270. doi: 10.1038/ejhg.2015.95, PMID: 25966632GadaletaR. M.Garcia-IrigoyenO.CarielloM.ScialpiN.PeresC.VetranoS.. (2020). Fibroblast growth factor 19 modulates intestinal microbiota and inflammation in presence of Farnesoid X receptor. EBioMedicine54:102719. doi: 10.1016/j.ebiom.2020.102719, PMID: 32259714GhavamiS. B.YadegarA.AghdaeiH. A.SorrentinoD.FarmaniM.MirA. S.. (2020). Immunomodulation and generation of tolerogenic dendritic cells by probiotic bacteria in patients with inflammatory bowel disease. Int. J. Mol. Sci.21:6266. doi: 10.3390/ijms21176266, PMID: 32872480GomollónF.DignassA.AnneseV.TilgH.Van AsscheG.LindsayJ. O.. (2017). 3rd European evidence-based consensus on the diagnosis and management of Crohn’s disease 2016: part 1: diagnosis and medical management. J. Crohns Colitis11, 3–25. doi: 10.1093/ecco-jcc/jjw168, PMID: 27660341GoyalA.YehA.BushB. R.FirekB. A.SieboldL. M.RogersM. B.. (2018). Safety, clinical response, and microbiome findings following fecal microbiota transplant in children with inflammatory bowel disease. Inflamm. Bowel Dis.24, 410–421. doi: 10.1093/ibd/izx035, PMID: 29361092GrootjansJ.KrupkaN.HosomiS.MatuteJ. D.HanleyT.SaveljevaS.. (2019). Epithelial endoplasmic reticulum stress orchestrates a protective IgA response. Science363, 993–998. doi: 10.1126/science.aat7186, PMID: 30819965HalfvarsonJ.BrislawnC. J.LamendellaR.Vázquez-BaezaY.WaltersW. A.BramerL. M.. (2017). Dynamics of the human gut microbiome in inflammatory bowel disease. Nat. Microbiol.2:17004. doi: 10.1038/nmicrobiol.2017.4, PMID: 28191884HowellK. J.KraiczyJ.NayakK. M.GasparettoM.RossA.LeeC.. (2018). DNA methylation and transcription patterns in intestinal epithelial cells from pediatric patients with inflammatory bowel diseases differentiate disease subtypes and associate with outcome. Gastroenterology154, 585–598. doi: 10.1053/j.gastro.2017.10.007, PMID: 29031501Inflammatory Bowel Disease Group, Chinese Society of Gastroenterology, Chinese Medical Association (2021). Chinese consensus on diagnosis and treatment in inflammatory bowel disease (2018, Beijing). J. Dig. Dis.22, 298–317. doi: 10.1111/1751-2980.12994, PMID: 33905603JädertC.PhillipsonM.HolmL.LundbergJ. O.BorniquelS. (2014). Preventive and therapeutic effects of nitrite supplementation in experimental inflammatory bowel disease. Redox Biol.2, 73–81. doi: 10.1016/j.redox.2013.12.012, PMID: 24494186KaserA.BlumbergR. S. (2009). Endoplasmic reticulum stress in the intestinal epithelium and inflammatory bowel disease. Semin. Immunol.21, 156–163. doi: 10.1016/j.smim.2009.01.001, PMID: 19237300KaserA.BlumbergR. S. (2011). Autophagy, microbial sensing, endoplasmic reticulum stress, and epithelial function in inflammatory bowel disease. Gastroenterology140, 1738–1747.e2. doi: 10.1053/j.gastro.2011.02.048, PMID: 21530740KhorB.GardetA.XavierR. J. (2011). Genetics and pathogenesis of inflammatory bowel disease. Nature474, 307–317. doi: 10.1038/nature10209, PMID: 21677747KimD.LangmeadB.SalzbergS. L. (2015). HISAT: a fast spliced aligner with low memory requirements. Nat. Methods12, 357–360. doi: 10.1038/nmeth.3317, PMID: 25751142KoC.-W.QuJ.BlackD. D.TsoP. (2020). Regulation of intestinal lipid metabolism: current concepts and relevance to disease. Nat. Rev. Gastroenterol. Hepatol.17, 169–183. doi: 10.1038/s41575-019-0250-7, PMID: 32015520KöktenT.GibotS.LepageP.D’AlessioS.HablotJ.NdiayeN. C.. (2018). TREM-1 inhibition restores impaired autophagy activity and reduces colitis in mice. J. Crohns Colitis12, 230–244. doi: 10.1093/ecco-jcc/jjx129, PMID: 28961797LarabiA.BarnichN.NguyenH. T. T. (2020). New insights into the interplay between autophagy, gut microbiota and inflammatory responses in IBD. Autophagy16, 38–51. doi: 10.1080/15548627.2019.1635384, PMID: 31286804LavelleA.SokolH. (2020). Gut microbiota-derived metabolites as key actors in inflammatory bowel disease. Nat. Rev. Gastroenterol. Hepatol.17, 223–237. doi: 10.1038/s41575-019-0258-z, PMID: 32076145LeeJ. M.WagnerM.XiaoR.KimK. H.FengD.LazarM. A.. (2014). Nutrient-sensing nuclear receptors coordinate autophagy. Nature516, 112–115. doi: 10.1038/nature13961, PMID: 25383539LiB.DeweyC. N. (2011). RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics12:323. doi: 10.1186/1471-2105-12-323, PMID: 21816040LiuY.YinF.HuangL.TengH.ShenT.QinH. (2021). Long-term and continuous administration of during remission effectively maintains the remission of inflammatory bowel disease by protecting intestinal integrity, regulating epithelial proliferation, and reshaping microbial structure and function. Food Funct.12, 2201–2210. doi: 10.1039/d0fo02786c, PMID: 33595001Lloyd-PriceJ.ArzeC.AnanthakrishnanA. N.SchirmerM.Avila-PachecoJ.PoonT. W.. (2019). Multi-omics of the gut microbial ecosystem in inflammatory bowel diseases. Nature569, 655–662. doi: 10.1038/s41586-019-1237-9, PMID: 31142855LongS. L.GahanC. G. M.JoyceS. A. (2017). Interactions between gut bacteria and bile in health and disease. Mol. Asp. Med.56, 54–65. doi: 10.1016/j.mam.2017.06.002, PMID: 28602676LoucaS.ParfreyL. W.DoebeliM. (2016). Decoupling function and taxonomy in the global ocean microbiome. Science353, 1272–1277. doi: 10.1126/science.aaf4507, PMID: 27634532LoveM. I.HuberW.AndersS. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol.15:550. doi: 10.1186/s13059-014-0550-8, PMID: 25516281MagočT.SalzbergS. L. (2011). FLASH: fast length adjustment of short reads to improve genome assemblies. Bioinformatics27, 2957–2963. doi: 10.1093/bioinformatics/btr507, PMID: 21903629MagroF.GionchettiP.EliakimR.ArdizzoneS.ArmuzziA.Barreiro-de AcostaM.. (2017). Third European evidence-based consensus on diagnosis and management of Ulcerative Colitis. Part 1: definitions, diagnosis, extra-intestinal manifestations, pregnancy, cancer surveillance, surgery, and Ileo-anal pouch disorders. J. Crohns Colitis11, 649–670. doi: 10.1093/ecco-jcc/jjx008, PMID: 28158501ManichanhC.Rigottier-GoisL.BonnaudE.GlouxK.PelletierE.FrangeulL.. (2006). Reduced diversity of faecal microbiota in Crohn’s disease revealed by a metagenomic approach. Gut55, 205–211. doi: 10.1136/gut.2005.073817, PMID: 16188921Matsuzawa-IshimotoY.ShonoY.GomezL. E.Hubbard-LuceyV. M.CammerM.NeilJ.. (2017). Autophagy protein ATG16L1 prevents necroptosis in the intestinal epithelium. J. Exp. Med.214, 3687–3705. doi: 10.1084/jem.20170558, PMID: 29089374Mirsepasi-LauridsenH. C.VallanceB. A.KrogfeltK. A.PetersenA. M. (2019). Pathobionts associated with inflammatory bowel disease. Clin. Microbiol. Rev.32, e00060–e00018. doi: 10.1128/CMR.00060-18, PMID: 30700431MorettiJ.RoyS.BozecD.MartinezJ.ChapmanJ. R.UeberheideB.. (2017). STING senses microbial viability to orchestrate stress-mediated autophagy of the endoplasmic reticulum. Cell171, 809–823.e13. doi: 10.1016/j.cell.2017.09.034, PMID: 29056340MorganX. C.KabakchievB.WaldronL.TylerA. D.TickleT. L.MilgromR.. (2015). Associations between host gene expression, the mucosal microbiome, and clinical outcome in the pelvic pouch of patients with inflammatory bowel disease. Genome Biol.16:67. doi: 10.1186/s13059-015-0637-x, PMID: 25887922NakagawaI.AmanoA.MizushimaN.YamamotoA.YamaguchiH.KamimotoT.. (2004). Autophagy defends cells against invading group A Streptococcus. Science306, 1037–1040. doi: 10.1126/science.1103966, PMID: 15528445NguyenH. T. T.LapaquetteP.BringerM.-A.Darfeuille-MichaudA. (2013). Autophagy and Crohn’s disease. J. Innate Immun.5, 434–443. doi: 10.1159/000345129, PMID: 23328432NighotP. K.HuC.-A. A.MaT. Y. (2015). Autophagy enhances intestinal epithelial tight junction barrier function by targeting claudin-2 protein degradation. J. Biol. Chem.290, 7234–7246. doi: 10.1074/jbc.M114.597492, PMID: 25616664OttS. J.MusfeldtM.WenderothD. F.HampeJ.BrantO.FölschU. R.. (2004). Reduction in diversity of the colonic mucosa associated bacterial microflora in patients with active inflammatory bowel disease. Gut53, 685–693. doi: 10.1136/gut.2003.025403, PMID: 15082587PascalV.PozueloM.BorruelN.CasellasF.CamposD.SantiagoA.. (2017). A microbial signature for Crohn’s disease. Gut66, 813–822. doi: 10.1136/gutjnl-2016-313235, PMID: 28179361PerteaM.PerteaG. M.AntonescuC. M.ChangT. C.MendellJ. T.SalzbergS. L. (2015). StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol.33, 290–295. doi: 10.1038/nbt.3122, PMID: 25690850PetersL. A.PerrigoueJ.MorthaA.IugaA.SongW.-M.NeimanE. M.. (2017). A functional genomics predictive network model identifies regulators of inflammatory bowel disease. Nat. Genet.49, 1437–1449. doi: 10.1038/ng.3947, PMID: 28892060PlichtaD. R.GrahamD. B.SubramanianS.XavierR. J. (2019). Therapeutic opportunities in inflammatory bowel disease: mechanistic dissection of host-microbiome relationships. Cell178, 1041–1056. doi: 10.1016/j.cell.2019.07.045, PMID: 31442399QuinnR. A.MelnikA. V.VrbanacA.FuT.PatrasK. A.ChristyM. P.. (2020). Global chemical effects of the microbiome include new bile-acid conjugations. Nature579, 123–129. doi: 10.1038/s41586-020-2047-9, PMID: 32103176R Core Team (2013). R: A Language and Environment for Statistical Computing.RobinsonM. D.McCarthyD. J.SmythG. K. (2010). edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics26, 139–140. doi: 10.1093/bioinformatics/btp616, PMID: 19910308RyanF. J.AhernA. M.FitzgeraldR. S.Laserna-MendietaE. J.PowerE. M.ClooneyA. G.. (2020). Colonic microbiota is associated with inflammation and host epigenomic alterations in inflammatory bowel disease. Nat. Commun.11:1512. doi: 10.1038/s41467-020-15342-5, PMID: 32251296SaijoF.MilsomA. B.BryanN. S.BauerS. M.VowinkelT.IvanovicM.. (2010). On the dynamics of nitrite, nitrate and other biomarkers of nitric oxide production in inflammatory bowel disease. Nitric Oxide22, 155–167. doi: 10.1016/j.niox.2009.11.009, PMID: 20005300ShenS.ParkJ. W.LuZ. X.LinL.HenryM. D.WuY. N.. (2014). rMATS: robust and flexible detection of differential alternative splicing from replicate RNA-Seq data. Proc. Natl. Acad. Sci. U. S. A.111, E5593–E5601. doi: 10.1073/pnas.1419161111, PMID: 25480548ShindeT.PereraA. P.VemuriR.GondaliaS. V.BealeD. J.KarpeA. V.. (2020). Synbiotic supplementation with prebiotic green banana resistant starch and probiotic Bacillus coagulans spores ameliorates gut inflammation in mouse model of inflammatory bowel diseases. Eur. J. Nutr.59, 3669–3689. doi: 10.1007/s00394-020-02200-9, PMID: 32067099SinhaS. R.HaileselassieY.NguyenL. P.TropiniC.WangM.BeckerL. S.. (2020). Dysbiosis-induced secondary bile acid deficiency promotes intestinal inflammation. Cell Host Microbe27, 659–670.e5. doi: 10.1016/j.chom.2020.01.021, PMID: 32101703SoltysK.StuchlikovaM.HlavatyT.GaalovaB.BudisJ.GazdaricaJ.. (2020). Seasonal changes of circulating 25-hydroxyvitamin D correlate with the lower gut microbiome composition in inflammatory bowel disease patients. Sci. Rep.10:6024. doi: 10.1038/s41598-020-62811-4, PMID: 32265456StackebrandtE.GoebelB. M. (1994). Taxonomic note: a place for DNA-DNA reassociation and 16S rRNA sequence analysis in the present species definition in bacteriology. Int. J. Syst. Evol. Microbiol.44, 846–849. doi: 10.1099/00207713-44-4-846StengelS. T.FazioA.LipinskiS.JahnM. T.AdenK.ItoG.. (2020). Activating transcription factor 6 mediates inflammatory signals in intestinal epithelial cells upon endoplasmic reticulum stress. Gastroenterology159, 1357–1374.e10. doi: 10.1053/j.gastro.2020.06.088, PMID: 32673694SudhakarP.JacominA.-C.HautefortI.SamavedamS.FatemianK.AriE.. (2019). Targeted interplay between bacterial pathogens and host autophagy. Autophagy15, 1620–1633. doi: 10.1080/15548627.2019.1590519, PMID: 30909843WangQ.GarrityG. M.TiedjeJ. M.ColeJ. R. (2007). Naive Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Appl. Environ. Microbiol.73, 5261–5267. doi: 10.1128/aem.00062-07, PMID: 17586664WardT.LarsonJ.MeulemansJ.HillmannB.LynchJ.SidiropoulosD.. (2017). BugBase predicts organism-level microbiome phenotypes. biorxiv [Preprint]. 133462. doi: 10.1101/133462WuY.WangB.XuH.TangL.LiY.GongL.. (2019). Probiotic attenuates oxidative stress- induced intestinal injury via p38-mediated autophagy. Front. Microbiol.10:2185. doi: 10.3389/fmicb.2019.02185, PMID: 31632359XieC.MaoX.HuangJ.DingY.WuJ.DongS.. (2011). KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res.39, W316–W322. doi: 10.1093/nar/gkr483, PMID: 21715386ZhangX.TongY.LyuX.WangJ.WangY.YangR. (2020). Prevention and alleviation of dextran sulfate sodium salt-induced inflammatory bowel disease in mice with -fermented milk inhibition of the inflammatory responses and regulation of the intestinal flora. Front. Microbiol.11:622354. doi: 10.3389/fmicb.2020.622354, PMID: 33519783Glossary
AKR1C2
Aldo-keto reductase family 1 member C2
AKT2
AKT serine/threonine kinase 2
ATF6
Activating transcription factor 6
ATG16L1
Autophagy related 16 like 1
ATG4A
Autophagy related 4A cysteine peptidase
ATG5
Autophagy related 5
ATG9A
Autophagy related 9A
BAG6
BAG cochaperone 6
BBC3
BCL2 binding component 3
BCL2L11
BCL2 like 11
BECN1
Beclin 1
BNIP3P
BCL2 interacting protein 3 pseudogene 1
CASP17P
Caspase 17, pseudogene
CD
Crohn’s disease
CHAC1
ChaC glutathione specific gamma-glutamylcyclotransferase 1