The detailed analyses process was shown in Figure 1. In the present research, we first performed the targeted NGS to the R-spondin/Wnt pathway genes in 400 Chinese participants (228 OP patients) and replicated the findings using SNPsacn technology among another independent cohort of 768 individuals (356 OP patients) in the same area. Bioinformatic analysis was then performed referring to the HaploReg, SNPinfo, 3DSNP, and Genotype-Tissue Expression (GTEx) data to provide more evidence for the OP-genotype associations. Third, the 16S rRNA gene high-throughput sequencing technology was adopted to explore OP-associated gut microbiota variations. Finally, integrative analysis was conducted combining the genotyping and microbiome data to clarify potential genome-microbiome associations in OP pathogenesis.

The study involved a total of 1,168 Han Chinese participants aged ≥60 years recruited at two communities in Wuhan city and the Wuhan Union Hospital during 2016–2018. The exclusion criteria were as follows: (1) with other endocrine diseases (e.g., hyperthyroidism, hypothyroidism, etc.) that influences bone metabolism; (2) of surgical menopause (i.e., hysterectomy and/or ovariectomy); and (3) taking medicines affecting bone health such as hormones. Moreover, individuals were further excluded before stool collection if they met the following criterion: (1) use of antibiotics within 1 month before fecal sample collection; and (2) with prevalent diseases of diabetes and gastrointestinal diseases.

The covariates have been described in our previous study (Li et al., 2020). Briefly, body mass index (BMI) and data of self-reported menstrual history, lifestyle factors, prevalent diseases, and medication history were collected. This study was approved by the Ethics Committee of Tongji Medical College of Huazhong University of Science and Technology. Informed consent was obtained from all participants before enrollment.

Dual-energy X-ray absorptiometry (Lunar Prodigy, GE, United States) was used to obtain BMD measurements, including values of BMD (expressed in g/cm²), T-score, and Z-score, at skeleton sites of the lumbar spine (LS) and total hip among all participants. A T-score ≤ −2.5 indicated prevalent osteoporosis referring to the World Health Organization criteria (Chen et al., 2020). Subjects with T-score ≤ −2.5 were grouped into OP patients group, and those with T-score > -2.5 were classified into the controls.

Details of venous blood sample collection, genome DNA extraction and gene sequencing methods as NGS and SNPscan were available in one of our studies (Li et al., 2020). Briefly speaking, the targeted NGS was used in the discovery stage (N = 400 and 228 OP patients) to reveal the OP-associated genetic variants located in the target regions of R-spondin/Wnt pathway genes including RSPO1, RSPO2, RSPO3, RSPO4, LRP5, LRP6, LGR4, LGR5, and LGR6. Then in the replication stage, the SNPscan^TM Kit (Genesky Biotechnologies Inc., Suzhou, China) was used for sequencing of the observed genetic variants in association with OP risk in an independent cohort (N = 768 and 356 OP patients). Since limited numbers of sample sizes in our study, we only involved common genetic variants (Minor Allele Frequency > 0.05) for analysis in the present study.

The identified OP-associated SNPs and their LD proxies (r² ≥ 0.8) were annotated for potential regulatory function by HaploReg v4.1 (Ward and Kellis, 2016). The FuncPred tool was used for functional SNPs prediction in the SNPinfo web server (Xu and Taylor, 2009). Three-dimensional (3D) chromatin looping data (Lu et al., 2017) were used to link promising SNPs to their three-dimensional interacting genes. Expression quantitative trait loci (eQTL) and splicing quantitative trait loci (sQTL) analysis were conducted on the SNPs using the GTEx project (Carithers and Moore, 2015).

