Dynamics of the gut microbiota in rats after hypobaric hypoxia exposure

Animal experiments

For the animal experiments, 10-week-old Sprague-Dawley (SD) male rats (n = 12) weighing 290 to 300 g were supplied by the Animal Center of the General Hospital of Chinese People’s Liberation Army (PLAGH). Animal experiments were performed strictly in accordance with the Guide for the Care and Use of Laboratory Animals published by the PLAGH. Rats were randomly assigned into the following two groups: the normal control group (NC, n = 5) and the hypobaric hypoxia group (HH, n = 7). Rats in the NC group were fed in normoxic conditions at an altitude of 43.5 m in Beijing and the samples from the NC group were collected only once. Rats in the HH group were fed in the hypobaric chamber (FLYDWC50-1C low pressure hypoxic experimental cabin produced by Guizhou Fenglei Air Ordance Co., Ltd. Guizhou, China) for 5 weeks (Inside group) and then returned to the normoxic conditions fed for 3 weeks (Outside group). The conditions of the hypobaric chamber were set to emulate an altitude of 5,500 m (50 KPa, 380 mmHg). In the Inside group, fecal samples at 1st day (In1d), 3rd day (In3d), 5th day (In5d), 1st week (In1w), 2nd week (In2w), 3rd week (In3w), 4th week (In4w) and 5th week (In5w) were collected sequentially. In the Outside group, fecal samples at 1st day (Out1d), 3rd day (Out3d), 5th day (Out5d), 1st week (Out1w), 2nd week (Out2w), 3rd week (Out3w) were collected sequentially. In the Inside group, fecal samples of three rats were removed due to contamination. In both groups, two rats were housed in a cage with a height of 30 cm and the floor of the cage covered with an adequate depth of dust free Aspen wood chip. To avoid coprophagy among rats, we isolated two rats housed in the same cage using wooden board. Playpens were provided in the cage for rats to exercise and socialize. The housed temperature and humidity were controlled at 25 ± 2 °C and 50 ± 10% with a 12 h light/12 h dark cycle. All rats have ad libitum access to standard pellet diet and water. Animals did need to be euthanized until the end of the experiment. If the rats are obviously unwell, the experiment should be terminated. We determined whether animals were unwell by observing behaviors such as silence, lethargy, inactivity, loss of appetite, and refusal to eat. After the experiment, all the animals were euthanized by intraperitoneal injection of barbiturates. All animal procedures and experimental protocols were reviewed and approved by the Animal Ethics Committee of PLAGH (2017-X13-05).

Sample collection and DNA extraction

All fecal samples collected were frozen immediately at −80 °C within 20 min and transported to the laboratory with dry ice to extract bacteria DNA. DNA extraction was conducted at Novogene Bioinformatics Technology Co., Ltd. according to the standard protocols of Tiangen kits. The magnetic soil and stool DNA kit (TIANGEN®, Beijing, China) was used for DNA extraction. Special grinding beads were applied in this kit to effectively process the lysis of various complex components in fecal samples. The 0.25−0.5 g fecal samples were added to two mL centrifuge tube, and then 500 µL buffer SA, 100 µL buffer SC and 0.25 g grinding beads were added and mixed evenly with vortex. The mix was heated at 70 °C for 15 min to improve lysis efficiency, and then centrifuged at 12,000 rpm (∼13400 × g) for 1 min and the supernatant was transferred to a new two mL tube. 200 µL buffer SH was added and mixed with vortex for 5 s, and stood at 4 °C for 10 min. The mix was centrifuged at 12,000 rpm (∼13400 × g) for 3 min, the supernatant was transferred to a new two mL centrifuge tube, then 500 buffer GFA was added and mixed evenly. 10 µL of MagAttract Suspension G was added and mixed for 5 min. The centrifuge tube was placed on the magnetic stand for 30 s. After the beads were completely attached, the liquid in the tube was removed. Then, 700 µL of buffer RD was added and mixed for 5 min. The operation of placing the tube on the magnetic stand was repeated. Then, 700 µL of buffer PWD was added and mixed for 3 min and repeated once. After removing the liquid, the tube was placed on the magnetic stand and dried at the room temperature. Then, 80 µL of buffer TB was added and incubated at 56 °C for 5 min, and mixed three times, four times each time. After the beads were completely attached on the magnetic stand, the DNA solution was transferred and stored at −80 °C for subsequent sequencing analysis.

Metagenomic sequencing and DNA library construction

