Alteration of gut microbiota in elderly postmenopausal women with osteoporosis


 Osteoporosis is a common bone disorder worldwide and causes bone fragility and fracture. Gut microbiota colonizes the gastrointestinal tract, and is associated with bone metabolism and osteoporosis. In our study, the alteration of gut microbiota in osteoporosis and its effects on bone metabolism were investigated. A total of 36 elderly postmenopausal osteoporotic women and 12 healthy controls were recruited, and their fecal samples were collected for 16S rRNA gene sequencing of gut microbiota on Illumina MiSeq platform. The venous blood and urine samples were also collected to determine the biochemical indexes. There was no obvious difference in Alpha diversity in the experiment group and control group, while differential Beta diversity was observed. The Partial Least Squares Discriminant Analysis (PLS-DA) model and variable importance in projection (VIP) scores showed that the osteoporotic women and healthy controls had different genera of Erysipelotrichaceae , Rothia , and Eubacterium . The metabolic function prediction of gut microbiota indicated that the experiment group had 634 unique functional categories, while the control group had 13 unique functional categories. The biochemical measurement revealed that the osteoclast activity indexes including urine N-terminal telopeptides of type I collagen (NTX), serum NTX, and serum C-terminal telopeptides of type I collagen (CTX) in the experiment group were higher than those in the control group, indicating that the osteoclast activity in osteoporotic women was increased. In addition, the correlation analysis of microbial metabolism with phenotypes showed the pathways in the significant modules of magenta, red, pink, and yellow were positively correlated with urine phosphate, urine creatinine, urine creatinine, serum calcium and other biochemical indexes. Collectively, our study identified the different genera between postmenopausal osteoporotic women and healthy controls, which might be potential targets for the treatment of osteoporosis.


Introduction
Osteoporosis is a common bone disorder worldwide with the characteristics of low bone strength, degeneration of bone tissues, decreased bone mass (Coughlan and Dockery, 2014). It is known that osteoporosis results in increased bone fragility and high risk of fracture, which affects more than 200 million people worldwide and remains to be a societal problem of health (Minisola et al., 2017;Kastner et al., 2010). As an aging-related disease, osteoporosis is more frequent in the elderly, especially among postmenopausal women (Wei et al., 2016;Wang et al., 2013;Gallagher and Levine, 2011). Approximately 50% of women aged over 50 years will have osteoporotic fracture (Sambrook and Cooper, 2006). The risk factors consist of age, gender, race, and other factors which may disturb the balance between bone loss and bone formation (e.g., lack of oestrogen, serum calcium level alteration, prolonged use of glucocorticoids). Currently, drugs are the main treatment for osteoporosis, however, the curative effect is limited, accompanied by side effects and high costs (Wu et al., 2013). Thus, huge efforts have been devoted to explore novel and effective targets for the treatment of osteoporosis.
The gut microbiota denotes the microbes colonizing the gastrointestinal tract and consists of anaerobes, aerobes and facultative anaerobes (Xu et al., 2017). Acquired at birth, gut microbiota coevolves with the host and is recognized as an organ which influences the physiology of host through multiple approaches (Ley et al., 2008). Some researchers have found that the gut microbiota is involved in bone mass and metabolism, indicating that gut microbiota may be associated with osteoporosis. Sjögren K et al showed that gut microbiota played a regulatory role in bone mass by acting on the immune status and resorption of bone in mice (Sjögren et al., 2012). Charles JF et al revealed that gut microbiota acted as anabolic stimulus to the skeleton in mice and enhanced the formation and growth of bone (Yan et al., 2016). It was found that variations of the gut microbiome in mice could weaken the bone strength and tissue material properties (Guss et al., 2017). Pacifici R et al demonstrated that gut microbiota was associated with the bone loss in mice with sex steroid deficiency, which could be improved by probiotics (Li et al., 2016). Although the studies mentioned above have explained the relationship between gut microbiota and osteoporosis, they were all performed based on animal experiments. The researches on human might unfold the association between gut microbiota and osteoporosis more directly and convincingly.
Our study investigated the influence of gut microbiota on the bone metabolism and osteoporosis among the elderly postmenopausal women through analysis of gut microbiota by next-generation sequencing technology (NGS) and biochemical detection. We aimed to elucidate the distribution and composition of gut microbiota which contributed to the improvement of bone metabolism and osteoporosis, and further implicate the possibility and feasibility of gut microbiota as promising targets for the treatment of osteoporosis based on current research.