The high-throughput 16S rRNA gene sequencing method was used to detect the microbiota of fecal samples collected from eligible participants in this study. Details of fecal sample collection, DNA extraction, and PCR amplification were described in one previous study (Li et al., 2019). Briefly, fecal samples were collected and stored at −80°C until further processing. The fecal microbial DNA was extracted using the QIAamp DNA Stool Mini Kit (Qiagen, Hilden, Germany). PCR amplification was carried out to the V3–V4 hypervariable region of the 16S rRNA gene with bacterial genomic DNA as template. Then, the PCR products were pooled and purified with Agencourt AMPure XP magnetic beads (Beckman Coulter, United States) using the TopTaq DNA Polymerase kit (Transgen, China). PCR products were sequenced using the Illumina Miseq platform with the 2 × 250 bp paired-end method after the library was quantified, mixed, and quality checked. To obtain clean reads, we firstly used TrimGalore¹ to filter raw reads at Q20 and adapter sequence and removed the reads with length <100 bp. Then, FLASH2 (Magoč and Salzberg, 2011) was used to merge pairs of reads from the original DNA fragments. We further removed the low-quality sequences after merging. After that, Mothur (Schloss et al., 2009) was used to remove primers in sequence, and sequences including N-base/homopolymer > 6 bp. Reads with an error rate of >2 and reads with length <100 bp were removed using USEARCH² to obtain clean reads for further analyses. Operational taxonomic units (OTUs) were assigned by clustering the sequences with a threshold of 97% pairwise identity and chimeras were removed using UPARSE (Edgar, 2013). OTUs were taxonomically assigned at a confidence threshold of 80% based on the Ribosomal Database Project database by Mothur.

The Pearson χ² test was used to analyze the between-group differences of categorical variables (e.g., sex, smoking, and drinking), which were expressed as percentages. Two independent-sample t-test was used to analyze the between-group differences of variables of a normal distribution, which were expressed as mean ± standard deviation (SD). Mann–Whitney U test was used to analyze the between-group differences of continuous variables not conforming to a normal distribution, which were presented as median [interquartile range (IQR)]. To assess the beta diversity, a principal coordinate analysis (PCoA) based on the weighted UniFrac was performed for the visualization of the microbiome structure separation across groups. Statistical significance was confirmed using Permutational MANOVA (PERMANOVA). The linear discriminant analysis (LDA) effect size (LEfSe) was performed to identify the species accounting for the differences between groups, where the species with LDA values > 2.5 at a P-value < 0.05 were considered significantly enriched. Logistic regression modeling was performed to assess the associations of sequenced genetic variants with OP risk. Spearman’s correlation analysis was used to investigate the correlation between SNPs and gut microbiota. A generalized linear model based on negative binomial distribution was conducted to explore the relationship between SNPs and gut microbiota by taking the attributes of the microbiome (like the relative abundance of the microbiome) as the dependent variables, the genotypes in the host loci dominant model as the independent variables, and the gender, age, BMI, smoking, and drinking as covariates. The R software and SPSS version 23.0 were used for statistical analyses. A two-sided P < 0.05 was considered to achieve statistical significance.

The demographic characteristics of all participants were presented in Table 1. In the discovery and replication stages, the mean age of individuals in the cases with prevalent osteoporosis (OP) was slightly younger than the controls with normal BMD values. Taking the idea that women have high risk of OP prevalence into account, we involved all 400 females (228 cases and 172 controls) in the discovery stage of genotyping. And in the replication stage among 768 individuals, we also included males. It seems that the OP patients were likely to have larger BMI compared to controls. For the women studied, the cases and controls had similar mean menarche age (about 14 years old), as well as menopause age (about 49 years old). Generally speaking, the cases had higher proportions of smoking, and alcohol drinking, while lower proportions of self-reporting history of fracture and osteoarthritis compared to the controls though the differences achieved significance or not.

