Intermittent Hypoxia and Hypercapnia, a Hallmark of Obstructive Sleep Apnea, Alters the Gut Microbiome and Metabolome

1Obstructive sleep apnea (OSA) afflicts nearly 12% of the adult population in the United States, with an annual cost burden of nearly $149.6 billion, per a recent study commissioned by the American Academy of Sleep Medicine. Timely diagnosis and treatment of OSA improves not only sleep and cognitive function but also management of comorbid cardiometabolic diseases (CMDs). Therefore, identifying the mechanisms underlying the downstream consequences of OSA would aide in development of effective treatment modalities, reducing overall health care utilization.

2OSA is marked by changes in oxygen and carbon dioxide-inspired concentrations, which impact the gut microbial community. Since the gut microbiota plays a key role in the metabolism of dietary precursors, including lipids, cholesterol, and choline, it impacts the cardiometabolic health of the host. Gut dysbiosis has already been linked to an array of cardiovascular and metabolic disorders, such as hypertension, type 2 diabetes, hepatic steatosis, and atherosclerosis. Additionally, previous work has identified specific gut bacteria to be significantly correlated with plasma cholesterol and apolipoprotein levels. Thus, probing this commensal ecosystem may provide a valuable avenue of investigation to understand the mechanisms underlying the cardiovascular consequences of OSA. In this study, we investigated the taxonomic and molecular alterations in the gut microbiome that potentially mediate the interplay between OSA and related CMDs.

3We used atherosclerosis-prone (Ldlr−/−) adult mice fed a high-fat diet (HFD) enriched in cholesterol and milk fat (resembling Western dietary practices) to evaluate atherosclerosis risk in OSA. We previously demonstrated that IHH increases atherosclerosis plaque formation in this model. As episodic hypoxia and hypercapnia mimic the changes in blood gases that occur in OSA-driven downstream consequences, these mice were exposed to IHH (treatment group; n = 8) or air (control group; n = 8) and examined longitudinally for 6 weeks (see Fig. S1 in the supplemental material). Fecal samples, representative of the gut ecosystem, were collected at baseline and twice each week thereafter, and the microbiome and metabolome were profiled using 16S rRNA amplicon sequencing and liquid chromatography-tandem mass spectrometry (LC-MS/MS)-based untargeted mass spectrometry, respectively. These data were processed to obtain relative abundances of microbial and molecular species per sample (referred to as “feature tables” henceforth), which were used for comparing mice mimicking humans with OSA and control mice.

4First, we performed principal-coordinate analysis (PCoA) on the microbiome and metabolome feature tables to identify major factors driving the clustering of samples. Figure 1 shows the PCoA results plotted against time to visualize the dynamics of clustering based on the gut microbiome (unweighted UniFrac distances are shown in Fig. 1a) and the metabolome (Gower distances are shown in Fig. 1b and c) as the duration of IHH exposure increases. Here, the first fecal sample represents the baseline gut composition before animals were switched to an HFD. There is a rapid shift in both microbial and molecular composition due to HFD alone, consistent with similar previous findings. Moreover, starting from a highly congruent gut composition, IHH-exposed mice significantly diverge from controls with increasing exposure duration (a permutational multivariate analysis of variance [PERMANOVA] test was performed per time point [Table S1]). This demonstrates that prolonged IHH exposure (analogous to chronic OSA) cumulatively perturbs the gut microbiome and metabolome. We tested the relationship between the two omics data sets by superimposing the principal coordinates computed from microbiome and metabolome data (procrustes analysis) (Fig. 1d and e). The ordination spaces are correlated (Mantel test r statistic = 0.36, P < 0.001), and changes in the metabolome and microbiome of samples within the treatment groups over time are proportional, suggesting microbe-dependent changes in intestinal metabolism with chronic OSA.

5We then tested for specific microbes and metabolites that changed with OSA. More than 80 (of ~730) microbial features differed significantly between the IHH-exposed group and controls (by permutation test with discrete false-discovery rate [FDR] correction). Figure S2a presents a global overview of these changes in gut microbiota per sample (sorted by duration of treatment). Table S2 provides a list of these differentially represented bacteria that potentially contribute to alterations in gut metabolism due to IHH. Figure 2a to f display trends in relative abundances of bacteria showing the largest differences, which belong to the Mogibacteriaceae (family), Oscillospira (genus), Lachnospiraceae (family), and Clostridiaceae (family). Previous studies have consistently associated these taxonomic groups with metabolic and inflammatory disturbances in the host, which suggests that related mechanisms might be at play in driving the consequences of hypoxic and hypercapnic stress.

