Characterization of Intestinal Flora in Osteoporosis Patients Based on 16S rDNA Sequencing

Introduction

Osteoporosis (OP) is a prevalent systemic bone disease characterized by reduced bone mass, damage to bone microstructure, increased fragility, and an elevated susceptibility to fractures. An osteoporotic fracture, occurring as a complication of OP, contributes significantly to mortality in the elderly, with a hip fracture-related mortality rate ranging from 20% to 30% within the first year. In addition, there is a critical deficiency in the treatment of OP patients after fracture, which leads to a great probability of re-fracture and high treatment costs.^1,2 Hence, intervening early during the progression of OP can alleviate the strain on families and society, enhance the quality of life for patients, and extend their long-term survival.

The gut microbiota consists of trillions of microorganisms encoding 150 times more genes than the human genome and has thus been called the “second largest human genome”. The intestinal flora plays a crucial role in preserving host homeostasis by facilitating food digestion and absorption, secreting metabolites, protecting the gastrointestinal mucosal barrier, and modulating the immune response.³ Intestinal flora can affect bone metabolism in several ways, including regulating nutrient absorption, modulating the immune system, and releasing metabolites.^4–6

Previous studies on the correlation between gut microbiota and bone metabolism usually used germ-free (GF) mice as a model of impaired gut microbiome function,⁷ or compared mice after antibiotic treatment with mice with normal gut microbiota. However, there are few studies on the comparison of gut flora characteristics between OP patients and the normal population. In this study, high-throughput sequencing of bacteria in fecal samples from healthy populations and OP patients was performed using 16S rDNA sequencing. Subsequently, a series of analyses were carried out based on the sequencing results. In summary, the outcomes of this study establish a theoretical foundation for investigating the commonality of gut microbiota among OP patients. Furthermore, these findings allow for hypotheses regarding the pivotal bacterial strains linked to OP development, opening avenues for novel strategies in OP treatment. This study was approved by the Ethics Committee of Huai’an First People’s Hospital, with the ethical approval number KY-2019-083-01, the study complies with the Declaration of Helsinki.

Material and Methods

Research Subjects

According to the inclusion and exclusion criteria, the fecal specimens and clinical data of 43 osteoporosis patients hospitalized in the Department of Endocrinology and Metabolism of Huai’an First Hospital Affiliated to Nanjing Medical University from October 2019 to May 2020 were collected and analyzed. At the same time, according to the principle of age and gender matching, 24 healthy people who underwent physical examination in our hospital were recruited as the control group. Inclusion criteria: A. Osteoporosis group: 1. In line with the latest internationally recommended diagnostic criteria for osteoporosis by WHO: postmenopausal women or men over 50 years old, with a bone mineral density T value of −2.5 or less measured by dual energy X-ray absorptiometry (DXA), or patients with fragile fractures of vertebral body and hip, or fragile fractures of proximal humerus, pelvis or distal forearm, and with bone loss measured by DXA (−2.5<T value≤-1); 2. Residents of Huai’an city. B. Healthy population group: 1. The bone mineral density T value measured by DXA>-1; 2. No previous history of fragile fractures; 3. Residents of Huai’an city.

Exclusion criteria: 1. Secondary osteoporosis caused by various causes; 2. A variety of acute and chronic diseases that can lead to gastrointestinal dysfunction, such as dysentery and acute gastroenteritis; 3. After various gastrointestinal diseases with intestinal loss and gastrointestinal anatomical changes; 4. Patients with hemorrhoids, diarrhea, and other factors that affect the morphology and main composition of feces within 4 weeks; 5. Patients who received antibiotic therapy within 6 weeks; 6. Patients who used probiotics or prebiotics within 2 weeks; 7. Patients who received long-term non-steroidal anti-inflammatory drug therapy; 8. Patients with immune deficiency.

Specimens and Data Collection

Clinical Data

