GWAS data sets
BMD measured by dual-energy X-ray absorptiometry (DXA) is the gold standard for OP diagnosis as recommended by the WHO. Since LS BMD is estimated to have the highest heritability among commonly measured skeletal sites [16] it was thus chosen for studying the relationship between BMD and BW in this genetic association analysis.
The LS BMD GWAS summary-statistic dataset was acquired from the Genetic Factors for Osteoporosis Consortium (GEFOS), and included of 53,236 individuals of European ancestry [17]. At the present time, it is the largest GWAS for DXA derived BMD measurements. The BW summary statistics were downloaded from the Early Growth Genetics (EGG) Consortium, and included 143,677 subjects of European ancestry [12]. To confirm that the variance estimated for each single-nucleotide polymorphism (SNP) was not inflated due to population stratification, standard genomic control procedures were applied to the two original GWAS studies.
SNP pruning and merging
We performed linkage disequilibrium (LD)-based pruning for each dataset using PLINK version 1.9 software. Firstly, the LD was computed for each pair of SNPs in a window containing 50 SNPs. For pairs with an r2 > 0.2, the SNP with the smaller minor allele frequency was removed. The window was moved 5 SNPs forward, and the procedure was repeated until no pairs of SNPs across the genome had r2 > 0.2. The pruning was based on the LD structure of the CEU HapMap 3 genotypes. After the pruning process was completed, there were 121,848 SNPs which overlapped between BMD and BW which were retained for the subsequent analysis.
Estimation of pleiotropic enrichment
Stratified Quantile-quantile plots (Q-Q plots) were constructed to visualize the pleiotropic enrichment between LS BMD and BW when conditioning on successively more stringent p-value thresholds of the conditional trait: p < 1 (all SNPs), p < 0.1, p < 0.01 and p < 0.001. The observed p-values of the principal trait, denoted as “nominal -log10 (p)”, were plotted against the cumulative distribution function of observed p-values, denoted as “empirical − log10 (q)”. The line x = y indicates the null hypothesis of no pleiotropic enrichment, and plots that deviate leftward from the null line indicate that a pleiotropic effect exists between the traits[13].
Calculation of cFDR and conjunction cFDR (ccFDR)
The calculation of the cFDR extends from the single phenotype case, where the unconditional false discovery rate (FDR) for a set of SNPs is characterized as the probability of a false positive association. The cFDR expands this idea to the two-phenotype case and is defined as the probability of a false positive association with the principal trait given that the association p-values with both the principal and conditional traits are at least as small as the observed p-values. Using the GWAS summary statistics, the cFDR for each SNP was separately computed for both orderings of the traits (BMD|BW and BW|BMD, where “|” indicates conditional upon). The SNPs were regarded as significantly related to the principal trait when cFDR < 0.05. The ccFDR value was defined as the maximum cFDR value of two trait orderings, and SNPs with ccFDR < 0.05 were interpreted as pleiotropic loci associated with both traits. As the cFDR method has been widely published and successfully utilized both by our team [10] and other teams [14, 15], the specific procedure and formulae for calculation will not be detailed in this article. Manhattan plots were constructed using R to graphically display the genomic locations of significant variants.
Annotation of novel SNPs and genes
The significant cFDR SNPs (cFDR < 0.05) were queried using the SNPinfo web server (https://snpinfo.niehs.nih.gov) to ascertain all corresponding SNPs with high LD (r2 > 0.8). All SNPs (including those with cFDR < 0.05 and those with high LD) were then compared to previous GWAS findings (p-value < 5 × 10− 8) on the European Bioinformatics Institute website. SNPs with a GWAS p-value > 5 × 10− 8 that have not been reported as having an association with BMD or/and BW were considered to be potential novel SNPs. We utilized both the SNP and CNV Annotation Database (SCAN, http://scandb.org/newinterface/about.html) and PubMed (https://www.ncbi.nlm.nih.gov/gene) to map the significant cFDR SNPs to nearby genes. Genes that were not previously identified in BMD or/and BW-related studies were deemed novel.
Functional enrichment analysis and protein–protein interaction analysis of identified genes
To establish the physiological role of genes of interest, we performed functional enrichment analysis using the Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/summary.jsp). To explore the functional interactions between proteins produced by cFDR-significant pleiotropic susceptibility genes, we performed protein–protein interaction analysis by the online tool STRING 10.0 (http://string-db.org/).
Fine-mapping
To discover putative pleiotropic causal SNPs and prioritize genes for the subsequent functional experiments, we performed multi-trait fine mapping analysis using PAINTOR[18]. In particular, we focused on a potential OP susceptibility locus located on chromosome 16 (14800000-16280000) which includes the gene PDXDC1.
Postmenopausal OP mouse models
Female C57BL/6 mice were obtained from the Animal Center of Southern Medical University, which either received sham-operated (SHAM) surgery or ovariectomy (OVX) under 1.2% tribromoethanol anesthesia at 8 weeks of age. The OVX group had bilateral ovary removal, while the SHAM group had a similar volume of adipose tissue removed from around the ovaries. The mice were sacrificed sixteen weeks post-surgery and hind-limb specimens were harvested for subsequent analyses. Four mice were randomly selected from each experimental group.
Micro-computed tomography analysis
The hindlimbs were fixed in 4% paraformaldehyde for 48 h then scanned in a micro-computed tomography (micro-CT) scanner (Viva CT40, Scanco Medical AG, Bassersdorf, Switzerland). Morphological analysis was conducted on trabecular bone. The region of interest in the trabecular bone began at a position 20 spongiosa slices (9 µm thick) from the lower growth plate of the femur and finished 160 slices later. Bone volume/total volume (BV/TV), trabecular thickness (Tb. Th), trabecular number (Tb. N), and trabecular separation (Tb. Sp) were computed using standard three-dimensional microstructural analysis.
Immunohistochemistry
The hindlimbs were fixed in 4% paraformaldehyde in phosphate-buffered saline (PBS) for 48 h at 4°C, then decalcified in 10% ethylenediaminetetraacetic acid (EDTA; pH 7.4) for 21 days at room temperature prior to dehydration in a rising gradient of ethanol and embedded in paraffin. The tissues were sliced into 3 or 4 micron thick sections for histological analysis. Tissue sections were incubated overnight at 4°C with PDXDC1 primary antibody (Proteintech, 21021-1-AP, 1:200) then for 1 h at room temperature with a secondary antibody (Arigo, ARG65351, 1:200). Diaminobenzidine (DAB, KGP1045/KGP1045-20/KGP1045-100, 1:1:1:20) was conjugated to the secondary antibody.
Both PDXDC1-positive and the total number of cells were enumerated (at 400× magnification) at the growth plate and within the trabecula bone in the femur or tibia. Four views from these regions picked at random were counted on each section, and three consecutive sections were selected for each mouse.
To distinguish whether the PDXDC1-positive cells were osteoclasts or osteoblasts, we conducted immunohistochemistry assay with three-micrometer-thick serial sections in the SHAM group. Tartrate-resistant acid phosphatase (TRAP) is a specific marker of osteoclasts and osteocalcin (OCN) is a specific marker of osteoblasts. The specific procedures of the immunohistochemistry assay with serial sections were as follows: 1) there were two thin and consecutive tissue sections on each glass slide labeled A and B in advance, 2) A-labeled sections were incubated with PDXDC1 primary antibody (Proteintech, 21021-1-AP, 1:200). B-labeled sections were incubated with anti-OCN primary antibody (abcam, ab93876, 1:500) or subjected to TRAP staining (Sigma-Aldrich), 3) sections not used for TRAP staining were incubated with suitable secondary antibody (Arigo, ARG65351, 1:200).
fat-1 transgenic (TG) mouse model
To preliminarily verify the mechanism of PDXDC1 affecting BW and BMD, a fat-1 TG mouse model was adopted. Male fat-1 TG mice were matched with wild-type (WT) female C57BL/6 mice to breed fat-1 gene-positive mice. The fat-1 gene-positive mice were identified using genomic DNA extracted from tail biopsies. Primer sequences were as follows: 5’-GGACCTGGTGAAGAGCATCCG-3’ and reverse, 5’-GCCGTCGCAGAAGCCAAAC-3’. We fed the mice until 16 months old when they were sacrificed and hind-limb specimens were harvested for subsequent analysis. We performed the immunohistochemical assay on both fat-1 TG mice and WT mice with PDXDC1 primary antibody (Proteintech, 21021-1-AP, 1:200).
Sections were imaged using a Zeiss microscope (Carl Zeiss, New York, USA). All experiments were performed three times or more for reproducibility.
Statistics
Data shown are means ± standard deviation (SD). Two sample t-tests were performed for comparison of experimental groups. Statistical significance was set at p-value < 0.05.
Mendelian randomization (MR) analysis
To investigate the potential causal relationship between BW and LS BMD, we performed two-sample MR analysis using the BW-associated SNPs as instrumental variables. [19]. The causal effects from multiple instruments were combined using several meta-analysis approaches including maximum likelihood estimation and inverse-variance weighting (IVW).