Four genetic models (i.e., the additive, dominant, co-dominant, and recessive models) were employed to analyze the association between the studied SNPs and OP risk, results of which were presented in Table 2. Combing the data from both stages, it was observed in the additive model that the rs10920362 T allele carriers had a significant increased risk of OP prevalence compared to the C allele carriers (OR = 1.78, 95% CI: 1.44–2.22, P-FDR = 1.19 × 10^–6). The dominant model revealed that the CT/TT genotypes were associated with larger OP risk than the CC genotypes (OR = 1.86, 95% CI: 1.44–2.41, P-FDR = 2.50 × 10^–5). And compared to the CT/CC genotypes, the TT genotype was significantly associated with an increased risk of OP prevalence (OR = 3.06, 95% CI: 1.71–5.49, P-FDR = 8.43 × 10^–4). As to the rs11178860, an increased risk of OP was observed in individuals with the A allele compared with the individuals carrying the G allele (OR = 1.48, 95% CI: 1.23–1.78, P-FDR = 1.51 × 10^–4). Importantly, the dominant and recessive models also showed results achieve significance. Moreover, the two common SNPs were also associated with the values of BMD measurements at skeleton sites of LS and total hip (Supplementary Table 1).

Referring to the predictions by HaploReg v4.1, the variant rs11178860 was a proxy (r² > 0.8) of genetic variants locating in a variety of regulatory elements of “Promoter histone marks,” “Enhancer histone marks,” “Motifs changed,” “GRASP QTL hits,” and “Selected eQTL hits” (Supplementary Table 2). Besides, the SNP rs10920362 was predicted to function in regulatory elements of “DNAse,” “Enhancer histone mark,” “Proteins bound,” “Motifs changed,” “GRASP QTL hits,” and “Selected eQTL hits” (Supplementary Table 2). According to the functional analysis by SNPinfo, the rs10920362 was identified as “possibly damaging” with features of an exonic splicing enhancer or exonic splicing silencer (Supplementary Table 3). Transcription factor binding sites (TFBSs) were identified for the rs2304269 of strong LD with the rs11178860 (r² = 0.814 in CHB) (Supplementary Table 3). The 3D chromatin looping data in bone demonstrated that the interacting genes of rs10920362 are LGR6 and SYT2, and rs10920362 was linked with 14 TFBSs (Supplementary Figure 1).

Based on the GTEx datasets, we analyzed the correlation between the two genetic variants (rs10920362 and rs11178860) and gene expressions in various tissues, results of which are shown in Supplementary Table 4. We observed potential associations between rs10920362 and the LGR6 gene expression in tissues of thyroid (P = 7.20 × 10^–39), brain hypothalamus (P = 1.10 × 10^–11), minor salivary gland (P = 2.10 × 10^–9), artery tibial (P = 1.00 × 10^–8), brain substantia nigra (P = 3.00 × 10^–7), artery aorta (P = 4.20 × 10^–7), and brain cortex (P = 4.30 × 10^–5). We failed to find the rs11178860 and LGR5 gene expression association, but the sQTL data showed a potential association between rs10879301 (of large LD with rs11178860: r² = 1.00 in CHB) and the splicing changes of the LGR5 gene in the muscle-skeletal tissue (Supplementary Figure 2).

A total of 180 fecal samples from 77 OP patients and 103 controls were included in the gut microbiota abundance analysis (Supplementary Table 5). The OP group had comparatively smaller numbers of bacterial taxa at species, genus, family, order, class, and phylum levels than the control group (Supplementary Table 6). A total of 1,699 OTUs were quantified in this study, specifically, 1,288 in the OP group and 1,556 in the control group with 1,145 shared in both groups.