The participants were provided with detailed information regarding the research objectives, content, and methodologies of this study, and they subsequently provided their informed consent by signing the necessary documentation. The age, gender, height, and weight variables of both groups were subjected to statistical analysis. All participants underwent dual-energy X-ray absorptiometry (DXA) testing to assess bone mineral density (BMD), with measurements recorded for the overall L1-4 lumbar spine region, left femoral neck area, and overall hip region.

Stool Collection

The aseptic fecal collection device was utilized to obtain 10g of freshly solid feces from the subjects in their first morning defecation. Subsequently, the collected fecal samples were divided into three portions and stored in 1.5mL aseptic cryopreserved tubes within a laboratory-grade −80°C refrigerator. Following the acquisition of fecal specimens from both subject groups, intestinal flora’s 16S rDNA sequencing and subsequent bioinformatics analysis were performed.

Sequencing of Intestinal Flora

DNA Extraction and Assay

Fecal DNA was extracted utilizing the Stool Extraction DNA Kit (OMEGA bio-tek, USA). The quality of extracted DNA was measured by agarose gel electrophoresis. In addition, DNA was quantified using a UV spectrophotometer.

PCR Amplification and Purification

The bacterial 16S rDNA V3-V4 region was selected as the target for PCR amplification. The following amplification primers were used: 341F (5’-CCTACGGGGNGGCWGCAG-3’) and 805R (5’-GACTACHVGGGGTATCTAATCC-3’). The reaction procedure was as follows: 98°C for 30s; 98°C for 10s; 54°C for 30s (35 cycles); 72°C for 45s; and 72°C for 10 min. PCR products were confirmed by 2% agarose gel electrophoresis, purified by AMPure XT beads (Beckman, USA), and quantified by Qubit dsDNA HS Assay Kit (Invitrogen, USA).

The size and number of amplicon libraries were evaluated using an Agilent 2100 Bioanalyzer (Agilent, USA) and a library quantification kit from Illumina (Kapa Biosciences, USA), respectively, and sequenced on the Illumina NovaSeq platform after meeting the standards.

Bioinformatics Analysis

Following the completion of on-board sequencing and the acquisition of raw data, the paired-end data were initially spliced using overlap. The spliced data were then subjected to quality control, chimera filtering, and a series of other processes to finally obtain high-quality valid (clean) data. The final validated data obtained were subjected to the class operational taxonomic unit (OUT) analysis using cluster analysis and species taxonomic analysis DADA2 (Divisive Amplicon Denoising Algorithm) to obtain the feature table and feature sequence. Subsequently, species composition was further described, species diversity was analyzed, and subsequent dissimilarity analyses were carried out.Subsequently, further analyses of species composition, species diversity, and differences were performed (Figure 1).

Figure 1 Information analysis flow.

Results

General Information

As presented in Table 1, no significant difference in gender, age, and body mass index (BMI) was observed between the OP and the healthy control categories (P > 0.05).

Table 1 Clinical Characteristics of the Study Population

Species Composition Analysis

The Venn diagram of feature distribution was plotted based on the feature list obtained from clustering as well as the feature value abundance table (Figure 2). Additionally, 1277 features common to both categories, 2995 unique to the OP category, and 2155 unique to the healthy controls were recorded.

Figure 2 Venn diagram: each circle in the figure represents a group, the number of overlapping parts represents the number of features shared between the groups, and the number without overlapping parts represents the number of features unique to the group.

Species Diversity Analysis

The diversity of the flora in the samples was further characterized. As shown in Figure 3, the dilution curves of both categories gradually flattened with increasing horizontal coordinates, confirming that the samples reached a reasonable sequencing depth. In addition, species were richer in the healthy population category than in the OP category at the same sequencing level. No significant difference was observed between the two categories in alpha diversity indices, including Chao1, Goods_coverage, Simpson, and Shannon (Table 2), confirming no significant differences in diversity, richness, and homogeneity of the intestinal flora between the two categories.

Table 2 Alpha Diversity Indices