All samples were paired-end sequenced on the Illumina platform at Novogene Bioinformatics Technology Co., Ltd. with insert size 350 bp and read length 150 bp. Preprocessing the raw data obtained from the Illumina HiSeq sequencing platform using Readfq (v8, https://github.com/cjfields/readfq) was conducted for quality control to acquire the clean data. Following quality control, removing reads containing low-quality bases (quality value ≤ 38) over 40 bp, ‘N’ bases over 10 bp and overlaps with the adaptor more than 15 bp. Also, the reads aligned to the human genome with Bowtie2 were removed (v2.2.4) (Langmead & Salzberg, 2012). The remaining clean reads were assembled by SOAP denovo (v2.04) (Luo et al., 2012) with several K-mer values (from 27 to 87), and the assembled results are interrupted from the ‘N’ junction to obtain Scaftigs without ‘N’ (Mende et al., 2012; Nielsen et al., 2014). All unmapped reads were combined for mixed assembly using the same methods and parameters as forward step. Scaftigs longer than 500 bp were retained for gene prediction using MetaGeneMark (v2.10). The unique initial gene catalogue was obtained by CD-HIT (v4.5.8) (Fu et al., 2012) with 95% identity and 90% minimum coverage. The abundance of each gene in each sample was calculated based on the number of aligned reads using Bowtie2.

Taxonomy annotation and abundance profiling

In total, 3,096,427 de-redundant genes were used to assign taxonomic groups. Taxonomic assignments for genes were obtained by aligning the genes to integrated NR database (v 2018-01-02) of NCBI including Bacteria, Fungi, Archaea and Viruses using DIAMOND (v 0.9.9) (Buchfink, Xie & Huson, 2015) with the parameter setting are blastp, -e 1e−5. Taxonomic annotation of each gene was determined by the lowest common ancestor-based algorithm (LCA). The results of which the e value ≤ the smallest e value ×10 were chosen to take the LCA algorithm which was applied to system classification of MEGAN software to make sure the species annotation information of sequences (Oh et al., 2014). The abundance of a taxonomic level in each sample was calculated by summing the abundances of all genes annotated to that level and the gene number of a taxonomic level in each sample equals to the number of genes whose abundance was not zero.

KEGG annotation and abundance profiling

The functional assignments of genes were acquired by aligning genes to Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2014) database using DIAMOND, and the highest-scoring annotated hits were retained for subsequent analysis. The abundance of each functional hierarchy was calculated by summing the abundance of genes annotated to that level.

Microbial diversity and statistical analysis

The Shannon and Simpson diversity indexes (α-diversity) were calculated based on normalized species profiles by R package vegan (v2.5-6). β-diversity was analyzed by principal coordinate analysis (PCoA) by R package ape (v5.3) based on the Bray-Curtis distance. Dissimilarity analysis among groups was conducted by analysis of similarities (ANOSIM) by R package vegan (v2.5-6). Differential species or KEGG modules between Inside and Outside groups were identified using the MaAsLin2 by R package Maasline2 (v1.0.0). To remove the individual effect in longitudinal study, rat id was used as the random effect in differential analysis. Fold changes were the ratio of the average relative abundance of species or KEGG modules in Outside group and Inside group. The Benjamin-Hochberg false discovery rate (FDR) correction was performed for multiple hypothesis testing to compute FDR corrected p-values (q-values). Differential species or KEGG modules with q-value <0.05 were retained for subsequent analysis. Trend lines fitted using LOESS smoothing and the curve was plotted with a 95% pointwise confidence interval band in gray. A paired t test was performed, plotted and only one set of paired t-test was conducted at a time using GraphPad Prism version 8.3.1 (GraphPad Software, San Diego, California USA, www.graphpad.com).

Correlation analysis

Correlation analysis was conducted for differential species and KEGG modules by R package psych (v2.0.7) based on Spearman’s rank correlation. Correlations with an absolute value of coefficient (r) >0.7 and q-value <0.05 were retained for further analysis and visualized by Cytoscape 3.7.2 (Shannon et al., 2003). Heat map was plotted by R package pheatmap (v1.0.12).

Alterations of microbial community diversity and structure