6Using the same statistical approach, we found that more than 380 out of ~1,700 molecular species (MS1 spectral features) differed significantly in relative abundance in animals exposed to IHH. Figure S2b provides a global representation of these differentially abundant molecules in samples belonging to treatment and control groups and sorted by treatment duration (Table S3 provides a comprehensive list of these molecules). To gain insight into the structures of these differentially abundant metabolites, we performed molecular networking using Global Natural Products Social Molecular Networking (GNPS). The molecular network is constructed using a cosine similarity measure between tandem mass spectral data and then visualized using cytoscape (Fig. S3). Each node in the network, which represents a consensus MS/MS spectrum, was searched against public libraries in this study. In total, we annotated about 400 molecular compounds in the GNPS database, including bile acids (BAs), fatty acids, and phytoestrogens. Additionally, all key compounds discussed in this work were defined to the highest level of annotation per the metabolomic standards initiative using commercial standards (Fig. S4, Fig. S5, and Table S4).

7Interestingly, the top differentially abundant features detected between IHH and control mice included molecules known to depend on gut microbes for their production. Below, we discuss some of these metabolites and their implications with respect to consequences of OSA.

8We observed significant alterations in BAs between IHH-exposed mice and control groups (Table S3). Figure 2g to n display these trends in primary (Fig. 2g to i) and secondary (Fig. 2j to l) bile acids with increasing IHH exposure duration. Primary BAs are amphipathic molecules synthesized in the liver from cholesterol. These are conjugated to glycine or taurine and released in the biliary tract. Together with other biliary components, these facilitate the emulsification and transportation of dietary fats, cholesterol, and fat-soluble vitamins. About 95% of the BAs are reabsorbed in the terminal ileum and recycled. The remaining 5% reach the colon and are deconjugated, dehydrogenated, and dehydroxylated by the intestinal bacteria to form secondary BAs. BAs, including microbially generated BAs, are potent signaling molecules that interact with the farnesoid X receptor (FXR) (expressed in the liver and intestine), which modulates BA synthesis by the liver. Perturbations in the gut microbial population disrupt normal signaling properties that regulate BA production and can profoundly alter the BA composition in the gut. A range of diseases, including cardiometabolic diseases, are characterized by aberrant BA profiles, and prolonged perturbations in the BA pool might also be a factor in mediating the consequences of OSA.

9The dietary hormones enterolactone (mammalian lignan) and enterodiol (oxidation product of enterolactone) were significantly elevated in the exposed mice compared to their levels in controls. Figure 2m and n show the trends in their abundances with increasing duration of IHH exposure. These molecules are phytoestrogens, i.e., plant-derived hormones that structurally mimic estrogen and are produced by intestinal microbiota upon bioconversion of dietary lignans. Owing to their affinity to estrogen receptors (producing estrogenic or/and antiestrogenic effects), they perturb many hormone-dependent systems in the body and have been linked to adverse metabolic, reproductive, and neurological outcomes. Sex-specific differences in OSA-diagnostic symptoms and risk factors suggest hormonal involvement. However, the contribution of microbes in maintaining hormonal homeostasis has not yet been investigated. Therefore, these findings motivate novel avenues of research for biomarkers and therapeutic targets to manage the metabolic consequences of OSA.

10In addition to detecting changes in bile acids and phytoestrogens, we detected differentially abundant fatty acid-related chemical families (Table S3). For example, we noted a significant reduction in a molecular feature matched to elaidic acid. Elaidic acid is an unsaturated fatty acid that increases plasma cholesteryl ester transfer protein (CETP) activity, which modulates systemic levels of low-density lipoprotein (LDL) and high-density lipoprotein (HDL) cholesterol. A decrease in elaidic acid in the IHH-exposed group suggests reduction in plasma CETP activity, a mechanism associated with adverse cardiovascular effects. Similarly, phytomonic, jasmonic, hexadecanoic, linoleic, and conjugated linoleic acids were also reduced in exposed mice compared to levels in controls. Of these, phytomonic acid and conjugated linoleic acid are known to be microbially produced, suggesting that changes in the microbiome contribute to these changes in metabolome.

11In summary, we demonstrate that IHH, a hallmark of OSA, changes the microbiota and the chemistry in the gut. We have highlighted changes in bile acids, phytoestrogens, and fatty acids under OSA-related conditions that could lead to CMDs. The present results reveal a previously unrecognized mechanistic link between OSA and gut microbes. It suggests that targeting gut microbiota and their metabolites may serve as a potential therapeutic approach for the treatment of cardiometabolic consequences of OSA patients.

12Atherosclerosis-prone 10-week-old male Ldlr−/− mice in the C57BL/6J background (stock number 002207; The Jackson Laboratory, Bar Harbor, ME) were used in this study. Ldlr deficiency was confirmed by PCR per the vendor’s instructions. All animal protocols were approved by the Animal Care Committee of the University of California, San Diego, and followed the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health.