Figure 3 Dilution curve, which is used to reflect the sequencing depth (Blue, OP; Red, Control).

Histogram Analysis of Flora

According to the order of relative abundance of flora, the distribution of flora at the phylum level in the OP category and the healthy population category was dominated by Firmicutes, Actinobacteria, Proteobacteria, and Bacteroidetes (Figure 4A). The top ten flora in relative abundance at the genus level were Bifidobacterium, Faecalibacterium, Escherichia-Shigella, Subdoligranulum, Bacteroides, Streptococcus, Agathobacter, Akkermansia, collinsella, and Roseburia (Figure 4B).

Figure 4 Cylindrical accumulation Map of relative abundance of species at Phylum level (A) and Genus level (B). The abscissa is grouping information; the ordinate represents Relative Abundance; others represents the sum of the relative abundances of all the phylums except the 10 phylums in the figure.

Between-Group Species Difference Analysis

Differential strains of intestinal flora were present at all levels between the two categories. There were four statistically different flora at the phylum level. Among the top five flora in terms of relative abundance, the abundance of Verrucomicrobia and Proteobacteria was higher in the OP category than in the controls, while the abundance of Actinobacteria was lower than in the controls (Table 3). A total of 33 significant differential flora were observed at the genus level. Among the top ten flora in terms of relative abundance, Akkermansia, Klebsiella, Megasphaera, Escherichia-Shigella, and Lactobacillus had lower abundance in the healthy population category than in the OP category. Additionally, Faecalibacterium, Ruminococcaceae_UCG-009, Alloprevotella, Clostridiales_unclassified, Parasutterella, Lachnospiraceae ND3007_group, and Negativibacillus had higher abundance in the healthy population than in the OP category (Table 3).

Table 3 Differential Flora Between the Two Categories at Phylum and Genus Levels

LEfSe and LDA Analyses

From the perspective of evolutionary branching, the differential flora that may play a role were Verrucomicrobia-Verrucomicrobiae-Verrucomicrobiales-Akkermansia; Gammaproteobacteria-Enterobacteriales-Enterobacteriaceae-Klebsiella; Enterobacteriales-Enterobacteriaceae-Escherichia-Shigella; and Lactobacillaceae-Lactobacillus (Figure 5). According to the LDA analysis, the flora with significantly higher abundance in the OP category were Enterobacteriales, Gammaproteobacteria, Enterobacteriaceae, Akkermansia, Proteobacteria, Escherichia-Shigella, Verrucomicrobiales, and Klebsiella. However, the predominant species in the intestinal samples of the healthy population category were actinobacteria, Faecalibacterium, and Subdoligranulum (Figure 6).

Figure 5 LEfSe analysis:Each small circle at different classification levels represents a classification at that level, and the diameter of the small circle is proportional to the relative abundance. The coloring principle: the species with no significant difference were uniformly colored yellow, the differential species Biomarker followed the group, and the red nodes represent the microbial group which played an important role in the red group. Green nodes represent microbial groups that play an important role in the green group.

Figure 6 The LDA value distribution histogram shows the species with an LDA score greater than the set value (the default setting is 4), that is, the Biomarker with statistical differences between groups. The length of the histogram represents the impact of the different species (LDA score).

Correlation Analysis of Clinical Data and Flora

To further validate the key flora that may affect bone mineral density (BMD), correlation analyses of differential flora and clinical data were performed at phylum and genus levels (Figures 7 and S1). At the phylum level, the abundance of Proteobacteria was negatively correlated with the BMD of the femoral neck and the whole hip (P < 0.01) and BMI (P < 0.05); the abundance of Verrucomicrobia was negatively correlated with the BMD of the whole L_1-4 lumbar vertebrae (P < 0.05). At the genus level, Akkermansia abundance was negatively correlated with the BMD of the whole L_1-4 lumbar vertebrae (P < 0.05); Klebsiella abundance was negatively correlated with the BMD of the femoral neck, the whole hip, and the whole L_1-4 lumbar vertebrae (P < 0.05); Escherichia-shigella abundance was negatively correlated with the BMI (P < 0.05) as well as the BMD of the femoral neck, the whole hip, and the whole L_1-4 lumbar vertebrae (P < 0.01). Faecalibacterium abundance was positively correlated with the BMD of the whole hip and the whole L_1-4 lumbar vertebrae (P < 0.01); Megasphaera abundance was negatively correlated with BMI (P < 0.05).