Ethics statement
This study was approved by the Ethics Committee of Zhongda Hospital Southeast University and performed according to the Declaration of Helsinki. Written informed consents were obtained from all participants.

Study subjects
A total of 36 elderly postmenopausal osteoporotic women (mean age: 74.42 ± 7.48) and 12 women without osteoporosis (mean age: 64.0 ± 4.86) were recruited as the experiment and control groups respectively. Bone mineral density (BMD, T-Score ≤ -2.5 SD) was used for diagnosis of osteoporosis according to the guidelines for the diagnosis and treatment of primary osteoporosis (Baum and Peters, 2008). The detailed clinical data for individual human subject was shown in Table S1.

Biochemical measurement
The venous blood and urine samples of patients in the experiment and control groups were collected.
Then the indexes of bone metabolism including serum and urine calcium, serum and urine creatinine, serum and urine phosphate, serum 25-Hydroxy Vitamin D, serum and urine N-and C-terminal telopeptides of type I collagen (CTX and NTX), serum osteocalcin, calcitonin, parathyroid hormone (PTH), urine microalbumin and urine creatinine were determined. The indexes above were analyzed by R software.

DNA extraction
The fecal samples were collected and stored at -80. Total DNA of the fecal samples was extracted using InviMag stool DNA kit (Invitek, Germany). The quality and quantity of DNA samples were detected by 0.8% agarose gel electrophoresis and NanoDrop spectrophotometer (Thermo, USA) respectively.
Sequencing of the V3 and V4 regions of 16S rRNA gene The V3 and V4 regions of 16S rRNA gene were amplified using the Q5 high-fidelity DNA polymerase (New England Biolabs), and the number of amplification cycles was controlled to ensure the low amplification cycles and consistent amplification conditions. The PCR products were detected by 2% agarose gel electrophoresis and purified by gel extraction kit (AXYGEN, USA). The real-time quantitative polymerase chain reaction (RT-qPCR) was performed according to the instructions of Quant-iT PicoGreen dsDNA Assay Kit (Invitrogen, USA) on Microplate reader (BioTek, FLx800). Then the samples were mixed proportionally based on the sequencing requirements and RT-qPCR results.
After preparation of the sequencing library of the V3 and V4 regions of 16S rRNA gene, the samples were sequenced using the Illumina MiSeq System (Illumina Inc., USA) with paired-end 300-bp reads (Ong et al., 2018).

Bioinformatics analysis of sequencing data
Quantitative Insights Into Microbial Ecology software (QIIME, v1.8.0, http://qiime.org/) was used to filter the sequencing data with the following criteria: the reads were longer than 150 bp; there was no fuzzy N base; reads with more than 8 continuous identical bases or 1 primer 5' base mismatch were excluded. Chimera was removed by USEARCH (v5.2.236, http://www.drive5.com/usearch/), and the length distribution of reads was analyzed by R language.
We determined the operational taxonomic units (OTUs) using UCLUST with a similarity threshold of 97%, and selected the most abundant sequence of each OTU as the representative one. Alpha diversity was calculated by QIIME, and rarefaction, Chao1 index, and ACE index were calculated. The species accumulation curves were plotted using R software. Principal component analysis (PCA) was conducted using R software to analyze the Beta diversity. Partial Least Squares Discriminant Analysis (PLS-DA) model was constructed using R software and variable importance in projection (VIP) scores were calculated to screen the different genera between experiment and control groups. The metabolic function of microbiota was predicted by Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) and the species abundance of the original data was adjusted.
Weighted gene co-expression network analysis (WGCNA) was performed to analyze the correlation of microbial metabolism with phenotypes by using WGCNA R package. The abundance matrix of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways was used as input matrix, and the biochemical indexes were used as clinical information.