13Starting at 10 weeks of age, male mice were provided with a high-fat diet (HFD) containing 1.25% cholesterol and 21% milk fat (4.5 kcal/g; TD96121; Harlan-Teklad, Madison, WI) for 6 weeks while being exposed to either IHH or room air.

14Intermittent hypoxia and hypercapnia (IHH) was maintained in a computer-controlled atmosphere chamber system (OxyCycler; Reming Bioinstruments, Redfield, NY) as previously described. IHH exposure was introduced to the mice in short periods (~4 min) of synchronized reduction of O2 (from 21% to 8%) and elevation of CO2 (from ~0.5% to 8%), separated by alternating periods (~4 min) of normoxia ([O2] = 21%) and normocapnia ([CO2] = ~0.5%), with 1- to 2-min ramp intervals for 10 h per day during the light cycle for 6 weeks. This treatment protocol mimics the severe clinical condition observed in obstructive sleep apnea patients. Mice on the same HFD but in room air were used as controls.

15As the experimental setup requires IHH-exposed mice in a controlled atmosphere chamber and controls in room air, we ensured that the effect of treatment is not confounded by the effect of distinct housing conditions. To do so, we used two cages per treatment group, and we compared the relative effect sizes of treatments and cages with redundancy analysis (RDA), which estimates the independent effect size of each covariate on microbiome composition variation based on unweighted UniFrac distance. The RDA results showed that treatment had a higher effect size than the cages, more specifically, that treatment contributed to 11.6% of the microbiome community variation, while cages had an independent effect size of around 9.8%; with respect to the metabolome, treatment contributed to 6.2% of the variation, while cages contributed to only about 0.7% of the variation.

16Prior to LC-MS/MS analysis, fecal samples were prepared using the following extraction procedures. For extraction, 500 µl of 50:50 methanol-H2O was added to all fecal samples (30 to 50 mg approximately) and vortexed. Fecal pellets in extraction solvent were placed in an ultrasonic bath and sonicated for 30 min to break apart the pellet and then allowed to incubate for an additional 30 min. Extracted samples were then centrifuged to separate insoluble material, and 450 µl of each liquid extract was subsequently transferred to a 96-well deep-well plate and dried completely using centrifugal evaporation (CentriVap centrifugal vacuum concentrator; Labconco, Kansas City, MO). The dried extracts were resuspended in 150 µl of methanol-H2O (1:1, vol/vol), with 1 µM amitriptyline included as an autosampler injection standard. After resuspension, the samples were transferred into 96-well plates and analyzed on a Vanquish ultrahigh-performance liquid chromatography (UPLC) system coupled to a Q Exactive orbital ion trap (Thermo Fisher Scientific, Bremen, Germany). For the chromatographic separation, a C18 core shell column (Kinetex column, 50 by 2 mm, 1.7-µm particle size, 100-Å pore size; Phenomenex, Torrance, CA) with a flow rate of 0.5 ml/min (solvent A, H2O-0.1% formic acid [FA]; solvent B, acetonitrile-0.1% FA) was used. After being injected, the samples were eluted from 0 to 0.5 min with a linear gradient of 5% solvent B, from 0.5 to 4 min with 5 to 50% solvent B, and from 4 to 5 min with 50 to 99% solvent B, followed by a 2-min washout phase at 99% solvent B and a 2-min reequilibration phase at 5% solvent B. For online MS/MS measurements, the flow was directed to a heated electrospray ionization (HESI) source. The ESI parameters were set to 35 liters/min for the sheath gas flow, 10 liter/min for the auxiliary gas flow, 2 liters/min for the sweep gas flow, and 400°C for the auxiliary gas temperature. The spray voltage was set to 3.5 kV, and the inlet capillary was set to 250°C. A 50-V S-lens radio frequency (RF) level was applied. Product ion spectra were recorded in data-dependent acquisition (DDA) mode. Both MS1 survey scans (m/z 150 to 1,500) and up to 5 MS/MS scans of the most abundant ions per duty cycle were measured with a resolution (R) of 17,500 with 1 microscan in positive mode. The maximum ion injection time was set to 100 ms. MS/MS precursor selection windows were set to m/z 3 with an m/z 0.5 offset. Normalized collision energy was stepwise increased from 20 to 30 to 40% with a z of 2 as the default charge state. MS/MS experiments were automatically triggered at the apex of a peak within 2 to 15 s from their first occurrence. Dynamic exclusion was set to 5 s.