Figure 7 Correlation analysis of flora abundance and clinical data:Spearman correlation analyses were performed to explore the correlation between differential flora in each level and the collected clinical data (including general information and laboratory examination results).Color represents the correlation coefficients between the two variables, where red indicates a positive correlation and blue a negative correlation; The depth of the color represents the strength of the correlation, the darker the color, the stronger the correlation, and the vice versa, the weaker the correlation; *P < 0.05, **P < 0.01.

Identification of Patients with OP Based on Intestinal Flora

To determine the potential of difference in the gut flora composition as a biomarker for differentiating OP patients from healthy populations, the dominant flora in the OP category at the genus level was selected to plot ROC curves to assess their discriminatory ability. As presented in Figure 8, the AUGs of Klebsiella, Akkermansia, and Escherichia-Shigella all ranged from 0.7 to 1.0, suggesting their certain diagnostic value for OP at appropriate thresholds.

Figure 8 ROC curves: abscissa: specificity scale, ordinate: sensitivity scale.ROC, receiver operating characteristic; AUC, area under the ROC curve; CI, confidence interval.

Prediction and Analysis of Genomic Function

To explore the potential mechanisms via which flora affects bone metabolism, the gene function of flora genes within the samples was predicted using two major databases, KEGG and COG. The COG functional prediction results (Figure 9) showed that the OP category had significantly more homologous genes in translational, ribosomal structural, and biogenesis signaling pathways than the healthy control category, but had significantly fewer homologous genes in signaling pathways related to inorganic ion transport and metabolism, intracellular transport, secretion and vesicular transport, and extracellular structure (P < 0.05). KEGG functional prediction results (Figure 10) showed that the OP category had significantly more homologous genes than the healthy controls in signaling pathways related to carbohydrate metabolism, unfavorable characteristics, and metabolism, but had significantly more homologous genes in signaling pathways related to cellular motility, excretory system, nucleic acid metabolism, translation, and infectious disease than the healthy population category.

Figure 9 COG functional prediction results.

Figure 10 KEGG functional prediction results. Based on the functional annotation results of COG and KEGG database, STAMP difference analysis was used to predict the function of different species between the two groups, in which color represented different groups (OP: yellow, control: blue). The Figure 9 and 10 show the functional prediction results with significant differences in level 2.

Discussion

Gut microbiota is closely related to the physiological activities of the host such as nutrient absorption, immunomodulation, and energy metabolism. In recent years, the role of intestinal flora in disease onset and progression has become a prominent and highly discussed topic in scientific research. Intestinal flora may affect bone metabolism through various pathways such as regulating the host immune system, controlling inflammation, and influencing calcium absorption and vitamin D levels. Therefore, regulating intestinal flora is expected to be a novel approach and strategy to prevent and treat bone loss, OP, and other bone metabolism-related diseases in the future.

The process of in vitro culture of bacteria is affected by a variety of complex factors, requiring a high level of stability in both the growth environment and experimental conditions. This process is time-consuming. Traditional experimental techniques can only perform qualitative and quantitative experiments on known bacterial strains, resulting in a very limited number of bacteria detected. The 16S rDNA sequencing technology, characterized by its culture-free nature, rapid output, and extensive data capabilities, plays a vital role in precisely identifying bacterial isolates and uncovering novel bacterial species. In this study, high-throughput sequencing of bacteria from 43 fecal samples of OP patients and 24 fecal samples of healthy people was performed using 16S rDNA sequencing. The results of microbe diversity, abundance, and compositional distribution were obtained. Subsequently, their correlation with BMD was explored in depth.