The PCoA of weighted UniFrac distance showed a significant difference between the OP and control groups (Supplementary Figure 3). The PC1 explained 46.64% of variation and the PC2 explained 10.53%. The PERMANOVA analysis also revealed that the two groups had a significant difference in beta diversity (F = 3.413, P = 1.0 × 10^–4). LEfSe analysis revealed that the Bacteroidetes (phylum), Porphyromonadaceae (family), Bacteroidaceae (family), Parabacteroides (genus), and Bacteroides (genus) were enriched in the OP individuals; and the Firmicutes (phylum), Proteobacteria (phylum), Actinobacteria (phylum), Lachnospiraceae (family), Bifidobacteriaceae (family), Butyricicoccus (genus) and Bifidobacterium (genus) were enriched in the controls (Supplementary Figure 4). The Firmicutes, Bacteroidetes, Proteobacteria, and Actinobacteria mainly dominated at phylum level in gut microbiota of all samples. The proportion of Bacteroidetes was significantly larger in the OP group than that in the control group (P < 0.05), while those of the Firmicutes, Proteobacteria, and Actinobacteria were smaller in the former than the latter (all P < 0.05) (Table 3). At family level, the OP patients had a higher proportion of Bacteroidaceae and Porphyromonadaceae than the controls (both P < 0.05). At genus level, the Bacteroides and Parabacteroides accounted for a larger proportion in the case group than the controls (both P < 0.05). Spearman’s correlation analysis observed that levels of the Firmicutes, Proteobacteria, and Actinobacteria phylum, the Bifidobacteriaceae, Lactobacillaceae, Ruminococcaceae, and Ruminococcaceae families, and the Bifidobacterium, Lactobacillus, and Gemmiger genus correlated positively with the BMDs and T-scores among all subjects (Supplementary Table 7). Of note, the ridge regression model still showed significant associations between gut microbiota and skeleton site-specific BMD measurements after adjusting for age, sex, smoking, and BMI (Supplementary Table 8).

A total of 113 participants (53 OP patients) simultaneously with genetic and gut microbiota data were included in this part. A total of 1,521 OTUs were identified. The Venn diagram showed that the rs10920362 CC and CT/TT carriers shared 973 (64.0%) OTUs with 348 found only in the CC carriers and 200 found only in the CT/TT carriers (Figure 2). The CC carriers had more OTUs (1,321) than the CT/TT carriers (1,173). A total of 864 (56.8%) OTUs were shared across the rs11178860 GG and GA/AA genotypes with 569 OTUs only found in the GG genotype individuals and 88 found only in the GA/AA genotype individuals. The number of OTUs was larger in the GG genotypes (1,433) than the GA/AA genotypes (952).

LEfSe analysis revealed that the bacterial species belonging to phylum Actinobacteria, family Bifidobacteriaceae, and family Streptococcaceae were enriched in the LGR6 rs10920362 CC genotype. Compared to the LGR6 rs10920362 CT/TT genotype individuals, CC genotype individuals had a higher level of Bifidobacterium and Streptococcus (Figure 3). As for rs11178860, no significant differences were found.

As data shown in further analyses, the proportion of Actinobacteria in the rs10920362 CT/TT genotype was significantly lower than that in the CC genotype at the phylum level (P = 0.005) (Table 4). At the family level, the individuals carrying CT/TT genotype had lower abundances of Bifidobacteriaceae compared to those with CC genotype (P = 0.002). At the genus level, individuals with CT/TT genotype had lower abundances of Bifidobacterium than individuals with CC genotype (P = 0.002). However, no significant differences in gut microbial community abundance were found between rs11178860 GG and GA/AA genotypes. Of note, the abundances of Actinobacteria, Bifidobacteriaceae, and Bifidobacterium were also positively correlated with BMD and T-score (Li et al., 2019).

Spearman’s correlation analysis showed that the abundances of Actinobacteria, Bifidobacteriaceae, and Bifidobacterium correlated negatively with the rs10920362 CT/TT genotype (all P < 0.05) (Table 5). No significant correlation was found between SNP rs11178860 and gut microbiota abundance.

Based on above results, we took the relative abundance of Actinobacteria, Bifidobacteriaceae, and Bifidobacterium as dependent variables, and the rs10920362 and rs11178860 genotypes as independent variables in the generalized linear modeling. Compared with the rs10920362 CC genotype, the CT/TT genotype associated with decreased relative abundance of Actinobacteria (β = −1.32, P < 0.001), Bifidobacteriaceae (β = −1.70, P < 0.001), and Bifidobacterium (β = −1.70, P < 0.001) controlling for sex, age, BMI, cigarette smoking, alcohol drinking, and self-reporting history of osteoarthritis (Table 6). No significant associations were found between the rs11178860 genotype and the relative abundance of above-mentioned gut microbiota.