The experimental design and analysis processes of this study are presented in Fig. 1. To study the alterations of within-sample diversity (α-diversity) over time under hypobaric hypoxia conditions, Shannon and Simpson diversity index of all samples were calculated, respectively. Consistently, the Shannon and Simpson indexes both displayed similar trends over time (Figs. 2A and 2B). After entering the hypobaric chamber, we found that α-diversity had a significant fluctuation within the first 5 days after entering the chamber and then increased gradually and decreased after four weeks. To further statistically explore the fluctuation of diversity within the first 5 days and over time after entering or leaving the chamber, paired t test was performed between groups within the first 5 days and groups after that. Paired t test result showed that there is a significant difference in Shannon index between In5d and In3w (p = 0.02, Fig. 2E). After leaving the chamber, a significant fluctuation of α-diversity also occurred within the first 5 days and then increased gradually. Paired t test results also showed a significant difference in Shannon index between Out5d and Out3d (p = 0.002, Fig. 2F). Also, there was a significant difference between Out3w and Out3d (p = 0.001, Fig. 2G). The individual dynamics of α-diversity was presented in Figs. S1A and S1B. The α-diversity results indicated that the gut microbiota of rats responded strongly to hypobaric hypoxia or normoxia within 5 days. Furthermore, the overall community structure (β-diversity) across groups were examined by principal coordinates analysis (PCoA) based on Bray-Curtis distance (Fig. 2C). Significant differences were observed in the overall community structure (ANOSIM test, R = 0.445, p < 0.001), which indicated that there was a clear separation of microbial structure over time with the change of oxygen concentration and atmospheric pressure. The first principal coordinate in the Bray-Curtis PCoA showed the significant fluctuations of gut microbiota over time from entering the hyperbaric chamber to returning to the normoxic conditions (Fig. 2D). Paired t test showed a significant difference in the first principal coordinate between In1d and In5d (p = 0.02, Fig. 2H). The individual dynamics of the first coordinate in the Bray-Curtis PCoA was presented in Fig S1.

Altered composition of gut microbiota over time

To explore the composition of gut microbiota over time, the top 10 abundant phyla and genera were screened out and visualized, respectively (Figs. 3A and 3B). In the phylum level, we observed that Firmicutes was more abundant than Bacteroidetes in the first three days after entering the hypobaric chamber, which was consistent with NC group. Paired t test result showed that the abundance of Bacteroidetes in In5d was significantly more abundant compared with In1d (p = 0.03, Fig. 3E). What’s more, from the second week to the fourth week, Firmicutes was more abundant than Bacteroidetes like initial state. Paired t test results showed that the abundance of Bacteroidetes in In5w was significantly more abundant compared with In3w (p = 0.04), while Firmicutes was reversed (p = 0.002, Fig. 3E). However, in the fifth week, the abundance of Bacteroidetes surpassed Firmicutes again. After leaving the chamber, the abundance of Bacteroidetes was the highest in the first three days. From the fifth day, the abundance of Firmicutes and Bacteroidetes reversed again, which was consistent with starting state and NC group. Paired t test result showed that the abundance of Bacteroidetes in Out3d was significantly more abundant compared with Out5d (p = 0.01), while Firmicutes was reversed (p = 0.01, Fig. 3E). The individual dynamics of relative abundance of Firmicutes and Bacteroidetes were displayed in Figs. S2A and S2B. These results indicated that normoxic conditions recovered the gut microbiota that destroyed by hypobaric hypoxia. Our results also revealed that Firmicutes was the most abundant phylum in normoxic conditions, and Bacteroidetes might was more adaptable to hypobaric hypoxia. Also, after entering or leaving the hypobaric chamber, the abundance of top 2 phyla both reversed on the fifth day, which implied that day 5 possibly was the time point when the composition of gut microbiota changed radically. In the genus level, we observed that the abundance of Prevotella was significantly increased after entering the hypobaric chamber and Lactobacillus was significantly enriched after leaving the chamber (Fig. 3B). The individual dynamics of relative abundance of top genera were displayed in Fig. S3. Paired t test result showed that the abundance of Prevotella was significantly increased In5d compared to In1d (p = 0.0005) and significantly decreased in Out5d compared to Out3d (p = 0.01, Fig. 3F). The abundance of Lactobacillus in Out5d was significantly increased compared with Out1d (p = 0.02, Fig. 3F). Previous study reported that Tibetan living at high altitude have a higher abundance of Prevotella whereas Bacteroides is more abundant in the Han living at the plain (Li et al., 2016). It implied that the Prevotella might be more adaptable to the high altitude with hypobaric hypoxia. Our results further prove the process of the abundance of Prevotella rising to the most abundant genus after exposure to the hypobaric hypoxia conditions. Although the abundance of Prevotella was decreased after leaving the hypobaric chamber five days later, it was still higher than Bacteroides. Previous study found that people living at high altitude maintained microbial signatures even after coming back to the plain three months later (Jia et al., 2020). Therefore, we speculated that it would take longer time to get back to the starting state.

To further investigate the relationships between dominant taxa, and their contributions to the separation of groups, loess regression was applied to points using the relative abundance and PCoA1 coordinates. We found that Firmicutes and Bacteroidetes showed an inverse trend (Fig. 3C). What’s more, the dominant genera, Prevotella and Bacteroides, had a similar relationship (Fig. 3D). These results further indicated that the competition between dominant taxa and the separation along the first principal coordinate was driven by a continuous variable that was a trade-off between Prevotella and Bacteroides (Gorvitovskaia, Holmes & Huse, 2016; Smits et al., 2017).