The presence of differential flora between the two categories at all levels was determined by intergroup difference analysis. According to LEfSe and LDA analyses, Proteobacteria, Gammaproteobacteria, Enterobacteriales, Verrucomicrobiales, Enterobacteriaceae, Akkermansia, Escherichia-Shigella, and Klebsiella were the flora that had significantly increased abundance in the OP category. In contrast, the dominant flora in intestinal samples from healthy populations were Actinobacteria, Faecalibacterium, and Subdoligranulum. Subsequently, clinical data such as differential flora between the two categories and BMD were subjected to correlation analysis. At the phylum level, Proteobacteria abundance was negatively correlated with BMD of the femoral neck and the whole hip; Verrucomicrobia abundance was negatively correlated with BMD of the whole L_1-4 lumbar vertebrae. At the genus level, Akkermansia abundance was negatively correlated with BMD of the whole L_1-4 lumbar vertebrae; Klebsiella and Escherichia-Shigella abundance was negatively correlated with BMD of the femoral neck, the whole hip, and the whole L_1-4 lumbar vertebrae; Faecalibacterium abundance was positively correlated with BMD of the whole hip and the whole L_1-4 lumbar vertebrae. To determine the potential of gut flora composition difference as a biomarker for distinguishing OP patients from the healthy population, the dominant flora in the OP category at the genus level were selected to plot the ROC curve to assess its discriminatory ability. The AUG values of Klebsiella, Akkermansia, and Escherichia-Shigella were all between 0.7–1.0, suggesting their certain diagnostic value in identifying OP patients. However, more studies are needed to explore the appropriate critical values. Overall, these flora play a role in bone health.

In this study, Akkermansia had a lower abundance in the healthy population than in the OP category, and its abundance was negatively correlated with the BMD of the whole L_1-4 lumbar vertebrae. In a study by Lawenius et al,⁸ Akkermansia-treated mice showed a reduced bone mass. Mechanistically, Akkermansia-treated mice had increased levels of PTH, which activated the mRNA expression of the calcium transporter Trpv5 in the kidney, increased calcium reabsorption, and promoted bone resorption by increasing RANKL and decreasing OPG expression. Based on the findings of previous animal studies and the data acquired in this research, the decrease in bone mass could be attributed to the enrichment of Akkermansia. Another differential flora that may play a role is the Gammaproteobacteria-Enterobacteriaceae-Klebsiella. This branch of the flora was more abundant in osteoporotic patients than in the healthy population and was negatively correlated with the BMD of the femoral neck and the whole hip and BMI. A previous study⁹ has confirmed that ectopic reproduction of Klebsiella in the gut can be mediated by dendritic cells, which activate TLRs and IL-18 signaling pathways, driving Th1 cell proliferation and inflammatory responses. The immune system plays an important regulatory role in the development of OP, and the inflammatory response induced by an increase in the number of T cells, as well as an increase in pro-osteoclastogenic cytokines, such as TNF-α and RANKL, will affect the homeostasis of bone reconstruction.^10,11 Thus, Klebsiella can indirectly affect bone metabolism through the immune system-inflammatory factor pathway. Therefore, Proteobacteria-Klebsiella plays a certain negative regulatory role in host bone metabolism. In the present study, another dominant flora in the OP category was Escherichia-Shigella. Similar to Klebsiella, enrichment of Escherichia-Shigella can lead to increased levels of circulating inflammatory factors, thereby inducing bone loss.¹² In addition, a decrease in the number of intestinal CD8+ T cells occurs in inflammatory bowel disease due to Escherichia-Shigella.¹³ This process reduces the expression of the 1α-hydroxylase CYP27B1 gene,^14,15 suppressing the conversion of vitamin D into the active 1.25-OH-VD, which in turn affects its multifaceted roles in promoting bone mineralization, modulating immunity, and protecting against pathogenic microbes. Therefore, negative regulation of bone mass by Escherichia-Shigella through influencing the activation process of host VD is also one of the possible mechanisms. Moreover, Faecalibacterium abundance was significantly higher in the healthy population than in the OP category and was positively correlated with the BMD of the whole hip and the whole L_1-4 lumbar vertebrae. Bacillus faecalis exerted a significant inhibitory effect on the expression of the inflammatory factor IL-23, thereby preserving the stability and activation of downstream Th17 cells. This effect can be attributed to the substantial anti-inflammatory effect of Bacillus faecalis metabolite, butyric acid.¹⁶ Consequently, this anti-inflammatory response plays a role in suppressing bone resorption. Butyric acid also promotes osteogenic differentiation in vitro by increasing the activity of ALP, a marker of osteoblast differentiation. Additionally, it promotes the accumulation of Ca2+ mineralized nodules during the induced differentiation of human mesenchymal stem cells through the up-regulation of TNFAIP3.¹⁷ Lucas et al¹⁸ showed a significant increase in bone mass and a significant improvement in inflammation-induced bone loss after the use of butyric acid supplementation treatment in mice. These findings further confirm that butyric acid is an effective regulator in maintaining bone homeostasis, which is consistent with the results of the present study. Therefore, it was hypothesized that the enrichment of Faecalibacterium is a positive factor affecting bone mass.