17Feature detection was as follows. Thermo raw data sets were converted to m/z extensible markup language (mzXML) in centroid mode using MSConvert (part of ProteoWizard). All mzXML files were cropped with an m/z range of 75.00 to 1,000.00 Da. MS1-based feature detection and MS2-based molecular networking was performed using the GNPS workflow (https://gnps.ucsd.edu/ProteoSAFe/static/gnps-splash.jsp). The parameters used are detailed at the following URL: http://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=c6438af750784d919dcd0ee0a783b4fc. Feature extraction parameters were optimized using MZmine2 (http://mzmine.sourceforge.net/) with a signal threshold of 2.0e5 and a 0.3-s minimum peak width. The mass tolerance was set to 10 ppm, and the maximum allowed retention time deviation was set to 10 s. For chromatographic deconvolution, the local minimum search algorithm was used with a minimum relative peak height of 1% and a minimum retention time range of 0.6 s. The maximum peak width was set to 1 min. After isotope peak removal, the peak lists of all samples were aligned with the above-mentioned retention time and mass tolerances. After the creation of a feature matrix containing the feature retention times and the exact mass and peak areas of the corresponding extracted ion chromatograms, the metadata of the samples (treatment type and duration) were added. The signal intensities of the features were normalized (probabilistic quotient normalization [PQN]) to an internal standard (m/z 278.189; real time, 3.81 min) for subsequent analysis.

18MS/MS annotations were as follows. Molecular features, in the form of MS/MS spectra, were putatively identified using MS2-based spectral library matches. The false-discovery rate (FDR) was estimated using a decoy database approach in GNPS and was found to be less than 1% above a cosine similarity score of 0.6 (GNPS job link, https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=feac48de4c9f45d485403e3feb7a470d). Therefore, we used a cosine score of 0.65 here. For level 1 annotation, as defined by the 2007 metabolomics standards initiative for differentially abundant metabolites, we purchased authentic standards of alpha-muricholic acid, chenodeoxycholic acid, cholic acid, lithocholic acid, deoxycholic acid, taurodeoxycholic acid, and the xenoestrogens enterodiol and enterolactone from Cayman Chemical (Ann Arbor, MI) and analyzed them using the same LC-MS/MS method described above. We then compared and verified the exact masses, fragmentation patterns, and retention times of those compounds to ensure correct annotations (Fig. S2 and Fig. S3).

19Statistical analysis was carried out as follows. QIIME 1.9.1 was used to perform principal-coordinate analysis (PCoA) (beta_diversity.py, a Gower dissimilarity metric) and a PERMANOVA test (compare_categories.py). The PCoA plots were visualized in EMPeror. Differential abundance analysis was performed using discrete FDRs.

20DNA extraction and 16S rRNA amplicon sequencing were done using Earth Microbiome Project (EMP) standard protocols (http://www.earthmicrobiome.org/protocols-and-standards/16s). In brief, DNA was extracted using the MO BIO PowerSoil DNA extraction kit (Carlsbad, CA). Amplicon PCR was performed on the V4 region of the 16S rRNA gene using the primer pair 515f to 806r with Golay error-correcting barcodes on the reverse primer. Amplicons were barcoded and pooled in equal concentrations for sequencing. The amplicon pool was purified with the MO BIO UltraClean PCR cleanup kit and sequenced on the Illumina HiSeq 2500 sequencing platform. Sequence data were demultiplexed and minimally quality filtered using the QIIME 1.9.1 script split_libraries_fastq.py, with a Phred quality threshold of 3 and default parameters to generate per-study FASTA sequence files.

21Feature detection and identification were performed as follows. The raw sequence data were processed using the Deblur workflow with default parameters in Qiita (https://qiita.ucsd.edu/). This generated a sub-operational taxonomic unit (sOTU) abundance per sample (BIOM format). Taxonomies for sOTUs were assigned using the sklearn-based taxonomy classifier (feature classifier plug-in) in QIIME 2. The sOTU table was rarefied to a depth of 2,000 sequences/sample to control for sequencing effort. A phylogeny was inferred using SATé-enabled phylogenetic placement, which was used to insert 16S Deblur sOTUs into Greengenes 13_8 at a 99% phylogeny.

22For statistical analysis, QIIME 2 was used to perform PCoA (unweighted UniFrac distances). QIIME 1.9.1 was used for the PERMANOVA test (compare_categories.py), Mantel test (compare_distance_matrices.py), and procrustes analysis (transform_coordinate_matrices.py). The PCoA and procrustes plots were visualized in EMPeror. Differential-abundance analysis was performed using discrete FDRs.

23The data generated in this study are available publicly under the following accession numbers: for metabolomics data, MSV000081482; for commercial standards, MSV000081853; and for microbiome data, ERP106495 (EBI database). Data analysis has been documented in Jupyter notebooks available on GitHub (https://github.com/knightlab-analyses/haddad_osa).