Statistical analysis
All data were expressed as mean ± standard deviation. Statistical analyses were performed using GraphPad Prism 6.0 statistical software (GraphPad Software Inc, La Jolla, CA, USA). Comparisons among multiple groups were accessed by two-way analysis of variance followed by Bonferroni test, while comparisons between two groups were evaluated by Student t test. p value < 0.05 was considered to indicate a statistically significant difference.

Sample Phylum
The rarefaction curves of the samples were plotted. As shown in Figure S1, the curves were flat, indicating that the OTU numbers did not increase significantly with the increase of sequencing depth, and the current sequencing depth was able to reflect the OTU number. The species accumulation curves of the samples were shown in Figure S2A. The result revealed that the species number did not significantly increase with the increase of sample size. Thus, the current sample size was able to reflect the abundance of species. The rank abundance curves were plotted based on the OTU abundance of samples and shown in Figure S2B.
After calculation of Chao1 and ACE indexes, it was found that there was no significant difference in Chao1 and ACE indexes of the experiment and control groups (Table 3, Fig. 2, p > 0.05). Difference in Beta diversity was observed between osteoporotic women and healthy controls PCA was performed for Beta analysis by transforming OTU abundance matrix into two and three dimensional images. As shown in Figure S3, the long distance between the experiment group and control group suggested that the two groups had obviously distinct gut microbiota.

Differential microbiota between osteoporotic women and healthy controls
The PLS-DA model was constructed based on the data of species abundance and grouping. The top 20 different genera were selected after calculation of the VIP scores (Table 4). Among them, Erysipelotrichaceae, Rothia, and Eubacterium showed higher VIP scores, indicating that there were obvious differences in these genera between osteoporotic women and healthy controls. As shown in Figure S4, the long distance between blue dots and orange dots, which denoted the control group and experiment group respectively, suggested that this model could classify the samples accurately. Table 4 The top 20 different genera between osteoporotic women and healthy controls.

Metabolic function prediction of gut microbiota
PICRUSt was used to predict the KEGG functional categories of the microbiota. The relative abundance of KEGG functional categories in experiment and control groups were shown in Fig. 3A.
Clustering analysis was performed on the top 50 functional categories with high abundance and the heatmap was shown in Fig. 3B. Then the common and unique function number of each group was calculated. As shown in Fig. 3C, the experiment group had 634 unique functional categories, while the control group had 13 unique functional categories. There were 5,456 common functional categories between the two groups.

Osteoporotic women presented increased osteoclast activity
The biochemical indexes of serum and urine of experiment and control groups were measured. As shown in Fig. 4, the osteoclast activity indexes including urine NTX (p < 0.05), serum NTX, and serum CTX in the experiment group were higher than those in the control group. The results indicated that compared with healthy controls, the osteoclast activity in osteoporotic women was increased. In addition, the serum calcium in the control group was significantly higher than that in the experiment group (P < 0.05).