In addition to analyzing the differential microbial communities, KEGG and COG gene function predictions were conducted using the PICRUSt tool. This analysis revealed significant functional disparities between the two groups, particularly in pathways related to metabolism, cellular motility, and genetic processes. Notably, patients with osteoporosis (OP) exhibited a reduced prevalence of homologous genes involved in inorganic ion transport and metabolism, as well as carbohydrate metabolism, compared to healthy controls. The absorption of inorganic ions, including calcium, magnesium, and zinc, plays a pivotal role in the regulation of bone mineral density (BMD) and the prevention of bone loss.^19,20 Carbohydrates, serving as substrates for gut microbiota, undergo anaerobic fermentation to produce short-chain fatty acids (SCFAs), which may indirectly modulate host bone density by influencing the secretion of endocrine factors such as GLP-1 and IGF-1.^21–25 Furthermore, butyric acid, a metabolite produced by gut microbiota, has the potential to impact bone health by altering osteoclast metabolism. Therefore, dysfunction in carbohydrate metabolism is linked to the development of OP.

Conclusion

In summary, this research elaborated on the compositional characteristics of the intestinal flora in OP patients and healthy people. Differential flora between the two categories at all levels were identified by intergroup difference analysis and LEfSe analysis. Based on the correlation analysis between key differential flora and clinical data such as host BMD, key flora with potential diagnostic value for OP were determined. In addition, this study further analyzed the possible pathways via which intestinal flora influence bone health through KEGG and COG gene functional predictions. In summary, the discoveries made in this research offer novel insights into the prevention and treatment of OP.

Limitations and Prospects

As a cross-sectional investigation, this study lacked a longitudinal assessment of intestinal flora status both before and after the onset or treatment of OP. Additionally, no bone metabolism-related indices were collected and correlated with the flora, limiting further exploration of the role of intestinal flora in regulating bone metabolism. Differences in the composition and structure of gut microbiota among different genders and age groups also exist. We agree that, although our study suggests the potential application of our findings in clinical practice, particularly in screening women at risk for osteoporosis, it is necessary to validate these results with larger sample sizes and functional human studies. We fully recognize the need for further research to confirm the clinical significance of our findings in the prevention and treatment of osteoporosis. Therefore, more prospective studies with stratified groups and larger sample sizes are needed. Moreover, the causal relationship between changes in flora composition and the development of OP should be clarified through in vitro or animal experiments, thereby exploring the potential of fecal microbiota transplantation in OP treatment.