Differential species between Inside and Outside groups

To identify the species associated with hypobaric hypoxia, differential abundance analysis was performed between Inside and Outside groups using MaAsLin2. In total, 444 differential species were identified (q-value <0.05), in which 276 species were increased in Outside group and 168 species were increased in Inside group (Table S1). Differential species with top 20 significance were visualized with heat map (Fig. 4A).

In the top 20 significance species, 19 species were significantly enriched in Outside group. The abundance of Prevotella loescheii, which can ferment sugars to produce succinate and acetate, was the most significant species increased in Outside group. Also, Butyrivibrio, Pseudobutyrivibrio, Porphyromonas endodontalis and Anaerophaga thermohalophila, producing short-chain fatty acids (SCFAs) (Denger et al., 2002; Shah & Collins, 1988), were significantly enriched in the Outside group. Probiotics like Bifidobacterium and Lactococcus were significantly increased in Outside. Probiotics usually offer a large range of health benefits and have been indicated as a marker for health, thus their decreased abundance in Inside may be associated with some diseases like diarrhea induced by hypoxia (Gamah et al., 2021; Zheng et al., 2009). Archangium gephyra producing argyrin A, as a potent antitumoral drug (Nickeleit et al., 2008), was more abundant in Outside than that of Inside. Veillonella atypica, which improves run time via its metabolic conversion of exercise-induced lactate into propionate (Scheiman et al., 2019), was significantly increased in Outside group. In addition, some aerobes or facultative anaerobes like Mycoplasma, Enterococcus and Paenibacillus were significantly increased in Outside, while anaerobes like Desulfovibrio and Alistipes were decreased. As a result of their oxygen dependence, their increased or decreased abundance might be associated with the oxygen concentrations of Outside and Inside. The abundance of Bacteroides vulgatus and Bacteroides dorei, reducing gut microbial lipopolysaccharide production and inhibiting atherosclerosis (Yoshida et al., 2018), was significantly enriched in Inside group. Oxalobacter formigenes, regulating host oxalate homeostasis (Liu et al., 2017), was significantly increased in Inside group. Bacteroides fragilis and Bacteroides caccae, both are obligate anaerobes, were significantly enriched in Inside group. Bacteroides massiliensis enriched in obese children (Hollister et al., 2018), was also increased in Inside group.

Microbial function analysis

To investigate the alterations of microbial function between Inside and Outside groups, differential analysis was performed for KEGG modules using MaAsLin2. In total, 96 KEGG modules with significant difference from 530 annotated modules were investigated (q-value <0.05, Table S2). Top 20 significantly differential modules were displayed with heat map (Fig. 4B). We found 32 modules were increased in Outside group, and 64 modules were decreased. In the top 20 differential modules, ‘beta-Carotene biosynthesis’ (M00097), and ‘Iron complex transport system’ (M00240) were significantly increased in Inside group. Pathways involved in biosynthesis of amino acids such as ‘Proline biosynthesis’ (M00015), ‘Histidine biosynthesis’ (M00026), ‘Lysine biosynthesis’ (M00016) and ‘Isoleucine biosynthesis’ (M00535) were significantly enriched in Inside group. ‘Two-component regulatory system’ (M00499, M00434, M00452, M00453, M00512), a basic stimulus–response coupling mechanism to sense and respond to changes of environmental conditions (Stock, Robinson & Goudreau, 2000), were significantly increased in Inside group. ‘Methanogenesis’ (M00567), an anaerobic respiration that generates methane as the final product of metabolism (Lyu et al., 2018), and ‘Fatty acid biosynthesis’ (M00083) were significantly increased in Inside group. In contrast, ‘cAMP signaling’ (M00695), ‘Reductive pentose phosphate cycle’ (M00166) and ‘Biotin biosynthesis’ (M00573) were significantly increased in Outside group.

To further investigate the correlation between top 20 significant species and microbial modules, Spearman rank correlation was conducted (Fig. 4C, Table S3). Correlation results showed that there are strongly significant correlations between 19 species and 12 modules (—r—>0.7, q-value <0.05). We found that many species enriched in Outside group were negatively correlated with ‘Proline biosynthesis’ (M00015) and ‘Histidine biosynthesis’ (M00026). However, ‘Biotin biosynthesis’ (M00573), ‘C1-unit interconversion’ (M00141) and ‘Pyrimidine ribonucleotide biosynthesis’ (M00052) were positively correlated with almost all species. Correlation results revealed that differential species influenced by hypobaric hypoxia conditions significantly facilitate the alterations of microbial functions.