Correlation of microbial metabolism with biochemical indexes
The input matrix included 48 samples and 6,909 KEGG pathways. Clustering analysis of the samples showed that there were no outlier samples (Fig. 5). Soft-threshold was selected as β = 16 to ensure it met the criterion of scale-free network, with the correlation coefficient between log(k) and log(p(k)) greater than 0.85 (Fig. 6A). k represented the intramodular connectivity and p(k) represented the frequency distribution of the connectivity. After transformation from expression matrix to adjacency matrix, the adjacency matrix was transformed to topological matrix. Based on topological overlap measure (TOM), hierarchical clustering analysis was performed. Dynamic branch cutting method was adopted to identity different modules, and the module eigengene (ME) was calculated as the representative gene of each module. Clustering analysis on the modules was carried out, with the highly correlated modules merged (height = 0.25). A total of 10 modules were obtained, and the grey module represented the pathways that could not be clustered into other modules (Fig. 6B).
Based on the ME, the module-phenotype relationship was calculated. The information of phenotype was shown in Table 5. As shown in Fig. 6C, the magenta, red, pink, and yellow modules were positively correlated with urine phosphate, urine creatinine, urine creatinine, serum calcium (correlation coefficients: 0.31, 0.33, 0.36, 0.39; p values: 0.03, 0.02, 0.01, 0.007, respectively) and other biochemical indexes. The pathways in these significant modules (N magenta module = 113; N red module = 153; N pink module = 128; N yellow module = 180) were exported for further analysis on their relationship with biochemical indexes of osteoporosis (Table S2). Osteoporosis, as an aging-related disease, affects many people and is considered as a societal problem of health (Kastner et al., 2010). Gut microbiota colonizes gastrointestinal tract and plays a key role in bone mass mediation through regulating several biological processes such as intestinal calcium absorption and osteoclastogenesis (Sjögren et al., 2012;Weaver and Connie, 2015;McCabe et al., 2015). In this study, we analyzed the gut microbiota in elderly postmenopausal osteoporotic women and healthy controls. The results showed that there was no obvious difference in Alpha diversity in the experiment group and control group, while Beta diversity showed significant difference. Moreover, the osteoporotic women and healthy controls had different genera of Erysipelotrichaceae, Rothia, and Eubacterium, indicating that these genera might be associated with the occurrence and development of osteoporosis.
Erysipelotrichaceae belongs to the Firmicutes phylum and is implicated in host physiology and diseases, such as metabolic disorders and gastrointestinal diseases (Kaakoush, 2010). It has been proved that Erysipelotrichaceae is associated with the lipid metabolism of the host (Martinez et al., 2013). The species of Erysipelotrichaceae presented a bloom in diet-induced obesity . Obese people showed elevated Erysipelotrichaceae level, while the hamster with hypercholesterolemia had decreased Erysipelotrichaceae level after treatment  weight loss, could suppress the Erysipelotrichaceae growth (Hurt and Wilson, 2012;Etxeberria et al., 2015). On the other hand, accumulating evidence has proved that lipid metabolism is associated with BMD, the key diagnostic index of osteoporosis. Fluctuation of lipid level was found to be related with BMD alteration, while the treatment of hyperlipidemia could lead to elevated BMD (Tian and Yu, 2015). Mundy G et al showed application of statins, which were used to reduce serum cholesterol, could promote bone formation and might be helpful for the treatment of osteoporosis (Mundy, 1999).
Kapitola J et al found women with osteoporosis and vertebral fractures showed increased cholesterol levels compared with the healthy controls (Broulik and Kapitola, 1993). Thus, we speculated that Erysipelotrichaceae might involve in osteoporosis via the lipid metabolism.
Eubacterium belongs to Firmicutes phylum and consists of diverse species (Kadam and Chuan, 2016).
Karlsson FH et al analyzed the gut microbiota in atherosclerosis patients and healthy controls, and found that different Eubacterium abundance between the two groups (Karlsson et al., 2012).
However, to our knowledge, the alterations of Rothia and Eubacterium abundance in osteoporosis were rarely reported. In addition, among the three different genera between osteoporotic patients and healthy controls of our study, both Erysipelotrichaceae and Eubacterium belong to the Firmicutes phylum. Ji YH el al analyzed the gut microbiota in osteoporosis patients and found that Firmicutes showed distinct proportion between osteoporosis patients and negative controls (Wang et al., 2017), which was consistent with our study. These results also indicated that Firmicutes phylum might be potential osteoporosis-related phylum.
In conclusion, our study identified the different genera between postmenopausal osteoporotic women and healthy controls, including Erysipelotrichaceae, Rothia, and Eubacterium. They might be potential targets for osteoporosis. However, further study on the specific effects of these genera on osteoporosis is needed to implicate the possibility and feasibility of them as promising targets for osteoporosis treatment.
Declarations samples. A total of 3,478 OTUs and 2,640 OTUs were identified for the experiment group and control group respectively. There were 2,541 common OUTs between the two groups.   The biochemical indexes in serum and urine of experiment and control groups (serum calcium, serum phosphate, parathyroid hormone (PTH), osteocalcin, serum N-terminal telopeptides of type I collagen (NTX), serum C-terminal telopeptides of type I collagen (CTX), Vitamin D, urine microalbumin, urine creatinine, urine CTX, calcitonin, urine calcium, urine phosphate, blood creatinine, and urine NTX). The osteoclast activity indexes including urine NTX (p < 0.05), serum NTX, and serum CTX in the experiment group were higher than those in the control group, while the serum calcium in the control group was significantly higher than that in the experiment group (p < 0.05).

Figure 5
Clustering analysis of the samples showed there were no outlier samples.