- 1Tulane Center of Biomedical Informatics and Genomics, Deming Department of Medicine, Tulane University School of Medicine, Tulane University, New Orleans, LA, United States
- 2Department of Endocrinology and Metabolism, The Third Affiliated Hospital of Southern Medical University, Guangzhou, China
- 3Center of Systems Biology, Data Information and Reproductive Health, School of Basic Medical Science, Central South University, Changsha, China
While the gut microbiome has been reported to play a role in bone metabolism, the individual species and underlying functional mechanisms have not yet been characterized. We conducted a systematic multi-omics analysis using paired metagenomic and untargeted serum metabolomic profiles from a large sample of 499 peri- and early post-menopausal women to identify the potential crosstalk between these biological factors which may be involved in the regulation of bone mineral density (BMD). Single omics association analyses identified 22 bacteria species and 17 serum metabolites for putative association with BMD. Among the identified bacteria, Bacteroidetes and Fusobacteria were negatively associated, while Firmicutes were positively associated. Several of the identified serum metabolites including 3-phenylpropanoic acid, mainly derived from dietary polyphenols, and glycolithocholic acid, a secondary bile acid, are metabolic byproducts of the microbiota. We further conducted a supervised integrative feature selection with respect to BMD and constructed the inter-omics partial correlation network. Although still requiring replication and validation in future studies, the findings from this exploratory analysis provide novel insights into the interrelationships between the gut microbiome and serum metabolome that may potentially play a role in skeletal remodeling processes.
Introduction
Osteoporosis is a progressive age-related condition associated with reduced bone mineral density (BMD) and increased susceptibility to low trauma fractures, which are the clinical endpoint of the disease. It represents the most prevalent metabolic bone disorder affecting >200 million people worldwide (Reginster and Burlet, 2006), and the burden is particularly large among postmenopausal women, which is mainly attributed to the reduced production of estrogen and other hormonal/metabolic changes that occur during menopause. It is estimated that at least one in three postmenopausal women have osteoporosis, and nearly half of those women will experience fragility fractures in their remaining lifetime (Melton et al., 1992). Dual-energy X-ray absorptiometry (DXA) derived BMD measurements of the hip and spine are the most frequently used metric for clinically diagnosing osteoporosis, as well as the most powerful known risk factor for predicting fracture risk (Kanis et al., 2005).
The gut microbiome, composed of the bacteria residing in the human gastrointestinal tract, is involved in a variety of diverse functions that are important for physiological wellbeing. There are several potential mechanisms through which the microbiome may impact bone metabolism, as previously reviewed (Hernandez et al., 2016; Chen et al., 2017). The microbiota can influence the intestinal absorption of essential minerals (e.g., calcium) that are important for maintaining skeletal homeostasis (Weaver, 2015), elicit immune responses which may alter the levels of inflammatory cytokines (e.g., TNF-α) that are important for bone health (Sjogren et al., 2012), produce metabolic byproducts (e.g., short chain fatty acids) which regulate critical cell signaling factors for bone remodeling processes (Lucas et al., 2018), and modulate the levels of hormones and neurotransmitters through the gut-brain axis (Cryan et al., 2019), including some (e.g., serotonin) that have been shown to interact with bone cells (Bliziotes, 2010). Although experimental animal models have provided compelling evidence that the gut microbiome may play a role in the regulation of bone mass (Sjogren et al., 2012), only a few limited studies have explored this relationship in humans (Wang et al., 2017; Das et al., 2019; Xu et al., 2020). While these early efforts reported significant differences in the microbial diversity between osteoporosis cases and healthy controls, they were generally limited by small sample sizes and the inability to reveal specific trait-associated bacteria or functional mechanisms.
Metabolomics enables the comprehensive profiling of the intermediate and end products of cellular metabolism. Since metabolites represent the downstream expression of genomic, transcriptomic, and proteomic factors, small changes in other omics may be amplified at the metabolomic level, enabling the detection of critical biomarkers or corresponding therapeutic target pathways closely related to disease risk (Johnson et al., 2016). However, the application of metabolomics for osteoporosis is rather limited. At present, most efforts are confined to animal experiments (Lv et al., 2016), although several early studies in humans have identified novel osteoporosis biomarkers involved in the metabolism of tryptophan, phenylalanine, lipids, and energy (Miyamoto et al., 2018; Moayyeri et al., 2018; Zhao Q et al., 2018; Gong et al., 2021). Notably, some studies demonstrated that the effects of the novel metabolites identified were more significant than classical bone turnover markers (Qi et al., 2016), supporting the crucial functions of small molecule metabolites in BMD regulation and osteoporosis prediction.
It is well established that changes in diet may be accompanied by shifts in the composition of the microbiome (Singh et al., 2017), but perhaps even more important is the resulting effect on the human metabolome. The diet contains many compounds that cannot be broken down by human digestive enzymes, and therefore pass to the gut where they are catabolized by the microbiota (Lamichhane et al., 2018). Some of the metabolic byproducts generated during these processes may then be absorbed into the circulating blood, where they can potentially impact human health. For instance, the gut metabolite Trimethylamine N-oxide (TMAO), regulated by dietary phosphatidylcholine intake, has been shown to promote the development of atherosclerosis (Wang et al., 2011; Tang et al., 2013). Based on these findings, a novel therapeutic approach was established to inhibit microbial production of TMAO (Wang et al., 2015). We hypothesized that there could be similar undiscovered mechanisms which contribute to the osteoporosis susceptibility.
Multi-omics integration analyses of microbiome and metabolite profiles collected from the same individuals are very much needed to elucidate the full range of interactions between these biological factors with respect to bone phenotypes. We integrated the paired gut microbiome and untargeted serum metabolite profiles from a large sample of peri- and early post-menopausal Chinese women to explore the crosstalk which may contribute to BMD variation at the femoral neck, the most common site for hip fracture, which is one of the most devastating types of osteoporotic fractures (LeBlanc et al., 2014). An overview of the study workflow is provided in Figure 1.
Figure 1 Overview of study workflow. 499 peri- and early post-menopausal women provided stool and blood samples for shotgun metagenomic sequencing and untargeted serum metabolomics profiling. Single omics association analyses were first conducted to identify microbes and metabolites that are associated with BMD. The paired microbiome and metabolite profiles were then integrated by performing a supervised feature selection with respect to BMD. The selected features were used to conduct inter-omics network analysis to explore the crosstalk between these biological factors.
Materials and Methods
Sample Recruitment
We randomly recruited 499 peri- and early post-menopausal Chinese women (aged 40 – 65) living in Guangzhou City, China. Perimenopausal refers to the menopause transition phase, characterized by irregular menstrual cycles, while postmenopausal is defined by the cessation of menstrual periods for >1 year (Lumsden, 2016). Women who had taken antibiotics or estrogens within three months of enrollment were excluded. We also excluded women with preexisting conditions relevant to bone mass development such as serious residual effects from cerebral vascular disease, diabetes mellitus, chronic renal failure, chronic liver failure, chronic lung disease, alcohol abuse, corticosteroid therapy for more than 6 months duration, evidence of other metabolic or inherited bone disease, rheumatoid arthritis, collagen disorders, and chronic gastrointestinal diseases. Each subject signed an informed consent, and the study protocol was approved by the Medical Ethics Committee of Southern Medical University.
BMD of the hip and spine were measured with DXA (Lunar, GE Healthcare, Madison, WI, USA) by trained and certified research staff. The machine was calibrated daily using a phantom scan for quality assurance, and the accuracy of BMD measurement was assessed by the coefficient of variation for repeated measurements, which was 0.89% for spine BMD. To minimize information loss from artificially dichotomizing individuals into low/high BMD groups, BMD was considered as a quantitative trait. BMD measurements were standardized to have a mean of zero and standard deviation of one, and the normalized values were used as the phenotype.
Each subject provided stool and blood samples for metagenomic and metabolomics analyses, respectively. Stool samples were frozen at -80°C after sample procurement until DNA extraction. To avoid variation due to circadian rhythm, which is known to affect the metabolome (Sahar and Sassone-Corsi, 2012), 10 ml of blood was drawn from each subject after >8 hours of overnight fasting. Serum was extracted from the blood samples according to the protein precipitation protocol (Bruce et al., 2009) developed for metabolomics analysis, aliquoted, and stored at -80°C until used for further analysis. The subjects also completed a questionnaire to collect relevant covariate information (e.g., demographic and lifestyle factors). Since sex hormones are involved in metabolism in general (Guarner-Lans et al., 2011), and bone metabolism more specifically (Drake et al., 2013), the serum levels of follicle stimulating hormone (FSH) and estradiol were measured using routine enzyme linked immunoassay ELISA kits (Immunodiagnostic Systems, Gaithersburg, MD, USA).
Metagenomic Sequencing
DNA was extracted from 200 mg of stool sample using the E.Z.N.A.® Stool DNA Kit (Omega, Norcross, GA, USA) following the manufacturer’s protocol. The total DNA was eluted in 50 μl of elution buffer (QIAGEN, Hilden, Germany) and stored at -80°C until metagenomic sequencing (LC-BIO Technologies Co. LTD., Hang Zhou, China). We constructed a fecal DNA library, and used HiSeq 4000 (Illumina, San Diego, CA, USA) with the paired end 150 bp strategy to conduct sequencing. Fecal DNA was fragmented using dsDNA Fragmentase (New England BioLabs, Ipswich, MA, USA) by incubating at 37°C for 30 min, and the DNA library was constructed by TruSeq Nano DNA LT Library Preparation Kit (Illumina, San Diego, CA, USA). Blunt-end DNA fragments were generated using a combination of fill-in reactions and exonuclease activity, and size selection was performed with the provided sample purification beads. An A-base was added to the blunt ends of the strands, preparing them for ligation to the indexed adapters. Each adapter contained a T-base overhang for ligating the adapter to the A-tailed fragmented DNA. The adapters were ligated to the fragments and the ligated products were amplified with PCR by the following conditions: initial denaturation at 95°C for 3 min, 8 cycles of denaturation at 98°C for 15 sec, annealing at 60°C for 15 sec, extension at 72°C for 30 sec, and then final extension at 72°C for 5 min.
The raw sequencing reads were then processed to obtain valid reads for further analysis by removing sequencing adapters with cutadapt v1.9 (Martin, 2011), trimming low quality reads using fqtrim v0.94 (Pertea, 2015), and aligning reads to the human reference genome (GRCh38/hg38) to remove host contamination with Bowtie2 v2.2.0 (Langmead and Salzberg, 2012). The quality filtered reads were de novo assembled to construct the metagenome for each sample using SPAdes v3.10.0 (Bankevich et al., 2012). All coding regions of metagenomic contigs were predicted using MetaGeneMark v3.26 (Zhu et al., 2010), and the coding sequences of all samples were clustered to obtain UniGenes with CD-HIT v4.6.1 (Fu et al., 2012). The UniGene abundances for a given sample were estimated by transcripts per million (TPM) based on the number of aligned reads, where k refers to the kth UniGene, r is the number of UniGene reads, and L is the UniGene length.
The DIAMOND+MEGAN approach was then applied for taxonomic annotation. The UniGenes were aligned against the NCBI non-redundant protein database with DIAMOND v0.9.20 (Buchfink et al., 2015). The quality of the alignments was determined based on the bit score, which represents the required size of a sequence database in which the current match could be found by chance, and E-value, which denotes the likelihood that a given sequence match is due purely to chance. The resulting alignments were then used as input for taxonomic binning using the lowest common ancestor (LCA) algorithm in MEGAN v6.12.3 (Huson et al., 2007), which places a read on the lowest taxonomic node in the NCBI taxonomy that lies above all taxa to which the read has a significant alignment. We note that while the limitations of this local sequence alignment approach have been documented (Koski and Golding, 2001), it is a standard protocol for taxonomic profiling (Bagci et al., 2021).
The microbiome data are relative abundances since the total number of read counts per sample is highly variable and constrained by the maximum number of reads the sequencer can provide (Gloor et al., 2017), and the data are considered compositional because the relative abundances of all bacteria species within each sample are proportions which have a unit sum. We eliminated the rare species with an average relative abundance <0.01% to reduce the extreme sparsity of the data and remove sequencing artifacts (Cao et al., 2020). The relative abundances were then normalized by the centered log-ratio (CLR) transformation, which has been shown to be effective in transforming the compositional data to be approximately multivariate normal (Gloor et al., 2017).
Serum Metabolomics Profiling
The serum samples were thawed on ice, and metabolites were extracted with 50% methanol buffer. 20 μl of sample was extracted with 120 μl of precooled 50% methanol, vortexed for 1 min, incubated at room temperature for 10 min, and stored overnight at –20°C. After centrifugation at 4,000g for 20 min, the supernatants were transferred into new 96-well plates and stored at –80°C prior to the liquid chromatography mass spectrometry (LC-MS) metabolomics analysis (LC-BIO Technologies Co. LTD., Hang Zhou, China). Pooled quality control samples were prepared by combining 10 μl of each extraction mixture. All chromatographic separations were performed using an ultra-performance liquid chromatography (UPLC) system (SCIEX, UK), and an ACQUITY UPLC BEH Amide column (100mm*2.1mm, 1.7μm, Waters, Wilmslow, UK) was used for the reversed phase separation. The column oven was maintained at 35°C, the flow rate was set to 0.4 ml/min, and the mobile phase consisted of solvent A (25mM ammonium acetate+25 mM NH4H2O) and solvent B (IPA: ACN=9:1+0.1% formic acid). Gradient elution conditions were set as follows: 0 ~ 0.5 min, 95% B; 0.5 ~ 9.5 min, 95% ~ 65% B; 9.5 ~ 10.5 min, 65% ~ 40% B; 10.5 ~ 12 min, 40% B; 12 ~ 12.2 min, 40% ~ 95% B; 12.2 ~ 15 min, 95% B.
A high-resolution tandem mass spectrometer Triple TOF5600plus (SCIEX, UK) was used to detect the metabolites eluted from the column in both positive and negative ion modes. The curtain gas was set to 30 PSI, ion source gas one and two were both set to 60 PSI, and the interface heater temperature was set to 650°C. The Ionspray voltage floating was 5000 V for positive ion mode and –4500 V for negative ion mode. The mass spectrometry data were acquired in IDA mode, and the TOF mass range was from 60 to 1200 Da. The survey scans were acquired in 150 ms, and as many as 12 product ion scans were collected if exceeding a threshold of 100 counts per second and with a 1+ charge-state. The total cycle time was fixed to 0.56 s. Four different time bins were summarized for each scan at a pulser frequency value of 11 kHz through monitoring of the 40 GHz multi-channel TDC detector with four-anode/channel detection. Dynamic exclusion was set for 4 s. During the acquisition, the mass accuracy was calibrated every 20 samples. To evaluate the stability of the LC-MS during the whole acquisition, a pooled quality control sample was acquired after every 10 samples.
The acquired MS data pretreatments including peak picking, peak grouping, retention time correction, and second peak grouping were performed using XCMS v3.16.1 (Smith et al., 2006). CAMERA v1.50 (Kuhl et al., 2012) was used to annotate the identified features with related isotopic peaks and adducts. Each ion was characterized by retention time and mass-to-charge ratios (m/z), and the intensities of each peak were recorded. Metabolite identification and data processing were performed using metaX v1.0.3 (Wen et al., 2017). The Human Metabolome Database (HMDB) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were used to annotate metabolites by performing a mass-based search with a weight tolerance of 10 ppm. To provide more confident and reproducible study findings, we retained metabolites with annotations that were validated using an in-house fragment spectrum library.
Metabolite features detected in <50% of quality control samples or <80% of biological samples were removed, and the remaining peaks with missing values were imputed with the k-nearest neighbor algorithm. Probabilistic quotient normalization was applied to minimize technical artifacts, and robust spline correction was used for the post-acquisition correction of batch effects. In addition, the relative standard deviations of the metabolite features were calculated across all quality control samples, and those >30% were removed. The remaining metabolite features were log transformed and scaled to have zero mean and unit variance, which is a common normalization technique (Zhao Q et al., 2018; Gong et al., 2021). The log transformation converts skewed data to symmetric, while scaling makes all metabolites of equal importance and enables comparison based on correlations (Li et al., 2016).
Microbiome Association Analysis
Individual microbes were tested for association with BMD using a constrained elastic net regression model, which is a commonly used feature selection approach with compositional covariates (Lin et al., 2014). The model imposes a sparsity penalty along with a constraint that the regression coefficients of the CLR-transformed relative abundances sum to zero, subject to . The elastic net regularization is a combination of both ridge and lasso penalty functions, where ridge results in a nonzero coefficient for every feature and lasso only assigns nonzero coefficients to the most strongly associated features. Since the penalized regression model does not provide conventional association p-values, partial Spearman correlation analysis was used to individually test each microbe selected in the initial feature screening.
Functional Profiling of Microbiota
The abundances of metabolic pathways in the microbiome community were profiled using the Human Microbiome Project Unified Metabolic Analysis Network (HUMAnN2) pipeline (Franzosa et al., 2018). HUMAnN2 first maps metagenomic reads to the pangenomes (Huang et al., 2014) of species identified by taxonomic profiling. The protein-coding sequences in these pangenomes have been pre-annotated to their respective UniRef90 families (Suzek et al., 2015), which serve as a non-redundant protein sequence database. Metagenomic reads that do not align to a known pangenome are subjected to a translated search against the full UniRef90 database. All hits are weighted by quality and sequence length to estimate the gene abundances. These genes are then annotated to metabolic enzymes and further analyzed to quantify the abundances of complete metabolic pathways obtained from MetaCyc (Caspi et al., 2016). HUMAnN2 assigns a coverage and abundance score for each pathway in each sample based on the detection of all its constituent genes. The coverage and abundance scores represent the number and abundance of complete copies of the pathway in each sample. Partial Spearman correlation analyses were used to test the associations between the pathway abundances and BMD.
Fecal Metabolite Imputation
The Model-based Genomically Informed High-dimensional Predictor of Microbial Community Metabolic Profiles (MelonnPan) approach (Mallick et al., 2019) was applied to predict the abundances of fecal metabolites from the microbiome gene abundances estimated by HUMAnN2. Elastic net prediction models were trained to select a sparse set of microbiome genes that are predictive for each fecal metabolite based on an independent set of 155 reference subjects for which both metagenomic and metabolomic profiling of the stool samples were both available (Franzosa et al., 2019). The fecal metabolite concentrations were then imputed as a linear combination of the microbiota gene abundances with weights learned from the training set. We retained the well predicted fecal metabolites, which had at least a moderate correlation (Spearman ρ >0.3) between the observed and imputed metabolite abundances in the training sample, as previously detailed (Mallick et al., 2019).
Metabolite Association Analysis
Partial least squares regression (PLS) is a multivariate approach which combines aspects of principal component analysis (PCA) and linear regression (Rohart et al., 2017). The principle is to extract a set of orthogonal components that have large covariance with the phenotype. PLS is well suited for the metabolomics analysis due to the high degree of correlation between functionally related metabolites (i.e., metabolites involved in the same metabolic pathways). A variable importance in projection (VIP) score is used to summarize the contribution of each feature to the model, which is computed as a weighted sum of the squared correlations between the PLS components and phenotype. Metabolites with VIP ≥2.0 were considered important for the phenotype. As a complementary approach, all metabolites were also individually tested using linear regression.
Coinertia Analysis
The global similarity between the gut microbiome and serum metabolome was investigated using coinertia analysis, which identifies successive axes of covariance between two datasets measured on a single group of subjects (Dray and Dufour, 2007). Principal coordinate analysis (PCoA) with Bray Curtis distance and PCA were applied to the microbiome and metabolite profiles, respectively, and the ordinations were used as input for the coinertia analysis. The coinertia analysis produces an RV coefficient, which is a multivariate extension of the squared Pearson correlation coefficient computed as where 0 < RV < 1. The coinertia between two hyperspaces is defined as the sum of the squared covariances between all variable pairs. Statistical significance of the RV coefficient was determined through a Monte-Carlo permutation test.
Supervised Multi-Omics Feature Selection
Canonical correlation analysis (CCA) has previously been proposed as a promising approach for performing integration analysis (Parkhomenko et al., 2009). Assuming two different data modalities measured on the same subjects, CCA seeks weighted linear combinations of the features from each dataset that have large correlation. However, the conventional CCA model assigns nonzero weights to every feature, which can result in overfitting for high dimensional data, and CCA is traditionally unsupervised since it does not take the phenotype information into consideration.
The overfitting issue can be addressed by introducing a sparsity penalty into the CCA model, which allows for the incorporation of feature selection. The sparse CCA model can then be further extended to be supervised (sCCA), such that the selected features are correlated across omics modalities with importance for a quantitative phenotype (Parkhomenko et al., 2009). The sCCA model is expressed as, uTXTYv subject to . X and Y denote the paired multi-omics datasets, u and v are the canonical vectors containing the weights for each feature, and Xu and Yv, taken to be the weighted linear combinations of features within each subject, are the canonical scores. The P1 and P2 represent lasso penalty functions on the canonical variates, and the resulting u and v are sparse for c1 and c2 sufficiently small. Q1 and Q2 denote subsets of features in X and Y that have large univariate correlation with the phenotype, and features that are not strongly associated with the phenotype are automatically assigned zero weights. The optimal tuning parameters for the model were selected by 10-fold cross validation.
Inter-Omics Network Analysis
The sCCA selected microbes and metabolites were used as input to construct the inter-omics Gaussian graphical model (GGM), where the edges represent partial correlations between features. The optimal GGM was selected by minimizing the extended Bayesian information criterion (EBIC) of unregularized GGM models (Foygel and Drton, 2010). We first selected the top 100 models by estimating a sparse inverse covariance matrix along a path of regularization parameters using the graph lasso penalty to select the significant edges. Each of these models was refit without regularization, and the model with the smallest EBIC was chosen as the optimal network.
Results
Sample Characteristics
The sample consisted of 499 peri- and early post-menopausal Chinese women that provided both stool and blood samples for metagenomic and metabolomic profiling (Table 1). 84% of these women were classified as postmenopausal (i.e., >1 year since final menstrual period), while 16% were still in the perimenopause transition period. The average time since menopause was 2.0 years (SD = 1.0), corresponding to the life stage when women typically begin to experience rapid bone loss (Melton et al., 1992). On average, the subjects were 53.0 years old (SD = 2.9) with body mass index (BMI) of 23.0 kg/m2 (SD = 2.9) and reported exercising approximately once per week (SD = 0.8). 62% of the women had undetectable estradiol levels (<18.35 pmol/L), an indicator of menopause, and the average level of FSH was 76.2 mIU/ml (SD = 32.2). BMI, exercise, time since menopause, estradiol, and FSH had significant bivariate associations with BMD (p-values <0.05).
Microbiome Association Analyses
After shotgun metagenomic sequencing of the stool DNA samples, we obtained approximately 7.35 giga base pairs of sequencing data per subject. Among >10,000 microbial features, there were 672 species with an average relative abundance >0.01%, which accounted for approximately 96% of the total microbiome across all subjects. 59.2% of these taxa belong to the Firmicutes phylum, 31.5% to Bacteroidetes, 6.1% to Proteobacteria, 2.2% to Actinobacteria, 0.5% to Fusobacteria, and 0.5% to Verrucomicrobia. On average, the Bacteroidetes and Firmicutes phyla accounted for 50% and 45% of the microbiome composition, respectively.
The SparCC approach (Friedman and Alm, 2012) was applied to explore the strength of relationships between the microbiota. SparCC accounts for the compositional nature of the data by approximating the correlations between the log-ratio transformed relative abundances. We observed that 8,697 pairs of taxa had strong positive correlations (ρ’s > 0.5 and p-values < 0.001), while 898 had strong negative correlations (ρ’s < –0.5 and p-values < 0.001). These relationships were visualized by the microbiome co-occurrence network (Figure 2). The underlying causes of these microbial interactions are complex. While mutualistic and phylogenetically related bacteria may sometimes co-occur, this is not always the case. Similarly, microbes with antagonistic relationships, such as those competing for the same niche, may sometimes have inverse associations, while in other circumstances they may actually co-occur due to variation in their shared environment (Levy and Borenstein, 2013). Co-exclusions can also arise due to incompatible abiotic factors in the microbiome community (Weiss et al., 2016).
Figure 2 Co-occurrence network showing strong interactions in the microbiome community (SparCC correlations >0.5 and <-0.5). Green/red edges correspond to positive/negative relationships. The color of each node corresponds to the phylum of the given microbe.
There were 44 taxa identified for potential association with BMD in the initial feature selection by the constrained elastic net regression model, including 22 which also had p-values <0.05 when tested individually using partial Spearman correlation analysis adjusted for relevant covariates (Table 2). Among the putative BMD associated microbes, 9 had FDR <0.05, and the remaining 13 had FDR <0.1. Several of the identified species including Bacteroides vulgatus, Bacteroides uniformis, Bacteroides fragilis, and Bacteroides massiliensis, all of which were negatively associated with BMD, were among the most abundant species in the microbiome with average relative abundances >1%. On the other hand, Firmicutes microbes, such as Clostridium leptum and Ruminococcus lactaris, were observed to be positively associated with BMD.
Functional profiling of the microbiome yielded pathway abundances for 516 metabolic pathways in the microbial community, all of which involve the microbiota producing metabolic byproducts from the catabolism of dietary components. Partial Spearman correlation analyses identified 22 pathways for putative association with BMD at a threshold of p-value <0.05 (Table 3). However, due to the number of pathways tested and the modest effect sizes, none of the pathway associations remained significant after multiple testing correction (FDR >0.2).
Predictive metabolomic profiling was performed to impute the fecal metabolite profiles based on the gene abundances in the microbiome communities. Among 80 predicted intestinal metabolites, 3 were identified for potential association with BMD based on VIP ≥2.0 in PLS, and 17 had p-values <0.05 with FDR of 0.2 when individually tested using linear regression (Table 4). Several of these compounds including butyrate, propionate, and valeric acid are short chain fatty acids (SCFAs), a special class of microbial byproducts that play an important role in gut and metabolic health (Blaak et al., 2020).
Serum Metabolite Association Analyses
Based on LC-MS untargeted metabolomics profiling of the serum samples, 3,202 unique metabolite features were identified in positive ion mode and 2,674 were detected in negative ion mode. Among the unique metabolite features, 381 had putatively confirmed identities. There were 12 serum metabolites identified for potential association with BMD based on VIP ≥2.0 in PLS, and 13 which had p-values <0.05 when tested individually by linear regression (Table 5). 8 serum metabolites were detected by both approaches, but none of the identified metabolites remained significant after the multiple testing correction (FDR >0.2). Notably, several putative BMD associated serum metabolites including 3-phenylpropanoic acid, which is primarily derived from the degradation of plant polyphenols (Trost et al., 2018), and glycolithocholic acid (Taylor and Green, 2018), a secondary bile acid, are intricately linked with the microbiota. While both these compounds were imputed in the fecal metabolite analysis, no significant associations were observed for the predicted fecal abundances.
Multi-Omics Integration
The coinertia analysis indicated that there were at least some statistically significant correlations between the microbiome and serum metabolome that should be further explored (RV = 0.064, p-value = 0.01). The sCCA model was then applied to the paired metagenomic and serum metabolite profiles, as has been done previously (Lee-Sarwar et al., 2019). We identified 14 bacteria species (Bacteroides vulgatus, Bacteroides fragilis, Bacteroides ovatus, Bacteroides xylanisolvens, Parabacteroides distasonis, Bacteroides sp. 4-3-47FAA, Bacteroides sp. 1-1-6, Bacteroides sp. 3-1-40A, Bacteroides sp. 2-2-4, Bacteroides sp. 3-1-23, Bacteroides sp. 1-1-30, Bacteroides sp. D22, Bacteroides sp. 2-1-22, and Fusobacterium ulcerans) and 6 serum metabolites (3-phenylpropanoic acid, hippuric acid, alpha-D-glucose, chenodeoxycholate, deoxycholic acid, and glycolithocholic acid) that were correlated across omics modalities with importance for BMD (Figure 3).
Figure 3 Correlation heatmap of sCCA selected microbes (rows) and metabolites (columns) that are correlated with each other and with BMD. Canonical loadings are provided for each feature, which represent the contributions to the inter-omics relationship. Positive correlations are represented by blue, negative correlations are shown in red, and intensity of the color represents the strength of association.
The metabolite and microbiome canonical scores, taken to be the weighted linear combination of features within each subject, were correlated with each other (ρ = 0.45, p-value < 0.001) and with BMD (βadj = 0.09 and 0.03 with p-values = 0.002 and 0.03, respectively), demonstrating the effectiveness of the supervised integrative feature selection. Based on the magnitude of the canonical loadings, which represent the contributions of each feature to the inter-omics relationship, the most important metabolites in the inter-omics relationship with respect to BMD were 3-phenylpropanoic acid and glycolithocholic acid, while the bacteria with the largest loadings were Fusobacterium ulcerans and Bacteroides fragilis, each of which were also detected in the single omics association analyses. The relationships between the abundances of these features and BMD after adjustment for covariates were visualized by added variable partial regression plots (Figure 4).
Figure 4 Added variable partial regression plots for (A) Bacteroides fragilis, (B) Fusobacterium ulcerans, (C) 3-Phenylpropanoic acid, and (D) Glycolithocholic acid. The plots illustrate the relationships between the abundance of a given microbe/metabolite and BMD after adjustment for age, BMI, exercise, years since menopause, FSH, and estradiol. The x-axes correspond to the residuals from regressing the microbe/metabolite on the covariates. The y-axes correspond to the residuals from regressing the phenotype on the covariates. The blue line represents the line of best fit from linear regression, and the corresponding confidence interval is shown in gray.
The GGM (Figure 5) had an edge density of 0.22, which represents the ratio of the number of edges and the number of possible edges, and a transitivity of 0.46, which is defined as the probability that adjacent nodes of a given node are connected. Fusobacterium ulcerans was positively connected to deoxycholic acid and negatively connected to both 3-phenylpropanoic acid and glycolithocholic acid. Bacteroides fragilis was negatively connected to glycolithocholic acid, while Bacteroides ovatus was negatively connected to 3-phenylpropanoic acid. Alpha-D-glucose and deoxycholic acid were positively and negatively connected to BMD, respectively.
Figure 5 Inter-omics Gaussian graphical model for sCCA selected features. The edges represent partial correlations, and significant edges were selected by the graph lasso penalty. Blue/red nodes correspond to microbe/metabolite features, BMD is shown in purple, and green/red edges correspond to positive/negative partial correlations. The edge width indicates the strength of association.
Discussion
In this systematic multi-omics analysis of a relatively large sample of peri- and early postmenopausal Chinese women, we characterized the microbiota, serum metabolites, and possible crosstalk between these biological factors that may influence bone physiology. To our knowledge, this is one of the first reports to integrate paired metagenomic and metabolomic profiles to provide novel insights into the molecular mechanisms of skeletal remodeling. The findings, although biologically plausible, still require replication and functional validation in future studies.
Many of the putative BMD associated microbes belong to the Bacteroides genus, which were inversely associated with bone mass. Bacteroides vulgatus and Bacteroides fragilis, which were identified in both the single omics and integrative analyses, have previously been reported to induce activation of the pro-inflammatory NF-κβ signaling pathway, which is associated with bone loss (Kim et al., 2002). Additionally, Bacteroides vulgatus was recently shown to increase serum levels of the bone resorption marker CTX-1 and decrease serum levels of the bone formation marker osteocalcin in vivo in an ovariectomized (OVX) mouse model (Lin et al., 2020). On the other hand, we observed a positive effect for the Firmicutes microbe Clostridium leptum, a probiotic species that is known to be an important producer of beneficial metabolic byproducts such as butyrate (Canani et al., 2012). We further observed a negative effect of Fusobacterium ulcerans, which also played an important role in the multi-omics integration analysis. Fusobacteria have been shown to promote M1 macrophage production via AKT2 signaling (Liu et al., 2019), which induces inflammation and has been associated with the development of osteoporosis (Yang and Yang, 2019).
To investigate the potential mechanisms by which the microbiota may influence BMD, we profiled the abundances of metabolic pathways in the microbial community and assessed their associations with BMD variation. Although the pathway associations were not significant at a stringent threshold accounting for multiple testing, several were still interesting due to their known roles in bone metabolism. We observed a positive association between BMD and several glycolytic pathways, such as glycogen biosynthesis/degradation, which are essential for cellular energy (Adeva-Andany et al., 2016). We further identified a negative association for urate biosynthesis, and it has been reported that uric acid induces intracellular oxidative stress and inflammatory cytokines that stimulate bone resorption by osteoclasts and inhibit bone formation by osteoblasts (Austermann et al., 2019). Lastly, we detected a positive association for L-methionine biosynthesis, an amino acid which has been shown to down-regulate NF- κβ signaling in osteoclast precursors to reduce bone loss (Vijayan et al., 2014).
Since microbial metabolites contribute to host-microbiome interactions (Rooks and Garrett, 2016), we performed metabolomic imputation to predict the intestinal metabolite profiles based on the observed genes in the microbiome community. Similar techniques are frequently used in genetic association studies to impute the gene expression levels based on genotype information (Gamazon et al., 2015). We identified positive associations between BMD and several SCFAs, which are exclusively produced by the microbiota through the breakdown of non-digestible dietary fiber. The SCFAs are potent signaling molecules that modulate host gene expression by interacting with various epigenetic factors such as DNA methylation and histone acetylation (Alenghat and Artis, 2014). Butyrate and propionate have previously been reported to induce metabolic alterations of osteoclasts that lead to down-regulation of crucial genes such as TRAF6 and NFATc1 (Lucas et al., 2018). Butyrate has also been shown to stimulate bone formation through regulatory T cell mediated regulation of Wnt10b expression (Tyagi et al., 2018). Furthermore, valeric acid was recently demonstrated to promote/inhibit osteoblast/osteoclast differentiation in vitro (Lin et al., 2020). We note that the precision of fecal metabolite imputation by MelonnPan has room for improvement (Yin et al., 2020). Additionally, the prediction models were trained in a different study population (mixed sexes, North American, some with irritable bowel disease), and therefore the accuracy of the imputation in this sample is unknown.
While the fecal metabolites are the most representative of the direct metabolic output of the gut, many of those compounds are excreted from the body without ever having any influence on human health. The serum metabolome, which includes both host and microbiota derived metabolites, provides a window into which gut metabolic byproducts are absorbed into the circulating blood to potentially impact host physiology. We observed that metabolites involved in energy metabolism, such as alpha-D-glucose, were positively associated with BMD. Energy metabolism is critical for bone remodeling, and it has also been demonstrated that there is a feedback loop where bone can act as an endocrine gland by secreting bone specific proteins such as osteocalcin and osteoprotegerin, which can regulate insulin function and glucose metabolism (Faienza et al., 2015). Additionally, we detected the amino acid histidinyl glycine, a conjugation of histidine and glycine. Targeted deletion of histidine decarboxylase in mice, which converts histidine to histamine, was found to increase bone formation and protect against bone loss (Fitzpatrick et al., 2003). Glycine is reported to improve bone health by increasing the production of collagen, which is a major building block of bone (Jennings et al., 2016).
Several of the BMD associated serum metabolites are involved in lipid metabolism, and accumulating evidence has demonstrated that alterations in lipid levels are associated with changes in bone metabolism (Tian and Yu, 2015). Lipid metabolism is regulated PPARγ, which inhibits the differentiation of osteoblasts and promotes the formation of adipocytes (Wan, 2010). Dodecanoic acid was recently reported to have a causal effect on BMD, which was validated using bone cells cultured in vitro as well as in vivo in an OVX mouse model (Gong et al., 2021). Palmitic acid is a saturated fatty acid which has been shown to increase bone loss by promoting osteoclast survival (Oh et al., 2010). The increase of LysoPC (18:0), a lysophoshpatidylcholine, is indicative of oxidative stress, and LysoPCs have been detected at elevated levels in the serum of osteoporotic mice (Zhao H et al., 2018).
Most notably, the serum metabolite analysis identified several microbiota-linked compounds for association with BMD. 3-phenylpropanoic acid, also known as hydrocinnamic acid, is mainly produced by the microbial catabolism of dietary polyphenols, which are acquired from plant-based food sources such as leafy greens, tea/coffee, wheat, berries, fruits, and other vegetables (Trost et al., 2018). Dietary polyphenols have been reported to reduce the risk of various age-related diseases and have a lot of promise for protecting against bone loss due to their antioxidant properties (Sato et al., 2011). One study found that phenolic acids in the serum of rats fed a blueberry diet stimulated osteoblast differentiation, resulting in elevated bone mass (Chen et al., 2010). Additionally, in vitro experiments with phenolic acids using bone marrow stroma cells demonstrated stimulation of osteoblast differentiation and inhibition of adipogenesis (Chen and Anderson, 2014). Plant polyphenols are also an established source of Hippuric acid, which was recently observed to inhibit osteoclast formation in vitro (Zhao et al., 2020). On the other hand, Hippuric acid can also be derived from aromatic organic acids such as phenylalanine and tryptophan (Wikoff et al., 2009).
Glycolithocholic acid and deoxycholic acid are secondary bile acids that are formed when primary bile acids produced by the liver, such as cholic acid and chenodeoxycholate, enter into the intestine via the bile duct and are acted on by the microbiota (Taylor and Green, 2018). It has previously been reported that bile acids are essential for the intestinal absorption of lipids and lipid-soluble compounds such as vitamin D (Nehring et al., 2007), and abnormal bile acid turnover has been linked with osteoporosis in postmenopausal women (Hanly et al., 2013). Lithocholic acid and deoxycholic acid have been shown to enhance and reduce calcium absorption, respectively (Marchionatti et al., 2018). Additionally, a growing body of evidence suggests that bile acids may regulate skeletal remodeling processes through direct interactions with osteoblasts and osteoclasts (Cho et al., 2013).
The network analysis revealed several inter-omics relationships that could potentially play a role in the regulation of BMD. First, we observed a positive connection between Fusobacterium ulcerans and deoxycholic acid. Sulfate esterification of bile acids in the liver enhances their excretion (Alnouti, 2009), and Fusobacteria have been reported to be involved in desulfatation, which keeps them in circulation (Jia et al., 2018). Second, we observed a negative connection between Fusobacterium ulcerans and 3-phenylpropanoic acid. It has previously been reported that polyphenol compounds have antimicrobial properties that inhibit the growth and virulence of Fusobacteria (Ben Lagha et al., 2017). Third, we observed a negative relationship between Bacteroides fragilis and glycolithocholic acid. Bacteroides are reported to promote deconjugation of bile acids, and individuals with higher abundance of Bacteroides have been shown to have lower plasma levels of secondary bile acid metabolites (Gu et al., 2017). Lastly, we observed a negative connection between 3-phenylpropanoic acid and Bacteroides ovatus. Previous studies have demonstrated that dietary polyphenols can inhibit the growth of Bacteroides microbes (Cardona et al., 2013). However, further analyses are needed to determine the precise mechanisms of these relationships and how they may be relevant for bone physiology.
Despite the novelty of this study in the bone field, there are several limitations that should be taken into consideration. First, among thousands of unique metabolite features, we were only able to produce high confidence annotations for a relatively small number, and there could be important compounds in the serum, especially gut metabolites, that were ignored. Second, we have only considered linear relationships between the molecular features. In the future, this could be addressed through nonlinear extensions of the applied methods, which may capture complex statistical dependencies that conventional correlation approaches fail to detect (Szekely et al., 2007). Third, we can only speculate about the directionality of the inter-omics relationships observed in the network based on previously reported findings. Lastly, it is unclear if the findings can be generalized to other populations. Different ethnicities have vastly different diets, and metabolomic patterns are influenced by various factors including age, genetics, and menopause (Auro et al., 2014).
In summary, we conducted a comprehensive multi-omics integration analysis and provided novel insights into the interactions between the gut microbiome and serum metabolome that may be relevant for the regulation of BMD. We hope that these findings will stimulate future studies to further explore the relationships between the microbiome and host omics factors that may be involved in bone health.
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 below: https://www.ebi.ac.uk/metagenomics/, PRJEB50761.
Ethics Statement
The studies involving human participants were reviewed and approved by Medical Ethics Committee of Southern Medical University. The patients/participants provided their written informed consent to participate in this study.
Author Contributions
JG conducted the data analysis and prepared the manuscript. XL, K-JS, and RG contributed to the data analysis plan. HS contributed to the manuscript revision. JS and H-MX managed and directed the study components conducted within their respective institutions. H-WD conceived, designed, and directed the whole project. All authors contributed to the article and approved the submitted version.
Funding
This research was benefited by the partial support by grants from the National Institutes of Health (P20GM109036, R01AR069055, U19AG055373, R01AG061917), the Science and Technology Program of Guangzhou, China (201604020007), the National Natural Science Foundation of China (81770878), and the National Key R&D Program of China (2016YFC1201805 and 2017YFC1001100).
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.
References
Adeva-Andany, M. M., González-Lucán, M., Donapetry-García, C., Fernández-Fernández, C., Ameneiros-Rodríguez, E. (2016). Glycogen Metabolism in Humans. BBA Clin. 5, 85–100. doi: 10.1016/j.bbacli.2016.02.001
Alenghat, T., Artis, D. (2014). Epigenomic Regulation of Host-Microbiota Interactions. Trends Immunol. 35, 518–525. doi: 10.1016/j.it.2014.09.007
Alnouti, Y. (2009). Bile Acid Sulfation: A Pathway of Bile Acid Elimination and Detoxification. Toxicol. Sci. 108, 225–246. doi: 10.1093/toxsci/kfn268
Auro, K., Joensuu, A., Fischer, K., Kettunen, J., Salo, P., Mattsson, H., et al. (2014). A Metabolic View on Menopause and Ageing. Nat. Commun. 5, 4708. doi: 10.1038/ncomms5708
Austermann, K., Baecker, N., Stehle, P., Heer, M. (2019). Putative Effects of Nutritive Polyphenols on Bone Metabolism In Vivo-Evidence From Human Studies. Nutrients 11, 1–14. doi: 10.3390/nu11040871
Bagci, C., Patz, S., Huson, D. H. (2021). DIAMOND+MEGAN: Fast and Easy Taxonomic and Functional Analysis of Short and Long Microbiome Sequences. Curr. Protoc. 1, e59. doi: 10.1002/cpz1.59
Bankevich, A., Nurk, S., Antipov, D., Gurevich, A. A., Dvorkin, M., Kulikov, A. S., et al. (2012). SPAdes: A New Genome Assembly Algorithm and Its Applications to Single-Cell Sequencing. J. Comput. Biol. 19, 455–477. doi: 10.1089/cmb.2012.0021
Ben Lagha, A., Haas, B., Grenier, D. (2017). Tea Polyphenols Inhibit the Growth and Virulence Properties of Fusobacterium Nucleatum. Sci. Rep. 7, 44815. doi: 10.1038/srep44815
Blaak, E. E., Canfora, E. E., Theis, S., Frost, G., Groen, A. K., Mithieux, G., et al. (2020). Short Chain Fatty Acids in Human Gut and Metabolic Health. Benef Microbes 11, 411–455. doi: 10.3920/BM2020.0057
Bliziotes, M. (2010). Update in Serotonin and Bone. J. Clin. Endocrinol. Metab. 95, 4124–4132. doi: 10.1210/jc.2010-0861
Bruce, S. J., Tavazzi, I., Parisod, V., Rezzi, S., Kochhar, S., Guy, P. A. (2009). Investigation of Human Blood Plasma Sample Preparation for Performing Metabolomics Using Ultrahigh Performance Liquid Chromatography/Mass Spectrometry. Anal. Chem. 81, 3285–3296. doi: 10.1021/ac8024569
Buchfink, B., Xie, C., Huson, D. H. (2015). Fast and Sensitive Protein Alignment Using DIAMOND. Nat. Methods 12, 59–60. doi: 10.1038/nmeth.3176
Canani, R. B., Di Costanzo, M., Leone, L. (2012). The Epigenetic Effects of Butyrate: Potential Therapeutic Implications for Clinical Practice. Clin. Epigenet. 4, 1–7. doi: 10.1186/1868-7083-4-4
Cao, Q., Sun, X., Rajesh, K., Chalasani, N., Gelow, K., Katz, B., et al. (2020). Effects of Rare Microbiome Taxa Filtering on Statistical Analysis. Front. Microbiol. 11, 607325. doi: 10.3389/fmicb.2020.607325
Cardona, F., Andres-Lacueva, C., Tulipani, S., Tinahones, F. J., Queipo-Ortuno, M. I. (2013). Benefits of Polyphenols on Gut Microbiota and Implications in Human Health. J. Nutr. Biochem. 24, 1415–1422. doi: 10.1016/j.jnutbio.2013.05.001
Caspi, R., Billington, R., Ferrer, L., Foerster, H., Fulcher, C. A., Keseler, I. M., et al. (2016). The MetaCyc Database of Metabolic Pathways and Enzymes and the BioCyc Collection of Pathway/Genome Databases. Nucleic Acids Res. 44, D471–D480. doi: 10.1093/nar/gkv1164
Chen, X., Anderson, J. J. (2014). Diet and the Bone Marrow Niche for Stem Cell Recruitment. J. Bone Miner Res. 29, 1041–1042. doi: 10.1002/jbmr.2234
Chen, Y. C., Greenbaum, J., Shen, H., Deng, H. W. (2017). Association Between Gut Microbiota and Bone Health: Potential Mechanisms and Prospective. J. Clin. Endocrinol. Metab. 102, 3635–3646. doi: 10.1210/jc.2017-00513
Chen, J. R., Lazarenko, O. P., Wu, X., Kang, J., Blackburn, M. L., Shankar, K., et al. (2010). Dietary-Induced Serum Phenolic Acids Promote Bone Growth via P38 MAPK/beta-Catenin Canonical Wnt Signaling. J. Bone Miner Res. 25, 2399–2411. doi: 10.1002/jbmr.137
Cho, S. W., An, J. H., Park, H., Yang, J. Y., Choi, H. J., Kim, S. W., et al. (2013). Positive Regulation of Osteogenesis by Bile Acid Through FXR. J. Bone Miner Res. 28, 2109–2121. doi: 10.1002/jbmr.1961
Cryan, J. F., O'riordan, K. J., Cowan, C. S. M., Sandhu, K. V., Bastiaanssen, T. F. S., Boehme, M., et al. (2019). The Microbiota-Gut-Brain Axis. Physiol. Rev. 99, 1877–2013. doi: 10.1152/physrev.00018.2018
Das, M., Cronin, O., Keohane, D. M., Cormac, E. M., Nugent, H., Nugent, M., et al. (2019). Gut Microbiota Alterations Associated With Reduced Bone Mineral Density in Older Adults. Rheumatol. (Oxf.) 58, 2295–2304. doi: 10.1093/rheumatology/kez302
Drake, M. A., Khosla, S., American Society for Bone and Mineral Research (2013). “The Role of Sex Steroids in the Pathogenesis of Osteoporosis,” in Primer on the Metabolic Bone Diseases and Disorders of Mineral Metabolism, 8th ed, vol. xxvi . Ed. Rosen, C. J. (Ames, Iowa: Wiley-Blackwell), 1078 p.
Dray, S., Dufour, A.-B. (2007). The Ade4 Package: Implementing the Duality Diagram for Ecologists. J. Stat. Software 22, 1–20. doi: 10.18637/jss.v022.i04
Faienza, M. F., Luce, V., Ventura, A., Colaianni, G., Colucci, S., Cavallo, L., et al. (2015). Skeleton and Glucose Metabolism: A Bone-Pancreas Loop. Int. J. Endocrinol. 2015, 758148. doi: 10.1155/2015/758148
Fitzpatrick, L. A., Buzas, E., Gagne, T. J., Nagy, A., Horvath, C., Ferencz, V., et al. (2003). Targeted Deletion of Histidine Decarboxylase Gene in Mice Increases Bone Formation and Protects Against Ovariectomy-Induced Bone Loss. Proc. Natl. Acad. Sci. U. S. A. 100, 6027–6032. doi: 10.1073/pnas.0934373100
Foygel, R., Drton, M. (2010). Extended Bayesian Information Criteria for Gaussian Graphical Models. Adv. Neural Inf. Process. Syst. 23, 2020–2028.
Franzosa, E. A., Mciver, L. J., Rahnavard, G., Thompson, L. R., Schirmer, M., Weingart, G., et al. (2018). Species-Level Functional Profiling of Metagenomes and Metatranscriptomes. Nat. Methods 15, 962–968. doi: 10.1038/s41592-018-0176-y
Franzosa, E. A., Sirota-Madi, A., Avila-Pacheco, J., Fornelos, N., Haiser, H. J., Reinker, S., et al. (2019). Gut Microbiome Structure and Metabolic Activity in Inflammatory Bowel Disease. Nat. Microbiol. 4, 293–305. doi: 10.1038/s41564-018-0306-4
Friedman, J., Alm, E. J. (2012). Inferring Correlation Networks From Genomic Survey Data. PloS Comput. Biol. 8, e1002687. doi: 10.1371/journal.pcbi.1002687
Fu, L., Niu, B., Zhu, Z., Wu, S., Li, W. (2012). CD-HIT: Accelerated for Clustering the Next-Generation Sequencing Data. Bioinformatics 28, 3150–3152. doi: 10.1093/bioinformatics/bts565
Gamazon, E. R., Wheeler, H. E., Shah, K. P., Mozaffari, S. V., Aquino-Michaels, K., Carroll, R. J., et al. (2015). A Gene-Based Association Method for Mapping Traits Using Reference Transcriptome Data. Nat. Genet. 47, 1091–1098. doi: 10.1038/ng.3367
Gloor, G. B., Macklaim, J. M., Pawlowsky-Glahn, V., Egozcue, J. J. (2017). Microbiome Datasets Are Compositional: And This Is Not Optional. Front. Microbiol. 8, 2224. doi: 10.3389/fmicb.2017.02224
Gong, R., Xiao, H.-M., Zhang, Y.-H., Zhao, Q., Su, K.-J., Lin, X., et al. (2021). Identification and Functional Characterization of Metabolites for Bone Mass in Peri-/Post Menopausal Chinese Women. J. Clin. Endocrinol. Metab. 106, e3159–e3177. doi: 10.1210/clinem/dgab146
Guarner-Lans, V., Rubio-Ruiz, M. E., Perez-Torres, I., Banos De Maccarthy, G. (2011). Relation of Aging and Sex Hormones to Metabolic Syndrome and Cardiovascular Disease. Exp. Gerontol. 46, 517–523. doi: 10.1016/j.exger.2011.02.007
Gu, Y., Wang, X., Li, J., Zhang, Y., Zhong, H., Liu, R., et al. (2017). Analyses of Gut Microbiota and Plasma Bile Acids Enable Stratification of Patients for Antidiabetic Treatment. Nat. Commun. 8, 1785. doi: 10.1038/s41467-017-01682-2
Hanly, R., Ryan, N., Snelling, H., Walker-Bone, K., Dizdarevic, S., Peters, A. M. (2013). Association Between Bile Acid Turnover and Osteoporosis in Postmenopausal Women. Nucl. Med. Commun. 34, 597–600. doi: 10.1097/MNM.0b013e3283608993
Hernandez, C. J., Guss, J. D., Luna, M., Goldring, S. R. (2016). Links Between the Microbiome and Bone. J. Bone Miner Res. 31, 1638–1646. doi: 10.1002/jbmr.2887
Huang, K., Brady, A., Mahurkar, A., White, O., Gevers, D., Huttenhower, C., et al. (2014). MetaRef: A Pan-Genomic Database for Comparative and Community Microbial Genomics. Nucleic Acids Res. 42, D617–D624. doi: 10.1093/nar/gkt1078
Huson, D. H., Auch, A. F., Qi, J., Schuster, S. C. (2007). MEGAN Analysis of Metagenomic Data. Genome Res. 17, 377–386. doi: 10.1101/gr.5969107
Jennings, A., Macgregor, A., Spector, T., Cassidy, A. (2016). Amino Acid Intakes Are Associated With Bone Mineral Density and Prevalence of Low Bone Mass in Women: Evidence From Discordant Monozygotic Twins. J. Bone Miner Res. 31, 326–335. doi: 10.1002/jbmr.2703
Jia, W., Xie, G., Jia, W. (2018). Bile Acid-Microbiota Crosstalk in Gastrointestinal Inflammation and Carcinogenesis. Nat. Rev. Gastroenterol. Hepatol. 15, 111–128. doi: 10.1038/nrgastro.2017.119
Johnson, C. H., Ivanisevic, J., Siuzdak, G. (2016). Metabolomics: Beyond Biomarkers and Towards Mechanisms. Nat. Rev. Mol. Cell Biol. 17, 451–459. doi: 10.1038/nrm.2016.25
Kanis, J. A., Borgstrom, F., De Laet, C., Johansson, H., Johnell, O., Jonsson, B., et al. (2005). Assessment of Fracture Risk. Osteoporos Int. 16, 581–589. doi: 10.1007/s00198-004-1780-5
Kim, J. M., Cho, S. J., Oh, Y.-K., Jung, H.-Y., Kim, Y.-J., Kim, N. (2002). Nuclear Factor-Kappa B Activation Pathway in Intestinal Epithelial Cells Is a Major Regulator of Chemokine Gene Expression and Neutrophil Migration Induced by Bacteroides Fragilis Enterotoxin. Clin. Exp. Immunol. 130, 59–66. doi: 10.1046/j.1365-2249.2002.01921.x
Koski, L. B., Golding, B. G. (2001). The Closest BLAST Hit Is Often Not the Nearest Neighbor. J. Mol. Evol. 52, 540–542. doi: 10.1007/s002390010184
Kuhl, C., Tautenhahn, R., Bottcher, C., Larson, T. R., Neumann, S. (2012). CAMERA: An Integrated Strategy for Compound Spectra Extraction and Annotation of Liquid Chromatography/Mass Spectrometry Data Sets. Anal. Chem. 84, 283–289. doi: 10.1021/ac202450g
Lamichhane, S., Sen, P., Dickens, A. M., Oresic, M., Bertram, H. C. (2018). Gut Metabolome Meets Microbiome: A Methodological Perspective to Understand the Relationship Between Host and Microbe. Methods 149, 3–12. doi: 10.1016/j.ymeth.2018.04.029
Langmead, B., Salzberg, S. L. (2012). Fast Gapped-Read Alignment With Bowtie 2. Nat. Methods 9, 357–359. doi: 10.1038/nmeth.1923
Leblanc, K. E., Muncie, H. L., Leblanc, L. L. (2014). Hip Fracture: Diagnosis, Treatment, and Secondary Prevention. Am. Family Physician 89945–951.
Lee-Sarwar, K. A., Kelly, R. S., Lasky-Su, J., Zeiger, R. S., O'connor, G. T., Sandel, M. T., et al. (2019). Integrative Analysis of the Intestinal Metabolome of Childhood Asthma. J. Allergy Clin. Immunol. 144, 442–454. doi: 10.1016/j/jaci.2019.02.032
Levy, R., Borenstein, E. (2013). Metabolic Modeling of Species Interaction in the Human Microbiome Elucidates Community-Level Assembly Rules. Proc. Natl. Acad. Sci. U.S.A. 110, 12804–12809. doi: 10.1073/pnas.1300926110
Lin, W., Shi, P., Feng, R., Li, H. (2014). Variable Selection in Regression With Compositional Covariates. Biometrika 101, 785–797. doi: 10.1093/biomet/asu031
Lin, X., Xiao, H.-M., Liu, H.-M., Lv, W.-Q., Greenbaum, J., Yuan, S.-J., et al. (2020). Gut Microbiota Impacts Bone via B.vulgatus-Valeric Acid-Related Pathways. medRxiv.
Li, B., Tang, J., Yang, Q., Cui, X., Li, S., Chen, S., et al. (2016). Performance Evaluation and Online Realization of Data-Driven Normalization Methods Used in LC/MS Based Untargeted Metabolomics Analysis. Sci. Rep. 6, 38881. doi: 10.1038/srep38881
Liu, L., Liang, L., Liang, H., Wang, M., Lu, B., Xue, M., et al. (2019). Fusobacterium Nucleatum Aggravates the Progression of Colitis by Regulating M1 Macrophage Polarization via AKT2 Pathway. Front. Immunol. 10, 1324. doi: 10.3389/fimmu.2019.01324
Lucas, S., Omata, Y., Hofmann, J., Bottcher, M., Iljazovic, A., Sarter, K., et al. (2018). Short-Chain Fatty Acids Regulate Systemic Bone Mass and Protect From Pathological Bone Loss. Nat. Commun. 9, 55. doi: 10.1038/s41467-017-02490-4
Lumsden, M. A. (2016). The NICE Guideline - Menopause: Diagnosis and Management. Climacteric J. Int. Menopause Soc. doi: 10.1080/13697137.2016.1222483
Lv, H., Jiang, F., Guan, D., Lu, C., Guo, B., Chan, C., et al. (2016). Metabolomics and Its Application in the Development of Discovering Biomarkers for Osteoporosis Research. Int. J. Mol. Sci. 17, 1–22. doi: 10.3390/ijms17122018
Mallick, H., Franzosa, E. A., Mclver, L. J., Banerjee, S., Sirota-Madi, A., Kostic, A. D., et al. (2019). Predictive Metabolomic Profiling of Microbial Communities Using Amplicon or Metagenomic Sequences. Nat. Commun. 10, 3136. doi: 10.1038/s41467-019-10927-1
Marchionatti, A., Rivoira, M., Rodriguez, V., Perez, A., Tolosa De Talamoni, N. (2018). Molecular Mechanisms Triggered by Bile Acids on Intestinal Ca2+ Absorption. Curr. Medi Chem. 25, 2122–2132. doi: 10.2174/0929867324666171116125131
Martin, M. (2011). Cutadapt Removes Adapter Sequences From High-Throughput Sequencing Reads. EMBnet.journal 17, 10–12. doi: 10.14806/ej.17.1.200
Melton, L., Chrischilles, E. A., Cooper, C., Lane, A. W., Riggs, L. B. (1992). Perspective. How Many Women Have Osteoporosis? J. Bone Mineral Res. 7, 1005–1010. doi: 10.1002/jbmr.5650070902
Miyamoto, T., Hirayama, A., Sato, Y., Koboyashi, T., Katsuyama, E., Kanagawa, H., et al. (2018). Metabolomics-Based Profiles Predictive of Low Bone Mass in Menopausal Women. Bone Rep. 9, 11–18. doi: 10.1016/j.bonr.2018.06.004
Moayyeri, A., Cheung, C. L., Tan, K. C., Morris, J. A., Cerani, A., Mohney, R. P., et al. (2018). Metabolomic Pathways to Osteoporosis in Middle-Aged Women: A Genome-Metabolome-Wide Mendelian Randomization Study. J. Bone Miner Res. 33, 643–650. doi: 10.1002/jbmr.3358
Nehring, J. A., Zierold, C., Deluca, H. F. (2007). Lithocholic Acid Can Carry Out In Vivo Functions of Vitamin D. Proc. Natl. Acad. Sci. U.S.A. 104, 10006–10009. doi: 10.1073/pnas.0703512104
Oh, S.-R., Sul, O.-J., Kim, Y.-Y., Kim, H.-J., Yu, R., Suh, J.-H., et al. (2010). Saturated Fatty Acids Enhance Osteoclast Survival. J. Lipid Res. 51.
Parkhomenko, E., Tritchler, D., Beyene, J. (2009). Sparse Canonical Correlation Analysis With Application to Genomic Data Integration. Stat. Appl. Genet. Mol. Biol. 8, 1–34. doi: 10.2202/1544-6115.1406
Qi, H., Bao, J., An, G., Ouyang, G., Zhang, P., Wang, C., et al. (2016). Association Between the Metabolome and Bone Mineral Density in Pre- and Post-Menopausal Chinese Women Using GC-MS. Mol. Biosyst. 12, 2265–2275. doi: 10.1039/C6MB00181E
Reginster, J. Y., Burlet, N. (2006). Osteoporosis: A Still Increasing Prevalence. Bone 38, S4–S9. doi: 10.1016/j.bone.2005.11.024
Rohart, F., Gautier, B., Singh, A., Le Cao, K. A. (2017). Mixomics: An R Package for 'Omics Feature Selection and Multiple Data Integration. PloS Comput. Biol. 13, e1005752. doi: 10.1371/journal.pcbi.1005752
Rooks, M. G., Garrett, W. S. (2016). Gut Microbiota, Metabolites and Host Immunity. Nat. Rev. Immunol. 16, 341–352. doi: 10.1038/nri.2016.42
Sahar, S., Sassone-Corsi, P. (2012). Regulation of Metabolism: The Circadian Clock Dictates the Time. Trends Endocrinol. Metab. 23, 1–8. doi: 10.1016/j.tem.2011.10.005
Sato, Y., Itagaki, S., Kurokawa, T., Ogura, J., Kobayashi, M., Hirano, T., et al. (2011). In Vitro and In Vivo Antioxidant Properties of Chlorogenic Acid and Caffeic Acid. Int. J. Pharm. 403, 136–138. doi: 10.1016/j.ijpharm.2010.09.035
Singh, R. K., Chang, H. W., Yan, D., Lee, K. M., Ucmak, D., Wong, K., et al. (2017). Influence of Diet on the Gut Microbiome and Implications for Human Health. J. Transl. Med. 15, 73. doi: 10.1186/s12967-017-1175-y
Sjogren, K., Engdahl, C., Henning, P., Lerner, U. H., Tremaroli, V., Lagerquist, M. K., et al. (2012). The Gut Microbiota Regulates Bone Mass in Mice. J. Bone Miner Res. 27, 1357–1367. doi: 10.1002/jbmr.1588
Smith, C. A., Want, E. J., O’maille, G., Abagyan, R., Siuzdak, G. (2006). XCMS: Processing Mass Spectrometry Data for Metabolite Profiling Using Nonlinear Peak Alignment, Matching, and Identification. Anal Chem. 78, 779–787. doi: 10.1021/ac051437y
Suzek, B. E., Wang, Y., Huang, H., Mcgarvey, P. B., Wu, C. H., Uniprot, C. (2015). UniRef Clusters: A Comprehensive and Scalable Alternative for Improving Sequence Similarity Searches. Bioinformatics 31, 926–932. doi: 10.1093/bioinformatics/btu739
Szekely, G. J., Rizzo, M. L., Bakirov, N. K. (2007). Measuring and Testing Dependence by Correlation of Distances. Ann. Stat 35, 2769–2794. doi: 10.1214/009053607000000505
Tang, W. H., Wang, Z., Levison, B. S., Koeth, R. A., Britt, E. B., Fu, X., et al. (2013). Intestinal Microbial Metabolism of Phosphatidylcholine and Cardiovascular Risk. N Engl. J. Med. 368, 1575–1584. doi: 10.1056/NEJMoa1109400
Taylor, S. A., Green, R. M. (2018). Bile Acids, Microbiota, and Metabolism. Hepatology 68, 1229–1231. doi: 10.1002/hep.30078
Tian, L., Yu, X. (2015). Lipid Metabolism Disorders and Bone Dysfunction–Interrelated and Mutually Regulated (Review). Mol. Med. Rep. 12, 783–794. doi: 10.3892/mmr.2015.3472
Trost, K., Ulaszewska, M. M., Stanstrup, J., Albanese, D., De Filippo, C., Tuohy, K. M., et al. (2018). Host: Microbiome Co-Metabolic Processing of Dietary Polyphenols - An Acute, Single Blinded, Cross-Over Study With Different Doses of Apple Polyphenols in Healthy Subjects. Food Res. Int. 112, 108–128. doi: 10.1016/j.foodres.2018.06.016
Tyagi, A. M., Yu, M., Darby, T. M., Vaccaro, C., Li, J. Y., Owens, J. A., et al. (2018). The Microbial Metabolite Butyrate Stimulates Bone Formation via T Regulatory Cell-Mediated Regulation of WNT10B Expression. Immunity 49, 1116–1131 e1117. doi: 10.1016/j.immuni.2018.10.013
Vijayan, V., Khandelwal, M., Manglani, K., Gupta, S., Surolia, A. (2014). Methionine Down-Regulates TLR4/MyD88/NF-kappaB Signalling in Osteoclast Precursors to Reduce Bone Loss During Osteoporosis. Br. J. Pharmacol. 171, 107–121. doi: 10.1111/bph.12434
Wan, Y. (2010). PPARgamma in Bone Homeostasis. Trends Endocrinol. Metab. 21, 722–728. doi: 10.1016/j.tem.2010.08.006
Wang, Z., Klipfell, E., Bennett, B. J., Koeth, R., Levison, B. S., Dugar, B., et al. (2011). Gut Flora Metabolism of Phosphatidylcholine Promotes Cardiovascular Disease. Nature 472, 57–63. doi: 10.1038/nature09922
Wang, Z., Roberts, A. B., Buffa, J. A., Levison, B. S., Zhu, W., Org, E., et al. (2015). Non-Lethal Inhibition of Gut Microbial Trimethylamine Production for the Treatment of Atherosclerosis. Cell 163, 1585–1595. doi: 10.1016/j.cell.2015.11.055
Wang, J., Wang, Y., Gao, W., Wang, B., Zhao, H., Zeng, Y., et al. (2017). Diversity Analysis of Gut Microbiota in Osteoporosis and Osteopenia Patients. PeerJ 5, e3450. doi: 10.7717/peerj.3450
Weaver, C. M. (2015). Diet, Gut Microbiome, and Bone Health. Curr. Osteoporosis Rep. 13, 125–130. doi: 10.1007/s11914-015-0257-0
Weiss, S., Van Treuren, W., Lozupone, C., Faust, K., Friedman, J., Deng, Y., et al. (2016). Correlation Detection Strategies in Microbial Data Sets Vary Widely in Sensitivity and Precision. ISME J. 10, 1669–1681. doi: 10.1038/ismej.2015.235
Wen, B., Mei, Z., Zeng, C., Liu, S. (2017). Metax: A Flexible and Comprehensive Software for Processing Metabolomics Data. BMC Bioinf. 18, 183. doi: 10.1186/s12859-017-1579-y
Wikoff, W. R., Anfora, A. T., Liu, J., Schultz, P. G., Lesley, S. A., Peters, E. C., et al. (2009). Metabolomics Analysis Reveals Large Effects of Gut Microflora on Mammalian Blood Metabolites. Proc. Natl. Acad. Sci. U.S.A. 106, 3698–3703. doi: 10.1073/pnas.0812874106
Xu, Z., Xie, Z., Sun, J., Huang, S., Chen, Y., Li, C., et al. (2020). Gut Microbiome Reveals Specific Dysbiosis in Primary Osteoporosis. Front. Cell Infect. Microbiol. 10, 160. doi: 10.3389/fcimb.2020.00160
Yang, D. H., Yang, M. Y. (2019). The Role of Macrophage in the Pathogenesis of Osteoporosis. Int. J. Mol. Sci. 20. doi: 10.3390/ijms20092093
Yin, X., Altman, T., Rutherford, E., West, K. A., Wu, Y., Choi, J., et al. (2020). A Comparative Evaluation of Tools to Predict Metabolite Profiles From Microbiome Sequencing Data. Front. Microbiol. 11, 595910. doi: 10.3389/fmicb.2020.595910
Zhao, H., Lazarenko, O. P., Chen, J. R. (2020). Hippuric Acid and 3-(3-Hydroxyphenyl) Propionic Acid Inhibit Murine Osteoclastogenesis Through RANKL-RANK Independent Pathway. J. Cell Physiol. 235, 599–610. doi: 10.1002/jcp.28998
Zhao, H., Li, X., Zhang, D., Chen, H., Chao, Y., Wu, K., et al. (2018). Integrative Bone Metabolomics-Lipidomics Strategy for Pathological Mechanism of Postmenopausal Osteoporosis Mouse Model. Sci. Rep. 8, 16456. doi: 10.1038/s41598-018-34574-6
Zhao, Q., Shen, H., Su, K. J., Zhang, J. G., Tian, Q., Zhao, L. J., et al. (2018). Metabolomic Profiles Associated With Bone Mineral Density in US Caucasian Women. Nutr. Metab. (Lond) 15, 57. doi: 10.1186/s12986-018-0296-5
Keywords: metagenomics, metabolomics, data integration, osteoporosis, bone
Citation: Greenbaum J, Lin X, Su K-J, Gong R, Shen H, Shen J, Xiao H-M and Deng H-W (2022) Integration of the Human Gut Microbiome and Serum Metabolome Reveals Novel Biological Factors Involved in the Regulation of Bone Mineral Density. Front. Cell. Infect. Microbiol. 12:853499. doi: 10.3389/fcimb.2022.853499
Received: 12 January 2022; Accepted: 21 February 2022;
Published: 16 March 2022.
Edited by:
Steven Gill, University of Rochester, United StatesReviewed by:
Nar Singh Chauhan, Maharshi Dayanand University, IndiaCecilia Noecker, University of California, San Francisco, United States
Copyright © 2022 Greenbaum, Lin, Su, Gong, Shen, Shen, Xiao and Deng. 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.
*Correspondence: Hong-Wen Deng, hdeng2@tulane.edu; Jie Shen, shenjiedr@163.com; Hong-Mei Xiao, hmxiao@csu.edu.cn
†These authors have contributed equally to this work