Association between SEMA3A signaling pathway genes and BMD/OP risk: An epidemiological and experimental study

Objective: This study aimed to explore the associations of genetic variants in the semaphorin 3A (SEMA3A) signaling pathway genes, including SEMA3A, NRP1, PLXNA1, PLXNA2 and PLXNA3 with osteoporosis (OP) risk and bone mineral density (BMD) in a Chinese Han older adult population.

Study design and method: A two-stage design was adopted. Total of 47.8kb regions in the 5 genes were sequenced using targeted next-generation sequencing (NGS) technology in the discovery stage, and the discovered OP-related single nucleotide polymorphisms (SNPs) were further genotyped using improved multiple linkage detection reaction technique in the validation stage. Methods of ALP/TRAP staining, real-time fluorescent quantitative PCR, and cell proliferation and apoptosis assays were performed with MC3T3-E1 and RAW 264.7 cell lines to clarify biological effects of observed functional variants in cell lines responsible for bone mass remodeling.

Results: Total of 400 postmenopausal women (211 OP cases) were involved in the discovery stage, where 6 common and 4 rare genetic variants were found to be associated with OP risk. In the validation stage among another 859 participants (417 women, 270 OP cases), the PLXNA2 rs2274446 T allele was associated with reduced OP risk and increased femoral neck (FN) BMD compared to the C allele. Moreover, significant associations of NRP1 rs2070296 with FN BMD/OP risk and of NRP1 rs180868035 with lumbar spine and FN BMDs were also observed in the combination dataset analysis. Compared to the osteoblasts/osteoclasts transfected with the wild-type NRP1 rs180868035, those transfected with the mutant-type had reduced mRNA expression of osteoblastic genes (i.e., ALP, RUNX2, SP7 and OCN), while elevated mRNA expression of osteoclastic genes (i.e., TRAP, NFATc1 and CTSK). Furthermore, mutant NRP1 rs180868035 transfection inhibited osteoblast proliferation and osteoclast apoptosis, while promoted osteoclast proliferation and osteoblast apoptosis in corresponding cell lines.

Conclusion: Genetic variants located in NRP1 and PLXNA2 genes were associated with OP risk and BMD. The NRP1 rs180868035 affects bone metabolism by influencing osteoblasts and osteoclasts differentiation, proliferation and apoptosis.

1 Introduction

Osteoporosis (OP) is a systemic skeletal disease characterized by low bone mineral density (BMD) and deteriorative bone microstructure. Fragility fracture is the most common and harmful complication of OP (1). Dual-energy X-ray absorptiometry (DXA)-derived BMD measurements at skeletal sites of total hip, lumbar spine (LS), and femoral neck (FN) were recommended for the diagnosis of OP (2). A systematic review and meta-analysis study demonstrated that the prevalence of OP was 18.3% worldwide, with a significantly higher prevalence in women (23.1%) than in men (11.7%) (3). It was reported the number of fragility fractures was estimated at 3.5 million in 2010, and it would increase 28%, from 3.5 million to 4.5 million, by 2025 in Europe (4). In a study involving 3157 American participants aged ≥50 years, 30% men and 49% women were diagnosed with osteopenia, and 2% men and 10% women were diagnosed with OP based on the BMD measurements at the FN (5). A recent nationwide epidemiological study in China estimated the age-standardized prevalence of OP was 6.46% and 29.13% for men and women aged ≥50 years, corresponding to 10.9 million men and 49.3 million women (6). Above-mentioned evidence suggests that the prevalence of OP stands at a high level, and OP not only causes high disability and mortality, but brings a heavy burden to society, becoming a serious public health issue.

Bone remodeling is a dynamically balanced process regulated by both osteoblasts and osteoclasts. OP occurs when osteoclast-induced bone mass loss exceeds osteoblast-induced bone mass formation. Semaphorin 3A (SEMA3A) is one of the semaphorin family members that are involved in various physiological and pathological activities in human body, e.g., embryonic and nervous system development, immune regulation, promotion of inflammation reduction and oncogenesis (7–10). Recently, SEMA3A was also found to possess dual role of promoting osteoblastogenesis and inhibiting osteoclastogenesis, suggesting that it may be a potential therapeutic target for OP (11). Animal studies found that the SEMA3A supplementation via transcutaneous injection into the center of the distraction zone elevated BMD, tissue mineral density (TMD) and vascular density of operative area (of tibia) in adult male mice underwent tibia osteotomy surgery (12). Hayashi et al. observed obvious OP phenotype in SEMA3A gene knocked out mice (13). Liu et al. found that fasting serum SEMA3A was positively correlated with bone formation marker—osteocalcin (OCN), while no significant association was found of fasting serum SEMA3A with BMD/osteoporotic fractures (14).

Neuropilin-1 (NRP1) is a type I transmembrane protein which possesses high affinity to SEMA3A, and it plays an indispensable role in the chemorepulsion mediated by SMEA3A (15). PlexinA, one class of the plexin proteins, is the major receptors for SEMA3A. SEMA3A binds to NRP1 to assemble and to active NRP-plexinA holoreceptor complexes (16). It was found that the NRP1-mutant mice exhibited a same osteoporotic phenotype as the SEMA3A-deficient mice, suggesting that the SEMA3A function through NRP1 (13). SEMA3A enhanced the competitive binding of NRP1 (against TREM2) with plexinA1 (PLXNA1), thereby blocking the PLXNA1-TREM2-DAP12 complex formation to inhibit osteoclast differentiation (13). Otherwise, SEMA3A can also bind to the NRP1-PLXNA1 receptor complex to promote osteoblast formation (17). In vitro, plexinA2 (PLXNA2) were found to influence the expression of RUNX2 and SP7—two major osteogenic transcription factors, further regulating osteogenic differentiation and mineralization (18). PlexinA3 (PLXNA3) may participate in the process of SEMA3A regulating bone innervation (19).

Taking all evidence above, the SEMA3A signaling pathway may be involved in the OP pathology, and relevant genes may be associated with the disease susceptibility. However, epidemiological and experimental evidence is still scarce. Herein, we conducted a two-stage study in 1259 Chinese Han elderly to explore OP-related genetic variants in SEMA3A signaling pathway genes, including SEMA3A, NRP1, PLXNA1, PLXNA2 and PLXNA3, using next generation sequencing (NGS) and improved multiple linkage detection reaction (iMLDR) technique. And in vitro experiments were also conducted to explore potential functions and biological mechanisms of revealed OP-related genetic loci.

2 Method and materials

2.1 Study design and participants

The work of this study consisted of four parts. First, we conducted NGS of 5 candidate genes (SEMA3A, NRP1, PLXNA1, PLXNA2, and PLXNA3) using Illumina Hiseq 2500 sequencing platform based on microarray capture sequencing method to discover OP-related genetic variants in 400 Chinese postmenopausal women (211 OP cases). Second, the discovered OP-related single nucleotide polymorphisms (SNPs) were further genotyped using iMLDR technique in another 859 Chinese participants (417 women, 270 OP cases) in the validation stage. Third, we conducted a combination dataset analysis of genotyping data in all participants in the discovery and validation stages to improve the test performance. After the three steps for OP-related genetic variant identification, bioinformatic tools were adopted for function annotations to the identified variants. Fourth, in vitro experiments were conducted to clarify the roles of disease-related and functional genetic variants in bone mass remodeling. Considering that postmenopausal women is a high-risk group for OP, we included postmenopausal women exclusively in the discovery stage to identify as many OP-related genetic variants as possible. Detailed illustration to the study design is shown in Figure 1.

FIGURE 1

Figure 1 Framework of the study design.

Participants were recruited from two communities in Wuhan and physical examination and rehabilitation center of Wuhan Union hospital during 2017–2018. The trained staff conducted an interview by phone to determine whether these participants were eligible for this study. We included Chinese Han participants aged ≥60 years old. We excluded participants: (i) with previous hysterectomy, ovarian and adnexal resection; (ii) with other prevalent endocrine diseases that may potentially influenced bone metabolism (i.e., hyperthyroidism and hypothyroidism); (iii) with hereditary bone disease, including osteogenesis imperfecta, osteochondrosis and multiple myeloma; (iv) who ever took hormonal medicines for a long time (>3 months) or in the past 6 months; and (v) with unavailable data on candidate genes and BMD. Finally, a total of 1259 participants were involved in this study. OP diagnosis is defined as the T value ≤-2.5 (BMD values is ≤2.5 standard deviation below the BMD values in normal adult with same sex and race) based on the World Health Organization OP diagnostic criteria (2). Thus, those passed the inclusion and exclusion criteria and met OP diagnostic criteria were assigned to the OP group, and the rest were assigned to the non-OP group.

This study was approved by the Ethics Committee of Tongji Medical College of Huazhong University of Science and Technology. Written informed consent was obtained from each participant before enrollment.

2.2 Data collection and BMD measurement

Questionnaire survey was conducted to collect information on general characteristics (name, gender and age), medical history (fracture, uterus and ovaries surgery and thyroid-related disease), medication history (hormones), lifestyle (smoking and drinking), and menopause status for women only (menarche age and menopause age) from each participant. Body mass index (BMI) was obtained from physical examination, and it was calculated as weight (kg) divided by the square of height (m). BMD at the skeletal sites of the LS and FN was measured by Dual-energy X-ray absorptiometry (Lunar Prodigy, GE, USA).

2.3 Genotyping and function annotation

Before the physical examination, 3 mL venous blood was drawn from each participant for genotyping. The Relax Gene Blood DNA System Kit (Tiangen, Beijing, China) was adopted to extract genomic DNA. The NGS was used to screen for genetic variants that were potentially associated with OP in the discovery stage. IMLDR technique was used for genotyping the discovered OP-related SNPs in the validation stage. During this process, we set negative controls for each plate, and 5% of double-blind samples were randomly selected for repeat typing verification. The consistency of typing results was 100%, indicating the reliability of iMLDR technique.

Based on minor allele frequency (MAF), genetic variants were classified into common variants (MAF≥0.01) and rare variants (MAF<0.01). The Sorting Intolerant From Tolerant (SIFT) (http://sift.jcvi.org/www/SIFT_chr_coords_submit.html), Polyphen2 (http://genetics.bwh.harvard.edu/pph/pph_help_text.html), and Mutation Taster software (http://sites.google.com/site/jpopgen/dbNSFP) were adopted to predict deleteriousness of mutations.

2.4 In vitro experiments

2.4.1 Plasmid construction and lentivirus packaging

Sequence of NRP1 gene was obtained from the NCBI website. The cDNA sequence of wild-type NRP1 (allele T at rs180868035, I140) and mutant NRP1 (allele G at rs180868035, L140) were inserted into the vector plasmid, named pLVX-NRP1-Puro and pLVX-NRP1-Mut-Puro, respectively. The correctness of the inserted target fragment was verified by sequencing. Both plasmid construction and validation were done by Wuhan Vapol Biological company. The pLVX-NRP1-Puro and pLVX-NRP1-Mut-Puro were introduced into 293 T cells by lentiviral packaging to produce high titer lentivirus containing the target fragment.

2.4.2 Construction of MC3T3-E1 and RAW 264.7 stable-transfer cell lines

MC3T3-E1 and RAW 264.7 cell lines were used as models for osteoblast and osteoclast, respectively. Both cell lines were purchased from Wuhan Procell Life Science & Technology Co., Ltd. Construction of each stable-transfer cell line included 5 steps: (1) cell recovery and culture: α-MEM+10%FBS+1% (Penicillin-Streptomycin Solution) and DMEM+10%FBS+1% (Penicillin-Streptomycin Solution) were suitable culture medium for MC3T3-E1 and RAW 264.7 cell lines, separately; (2) cell passaging; (3) determine the optimal screening concentration of puromycin; (4) lentiviral infection; and (5) cell screening.

2.4.3 Measurement of osteoblast differentiation ability

Alkaline phosphatase (ALP) staining test was used to measure the activity of ALP, which was reckoned as the hallmark enzyme of mature osteoblast. The HiScript II Q Select RT SuperMix for qPCR (+ gDNA wiper) kit was used to reverse transcribe RNA into DNA, and real-time fluorescent quantitative PCR was performed to detect the mRNA expression levels of osteoblast marker genes, including ALP, RUNX2, SP7 and osteocalcin (OCN). In the process, the cDNA synthesized by reverse transcription, as a template, was amplified on an ABI QuantStudio 6 real-rime fluorescence quantitative PCR instrument using SYBR Green Master Mix reagent. The GAPDH gene was used as an internal reference for the relative quantification of osteoblast marker genes mRNA levels, and the relative quantification was calculated using the 2^-ΔΔCt method.

2.4.4 Measurement of osteoclast differentiation ability

Tartrate-resistant acid phosphatase (TRAP) staining test was used to measure the activity of TRAP, which was reckoned as the hallmark enzyme of mature osteoclast. The HiScript II Q Select RT SuperMix for qPCR (+ gDNA wiper) kit was used to reverse transcribe RNA into DNA, and real-time fluorescent quantitative PCR was performed to detect the mRNA expression levels of osteoclast marker genes, including TRAP, nuclear factor of activated T cells c1 (NFATc1) and cathepsin K (CTSK). Other processes were consistent with above.

2.4.5 Cell proliferation assay

The proliferation ability of osteoblast and osteoclast was detected by CCK8 cell proliferation and cytotoxicity assay kit. MC3T3-E1 and RAW 264.7 cells in good logarithmic growth phase were inoculated into 96-well plates at a concentration of 5×10³ and 2×10³ cells, respectively, per well, and three replicate wells were set up for each group. Then they were incubated at 37 °C in a 5% CO₂ incubator for 48 hours. And 10 μL CCK8 solution was added to each well and incubated for 4 hours at 37 °C, then the optical density (OD) value at 450 nm of each well was measured by enzyme-labeled instrument.

2.4.6 Cell apoptosis assay

The AnnexinV-APC/7-AAD apoptosis detection kit was adopted for apoptosis assay. The kit uses fluorescein APC-labeled Annexin V as a probe to reflect early apoptosis by binding phosphatidylserine exposed on the outside of the cell to the cytosol of early apoptosis cells. The 7-AAD, provided in this kit, is a nucleic acid dye that cannot penetrate the intact cell membrane of normal or early-stage apoptotic cells, but it can penetrate and bind to DNA within late-stage apoptotic or necrotic cells. Thus, it can be used to distinguish surviving early-stage cells from necrotic or late-stage apoptotic cells. When Annexin V-APC was used in combination with 7-AAD, 7-AAD was excluded from live cells and early apoptotic cells, while late apoptotic cells and necrotic cells were stained, presenting double positive status.

2.5 Statistical analyses

Continuous variables were presented as mean ± standard deviation (SD) or median (interquartile range, IQR) where appropriate, and Student’s t test or Mann-Whitney U test was used to exam between-group difference. Categorical variables were presented as number (percentage), and Pearson χ² test was adopted to exam difference between two groups. The χ² Goodness-of -Fit test was applied to assess whether the genotype distribution of genetic variants in the control group conformed to the Hardy-Weinberg Equilibrium (HWE). The associations between common variants and BMD/OP risk were analyzed using four kinds of multiple linear regression models or multivariate-adjusted logistic regression models, including co-dominant, additive, dominant and recessive models. For rare variants, considering the small sample size of control group, the East Asian population in the GnomAD database was also used as the control. Fisher’s exact test was adopted to compare the allele frequency and genotype distribution of rare variants between OP and non-OP groups, and linear regression was used to analyze the association of rare variants with BMD. The discovery stage was designed to investigate potentially OP-related genetic variant therefore the results of this stage were not corrected for false discovery rate (FDR). The other results were corrected for FDR to control the false positive rate caused by multiple comparisons. Statistical analyses were performed using IBM SPSS v22.0 and R v4.0.3 (R Core Team, Vienna, Austria), and picture production was achieved by GraphPad Prism v5.01. All statistical analyses were two-sided, and a P value < 0.05 indicated statistical significance.

3 Results

3.1 General characteristics

General characteristics of all participants were presented in Table 1. In the discovery stage involving 400 postmenopausal women (211 OP cases), the mean (SD) age of OP cases was 67.67 (5.87) years, comparatively larger than the non-OP individuals aged 66.26 years (4.63) (P=0.008). The mean (SD) BMI of OP cases was 23.22 (3.14) kg/m², significantly lower than the counterpart values (25.24 kg/m² (SD: 3.05)) of non-OP group (P<0.001). The proportions of smoker and drinker in both OP and non-OP groups were no more than 3.3% and 2.4%, and both two proportions were not significantly different between two groups. The mean (SD) BMD measurement and T value at the skeletal sites of LS were 0.80 (0.08) g/cm² and -3.12 (0.63), and were 0.66 (0.08) g/cm² and -2.66 (0.67) of FN in OP group, significantly lower than the counterpart values, sequentially were 1.08 (0.15) g/cm², -0.80 (1.09), 0.83 (0.09) g/cm² and -1.25 (0.77), in non-OP group (all P<0.001).

TABLE 1

Table 1 General characteristics of involved participants.

In the validation stage involving 859 older adults (270 OP cases), the proportion of women was 60.7% in the OP group, significantly higher than that in the non-OP group (43.0%) (P<0.001). The mean (SD) age of OP cases was 68.36 years (6.02), higher than the counterpart values (67.36 years (6.29)) of the non-OP group (P=0.029). The OP individuals had lower BMI (mean of 23.15 kg/m², SD of 3.15) than the non-OP individuals (mean of 24.65 kg/m², SD of 3.14) (P<0.001). Nearly similar proportion of OP (16.2%) and non-OP (17.5%) individuals smoke, and more reported drinking in the latter (18.5%) than the former (13.9%) without significance though. The mean (SD) BMD and T values at the skeletal sites of LS were 0.89 g/cm² (0.13) and -2.51 (1.03), and were 0.70 g/cm² (0.16) and -2.62 (0.69) of FN in OP group, significantly lower than the counterpart values (1.13 (0.18) g/cm² and -0.52 (1.43) at LS1-4, and 0.85 g/cm² (0.11) and -1.35 (0.81) at FN) of non-OP group (all P<0.001).

3.2 Results in the discovery stage

A total of 794 genetic variants were genotyped in the 5 SEMA3A signaling pathway genes. After excluding genetic variants (i) with low quality and possible repetitive sequences, (ii) with average sequencing depth ≤30×, (iii) with detection rate ≤90%, and (iv) not met HWE equilibrium, 586 variants were retained for subsequent analysis. Among them, 119 variants were new, the rs numbers of whom were unavailable based on databases of 1000G, ExAC, ESP6500 and GnomAD. Among the newly found variants, 12 were nonsynonymous mutations locating in the SEMA3A (n=2), NRP1 (n=1), PLXNA1 (n=3), PLXNA2 (n=3) and PLXNA3 (n=3) genes, respectively. It was predicted that two of them were deleterious to protein function by at least two of three software i.e., the SIFT, POLYphen2 and MutationTaster (Table S1).

Results of associations between the sequenced common variants and OP risk are presented in Table S2. After adjusting for age, BMI and menopause age, 6 common variants were associated with OP risk, including NRP1 rs2070296, PLXNA1 rs4679323 and rs73861745, and PLXNA2 rs1664227, rs2274446 and rs3748735. Compared to the wild genotype carriers (i.e., the rs4679323 CC, rs1664227 GG, rs2274446 CC and rs3748735 CC), carriers of rs4679323 CA (OR=0.53, 95% CI: 0.31, 0.90), rs2274446 CT (OR=0.55, 95% CI: 0.34, 0.88) and rs3748735 CT (OR=0.55, 95% CI: 0.32, 0.93) genotypes were less likely to have prevalent OP, while the rs1664227 CC genotype carriers were more likely to have prevalent OP (OR=2.04, 95% CI: 1.07, 3.86). Results of dominant, recessive and additive models were also presented (Table S2).

As shown in Table S3, the MAF distribution of four rare variants, i.e., the rs180868035 and rs767142032 in NRP1, and rs369477952 and rs146550621 in PLXNA1 were found to be significantly different between OP group and the East Asian population in the GnomAD database, with P values of 0.019, 0.022, 0.005 and 0.022 sequentially (Table S3).

3.3 Results in the validation stage

Above-mentioned 10 OP-related genetic variants (6 common and 4 rare ones) were further genotyped in another 859 subjects in the validation stage. Characteristics of the 10 variants are shown in Table S4. Compared to the NRP1 rs2070296 CC/CT genotype, the TT genotype was associated with an increased OP risk (OR=1.46, 95%CI: 1.01, 2.11), while the association did not achieve an FDR-corrected significance (Table 2). The PLXNA2 rs2274446 CT/TT genotype carriers had a 37% reduction in OP risk (OR=0.63, 95%CI: 0.46, 0.87, P_FDR=0.027) compared to the CC genotype carriers; and OP risk decreased by about 32% on average in response to a T to C mutant in the PLXNA2 rs2274446 (OR=0.68, 95%CI: 0.52, 0.89, P_FDR=0.029).

TABLE 2

Table 2 Multivariate logistic regression results of associations between common genetic variants and OP risk in the validation stage.

Results of associations of the 6 common variants and BMD measurements at skeletal sites of LS and FN are presented in Table 3. We found the PLXNA2 rs2274446 CT/TT genotype was associated with increased FN BMD (β=0.030 g/cm², 95%CI: 0.010, 0.050, P_FDR=0.004) compared to the CC genotype, and the FN BMD increased by 0.030 g/cm² (95% CI: 0.010, 0.040, P_FDR=0.006) on average in response to a T allele increase.

TABLE 3

Table 3 Multivariate linear regression results of associations between common genetic variants and BMDs in the validation stage.

As to the rare variant-BMD association, it was found that the BMD decreased in response to per rs180868035 G allele mutation at skeletal sites of LS (β=-0.282 g/cm², 95%CI: -0.507, -0.057) and FN (β=-0.208 g/cm², 95%CI: -0.375, -0.035), respectively.

3.4 Results in combination dataset analysis

In general, the significant associations found in the combination dataset analysis were consistent with those found in the validation stage. Notably, we also found the NRP1 rs2070296 TT genotype was associated with higher OP risk (OR=1.63, 95%CI: 1.21, 2.19, P_FDR=0.007) and lower FN BMD (β=-0.030 g/cm², 95% CI: -0.050, -0.010, P_FDR=0.004) compared to the CC/CT genotype (Tables 4–6).

TABLE 4

Table 4 Multivariate logistic regression results of associations between common genetic variants and OP risk in the combination dataset analysis.

TABLE 5

Table 5 Multivariate linear regression results of associations between common genetic variants and BMDs in the combination dataset analysis.

TABLE 6

Table 6 Multivariate linear regression results of associations between rare genetic variants and BMDs in the combination data analysis.

3.5 Effect of NRP1 rs180868035 mutation on osteoblast and osteoclast differentiation

Considering our findings from above work and the bioinformatics analysis of function annotations, the NRP1 rs180868035 was preferentially selected for function research via in vitro experiments. The NRP1 rs180868035 is a nonsynonymous mutation located in the exonic region. It was predicted by the SIFT and Mutation Taster software that the NRP1 rs180868035 T>G mutation can cause the coded amino acid (p.I140L) change (Table S4).

The ALP staining was positive in osteoblasts transfected with an empty vector, and vectors with either the wild-type (allele T at rs180868035, I140) or mutant (allele G at rs180868035, L140) NRP1 gene sequence. Compared to the osteoblasts transfected with empty vector, wild-type transfected osteoblasts apparently enhanced the depth of ALP staining (Figure S1), as well as mRNA expression of the osteoblast marker genes, i.e., ALP, RUNX2, SP7 and OCN (Figure S2). Moreover, compared to the wild-type transfected cells, the depth of ALP staining was apparently reduced (Figure S1), and osteoblast marker genes all were lower in the mutant transfected cells (Figure S2).

The TRAP staining was positive in osteoclasts transfected with an empty vector, and vectors with either the wild-type or mutant NRP1 gene sequence. Compared to the osteoclasts transfected with empty vector (Figure S3), the depth of TRAP staining and osteoclast marker genes expression level, i.e., TRAP, NFATc1 and CTSK, all were reduced in the wild-type transfected osteoclasts (Figure S4). Besides, compared to the wild-type transfected cells, mutant transfection significantly enhanced the depth of TRAP staining (Figure S3), and mRNA expression level of osteoclast marker genes (Figure S4).

3.6 Effect of NRP1 rs180868035 mutation on osteoblast and osteoclast proliferation

The average proliferation rates of osteoblasts transfected an empty vector, and vectors with either the wild-type and mutant NRP1 were 97.3%, 136.8% and 103.5%, respectively. As shown in Figure 2A, compared to those transfected with the empty vector, the proliferation rate was significantly promoted in osteoblasts transfected with the wild-type NRP1; and compared to those transfected with the wild-type, the proliferation rate was significantly inhibited in osteoblasts transfected with the mutant NRP1.

FIGURE 2

Figure 2 Proliferation rate comparisons among osteoblasts (A) and osteoclasts (B) transfected with empty and wild-type and mutant of NRP1 rs180868035 vectors. ** indicated P < 0.001.

The average proliferation rates of osteoclasts transfected an empty vector, and vectors with either the wild-type and mutant NRP1 were 92.8%, 61.5% and 86.1%, respectively. As shown in Figure 2B, compared to those transfected with the empty vector, the proliferation rate was significantly inhibited in osteoclasts transfected with the wild-type NRP1; and compared to those transfected with the wild-type, the proliferation rate was significantly promoted in osteoclasts transfected with the mutant NRP1.

3.7 Effect of NRP1 rs180868035 mutation on osteoblast and osteoclast apoptosis

The average apoptosis rates of osteoblasts transfected with an empty vector, and vectors with either the wild-type and mutant NRP1 were 6.2%, 3.8% and 5.8%, respectively. Compared to those transfected with the empty vector, the apoptosis rate was significantly reduced in osteoblasts transfected with the wild-type NRP1; and compared to those transfected with the wild-type, the apoptosis rate was significantly increased in osteoblasts transfected with the mutant NRP1 (Figure 3).

FIGURE 3

Figure 3 Cellular apoptosis comparisons among osteoblasts transfected with empty, and wild-type and mutant of NRP1 rs180868035 vectors. (A) flow cytometry. Representative APC-A-area versus 7-AAD-A-area density plots of different treatment groups with percentages of cells at different stage of apoptosis shown in the relevant quadrants. Q1-UR quadrant is positive for both APC-A and 7-AAD-A and represents cells in late apoptosis; Q1-LR is APC-A positive/7-AAD-A negative and represents cells in early apoptosis; Q1-LL is negative for both APC-A and 7-AAD-A and contains viable cells; Q1-UL is 7-AAD-A positive/APC-A negative and contains cells undergoing necrosis. (B) histogram of osteoblast apoptosis. ** indicated P < 0.001.

The average apoptosis rates of osteoclasts transfected with an empty vector, and vectors with either the wild-type and mutant NRP1 were 8.1%, 27.0% and 10.2%, separately. Compared to those transfected with the empty vector, the apoptosis rate was significantly elevated in osteoclasts transfected with the wild-type NRP1; and compared to those transfected with the wild-type, the apoptosis rate was significantly decreased in osteoclasts transfected with the mutant NRP1 (Figure 4).

FIGURE 4

Figure 4 Cellular apoptosis comparisons among osteoclasts transfected with empty, and wild-type and mutant of NRP1 rs180868035 vectors. (A) flow cytometry. Representative APC-A-area versus 7-AAD-A-area density plots of different treatment groups with percentages of cells at different stage of apoptosis shown in the relevant quadrants. Q1-UR quadrant is positive for both APC-A and 7-AAD-A and represents cells in late apoptosis; Q1-LR is APC-A positive/7-AAD-A negative and represents cells in early apoptosis; Q1-LL is negative for both APC-A and 7-AAD-A and contains viable cells; Q1-UL is 7-AAD-A positive/APC-A negative and contains cells undergoing necrosis. (B) histogram of osteoclast apoptosis. ** indicated P < 0.001.

4 Discussion

In the present study, we conducted a two-stage designed epidemiological research in combination with bioinformatics analysis to identify functional genetic variants associated with OP risk/BMD in 5 SEMA3A pathway genes (SEMA3A, NRP1, PLXNA1, PLXNA2 and PLXNA3). In vitro experiments were further performed to investigate roles of identified variants in bone mass remodeling. The NRP1 rs2070296 and rs180868035, and PLXNA2 rs2274446 were revealed to be associated with OP risk/BMD. Based on function annotations to those variants, in vitro experiments were conducted and showed that the NRP1 rs180868035 mutation affect differentiation, proliferation and apoptosis of osteoblasts and osteoclasts.

Our findings showed that the NRP1 rs2070296 TT genotype was associated with increased OP risk and decreased FN BMD, and LS and FN BMDs decreased in response to the NRP1 rs180868035 G allele mutation. The rs2070296 C>T (V179V) is a synonymous mutation at physical position of chr10:33552695, located in exon 4 of the NRP1 gene. According to the rSNPBase database (http://rsnp3.psych.ac.cn/index.do), the rs2070296 is either a regulatory SNP or an expression quantitative trait locus, indicating it may be involved in the regulation of NRP1 gene expression. We found little evidence on the association of the rs2070296 variant with specific outcomes other than one literature observed a cumulative effect of SNP rs2070296 on visual acuity underwent ranibizumab treatment (20). The rs180868035 T>G (I140L) is a nonsynonymous mutation at physical position of chr10:33559615 located in exon 3 of the NRP1 gene, corresponding to the change of Ile amino acid to Leu amino acid (p.I140L). It was predicted by the SIFT algorithm to be deleterious to protein function (21). No study on the association between genetic variant of NRP1 rs180868035 and bone metabolism was reported before ours.

We further performed in vitro experiments to explore the influence of NRP1 mutation on bone metabolism. Our findings revealed that overexpression of a rs180868035 T>G mutant decreased the expression levels of osteoblast marker genes, i.e., ALP (22), RUNX2 (23), SP7 (24) and OCN (25), inhibited osteoblast proliferation and promoted osteoblast apoptosis. In the meanwhile, the overexpression of the mutant was found to enhance the expression of osteoclast marker genes, i.e., TRAP (26), NFATc1 and CTSK (27), promote osteoclast proliferation and inhibit osteoclast apoptosis. Findings from in vitro experiments corroborated the association of NRP1 rs180868035 with bone phenotype found in our epidemiological study. The N-terminal of SEMA3A contains a seven-bladed β-propeller structure, which was reckoned as the signature ‘Sema’ domain, determining its unique ability to bind to NRP1, and the NRP1 rs180868035 genetic variant is located in the functional domain in extracellular region of NRP1 that binds to the SEMA3A (28, 29). The Npn1-mutant mice (NRP1^Sema- mice), lacking the Sema-binding site, exhibited a same osteoporotic phenotype as SEMA3A-deficient mice, indicating that the NRP1 functions through binding to SEMA3A (13). In addition, the binding of SEMA3A to NRP1 could inhibit osteoclast differentiation induced by RANKL signaling and enhance osteoblast differentiation through classic WNT/β catenin pathway (9, 13). We therefore speculated that the NRP1 mutation caused amino acid sequence alteration that resulted in the inability of NRP1 to bind to SEMA3A, further disrupting the differentiation, proliferation and apoptosis of osteoblasts and osteoclasts. The NRP1 served as a co-receptor and interacted with the discoidin domain receptor 2 (DDR2)—a receptor tyrosine kinase possessing a role of promoting bone-depositing cells and osteoblasts differentiation. Animal experiments suggested that NRP1 enhanced the bone formation induced by DDR2, and knockdown of NRP1 or ectopic expression of NRP1 can influence osteoclastogenesis in cells with different DDR2 expression levels (30). Thus, NRP1 mutation might inhibit the interaction between NRP1 and DDR2. It is an intriguing idea for the future studies. Besides, NRP1 expression was tightly regulated by growth factors, transcription factors, and injury, e.g., tumor necrosis factor-alpha, Prox-1, and ischemia (31). But evidence on how these modulators regulate NRP1 in the process of OP development is still scarce. More studies are needed to investigate potential mechanisms.

We found genetic variant of PLXNA2 rs2274446 (C > T) was associated with reduced OP risk and increased FN BMD. PLXNA2 is one of transmembrane receptors and can mediate SEMA3A signals through binding to NRP1 as a co-receptor for SMEA3A (32). Oh et al. found that the PLXNA2 expression could be upregulated by the bone morphogenic protein 2, subsequently activating the expression of two major osteogenic transcription factors, Runx2 and SP7, and reinforcing osteogenic mineralization and differentiation, indicating that the PLXNA2 possesses potential ability to affect bone remodeling and is a noteworthy target for OP treatment (18). PLXNA2 is a receptor of SEMA6A in osteoclasts. The SEMA6A-PLXNA2 axis participated in osteoclastogenesis induced by RANKL, and this function is dependent on the activation of NFATc1 induced by PLCγ (33). It offers potential therapeutic targets in the intervention of OP. In a study including 560 postmenopausal women of Korean ethnicity, ten SNPs in the PLXNA2 gene were selected for genotyping and the rs3748735 variant was found to be significantly associated with FN BMD (32). Genetic variant of PLXNA2 rs3748735 was found to be associated with FN BMD at the discovery stage in our study, while the significant association did not maintain at the further validation stages. This may be partly due to different sample sizes and study populations.

As far as we known, our study is one of few studies to explore the association of genetic variants in the SEMA3A signaling pathway genes with OP risk. Our work consisted of population-based epidemiological study and in vitro experiments, concluding genetic variant of NRP1 rs180868035 may affect BMD by increasing bone resorption and reducing bone formation. Results from these two parts corroborate and complement each other. Thus, our study can provide stable and valid evidence for OP susceptibility loci identification. Meanwhile, several limitations should be mentioned. First, though we found genetic variant of NRP1 rs180868035 could affect osteoblast and osteoclast proliferation and apoptosis, the specific biological mechanism through which NRP1 rs180868035 influence bone metabolism needs to be elucidated in the future. Second, the extensibility of our findings will be limited because we only included Chinese Han residents living in Wuhan, China. Other population-based studies are needed to confirm our findings. Third, in the in vitro experiment part, the early stage of osteoblasts and osteoclasts differentiation was missed to be measured in this study. At last, several potential confounders, e.g., physical activity (34), serum calcium and vitamin D levels (35), were not adjusted in the present study due to insufficient data, which would bias the conclusions.

5 Conclusion

This study firstly and systematically analyzed the associations between genetic variants in the SEMA3A signaling pathway genes and OP risk/BMD in a Chinese population. We concluded genetic variants of NRP1 rs2070296 and rs180868035, and PLXNA2 rs2274446 were associated with OP risk and BMDs. The NRP1 rs180868035 influences bone metabolism by regulating osteoblast and osteoclast differentiation, proliferation and apoptosis. More researches are needed to prove our findings.

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: www.ncbi.nlm.nih.gov/bioproject/, PRJNA884817.

Ethics statement

The studies involving human participants were reviewed and approved by The Ethics Committee of Tongji Medical College of Huazhong University of Science and Technology. The patients/participants provided their written informed consent to participate in this study.

Author contributions

Conceptualization: QW, QH, H-LZ and M-HW. Study conduct: H-LZ and M-HW. Data collection and preparation: D-SD, R-YZ, J-LZ, T-TY, QL and T-TZ. Data analysis: M-HW and H-LZ. Writing - Original draft preparation: H-LZ and M-HW. Wring – Review and Editing: QW and QH. All authors read, reviewed and approved the final manuscript. QW and QH had primary responsibility for final content. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Natural Science Foundation of China (Grant nos.82273711 and 72061137006) and Wuhan Municipal Health Commission (Grant nos. WY22B06 and WG20C09).

Acknowledgments

We would like to thank the staff of the two communities in Wuhan and physical examination and rehabilitation center of Wuhan Union hospital for their administrative and technology assistance.

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.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fendo.2022.1014431/full#supplementary-material

Abbreviations

ALP, alkaline phosphatase; BMD, bone mineral density; BMI, body mass index; CTSK, cathepsin K; DXA, dual-energy X-ray absorptiometry; FDR, false discovery rate; FN, femoral neck; HWE, Hardy-Weinberg equilibrium; iMLDR, improved multiple linkage detection reaction; IQR, interquartile range; LS, lumbar spine; MAF, minor allele frequency; NFATc1, nuclear factor of acid phosphatase; NGS, next generation sequencing; NRP1, neuropilin 1; OCN, osteocalcin; OD, optical density; OP, osteoporosis; PLXNA1, plexinA1; PLXNA2, plexinA2; PLXNA3, plexinA3; SD, standard deviation; SEMA3A, semaphoring 3A; SIFT, sorting intolerant from tolerant; SNP, single nucleotide polymorphisms; TRAP, tartrate-resistant acid phosphatase.

References

1. Baccaro LF, Conde DM, Costa-Paiva L, Pinto-Neto AM. The epidemiology and management of postmenopausal osteoporosis: a viewpoint from Brazil. Clin Interv Aging (2015) 10:583–91. doi: 10.2147/CIA.S54614

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Camacho PM, Petak SM, Binkley N, Diab DL, Eldeiry LS, Farooki A, et al. American Association of clinical Endocrinologists/American college of endocrinology clinical practice guidelines for the diagnosis and treatment of postmenopausal osteoporosis-2020 update. Endocr Pract (2020) 26:1–46. doi: 10.4158/GL-2020-0524SUPPL

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Salari N, Ghasemi H, Mohammadi L, Behzadi MH, Rabieenia E, Shohaimi S, et al. The global prevalence of osteoporosis in the world: A comprehensive systematic review and meta-analysis. J Orthop Surg Res (2021) 16:609. doi: 10.1186/s13018-021-02772-0

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Kanis JA, Cooper C, Rizzoli R, Reginster JY, Scientific Advisory Board of the European Society for Clinical and Economic Aspects of Osteoporosis (ESCEO), Committees of Scientific Advisors and National Societies of the International Osteoporosis Foundation (IOF). European Guidance for the diagnosis and management of osteoporosis in postmenopausal women. Osteoporos Int (2019) 30:3–44. doi: 10.1007/s00198-018-4704-5

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Looker AC, Melton LJ 3rd, Harris TB, Borrud LG, Shepherd JA. Prevalence and trends in low femur bone density among older US adults: NHANES 2005-2006 compared with NHANES III. J Bone Miner Res (2010) 25:64–71. doi: 10.1359/jbmr.090706

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Zeng Q, Li N, Wang Q, Feng J, Sun D, Zhang Q, et al. The prevalence of osteoporosis in China, a nationwide, multicenter DXA survey. J Bone Miner Res (2019) 34:1789–97. doi: 10.1002/jbmr.3757

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Xu R. Semaphorin 3A: A new player in bone remodeling. Cell Adh Migr (2014) 8:5–10. doi: 10.4161/cam.27752

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Rienks M, Carai P, Bitsch N, Schellings M, Vanhaverbeke M, Verjans J, et al. Sema3A promotes the resolution of cardiac inflammation after myocardial infarction. Basic Res Cardiol (2017) 112:42. doi: 10.1007/s00395-017-0630-5

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Kim JM, Lin C, Stavre Z, Greenblatt MB, Shim JH. Osteoblast-osteoclast communication and bone homeostasis. Cells (2020) 9:2073. doi: 10.3390/cells9092073

CrossRef Full Text | Google Scholar

10. Liu F, Wang C, Huang H, Yang Y, Dai L, Han S, et al. SEMA3A-mediated crosstalk between prostate cancer cells and tumor-associated macrophages promotes androgen deprivation therapy resistance. Cell Mol Immunol (2021) 18:752–4. doi: 10.1038/s41423-021-00637-4

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Hayashi M, Nakashima T, Yoshimura N, Okamoto K, Tanaka S, Takayanagi H. Autoregulation of osteocyte Sema3A orchestrates estrogen action and counteracts bone aging. Cell Metab (2019) 29:627–37.e625. doi: 10.1016/j.cmet.2018.12.021

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Zhang N, Hua Y, Li Y, Pan J. Sema3A accelerates bone formation during distraction osteogenesis in mice. Connect Tissue Res (2021) 1-11:382–92. doi: 10.1080/03008207.2021.1974850

CrossRef Full Text | Google Scholar

13. Hayashi M, Nakashima T, Taniguchi M, Kodama T, Kumanogoh A, Takayanagi H. Osteoprotection by semaphorin 3A. Nature (2012) 485:69–74. doi: 10.1038/nature11000

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Liu DM, Lu N, Zhao L, Zhang MJ, Tao B, Xuan Y, et al. Serum Sema3A is in a weak positive association with bone formation marker osteocalcin but not related to bone mineral densities in postmenopausal women. J Clin Endocrinol Metab (2014) 99:E2504–2509. doi: 10.1210/jc.2014-1443

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Gu C, Rodriguez ER, Reimert DV, Shu T, Fritzsch B, Richards LJ, et al. Neuropilin-1 conveys semaphorin and VEGF signaling during neural and cardiovascular development. Dev Cell (2003) 5:45–57. doi: 10.1016/S1534-5807(03)00169-2

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Zhou Y, Gunput RA, Pasterkamp RJ. Semaphorin signaling: progress made and promises ahead. Trends Biochem Sci (2008) 33:161–70. doi: 10.1016/j.tibs.2008.01.006

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Worzfeld T, Offermanns S. Semaphorins and plexins as therapeutic targets. Nat Rev Drug Discov (2014) 13:603–21. doi: 10.1038/nrd4337

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Oh JE, Kim HJ, Kim WS, Lee ZH, Ryoo HM, Hwang SJ, et al. PlexinA2 mediates osteoblast differentiation via regulation of Runx2. J Bone Miner Res (2012) 27:552–62. doi: 10.1002/jbmr.1471

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Gomez C, Burt-Pichat B, Mallein-Gerin F, Merle B, Delmas PD, Skerry TM, et al. Expression of semaphorin-3A and its receptors in endochondral ossification: potential role in skeletal development and innervation. Dev Dyn (2005) 234:393–403. doi: 10.1002/dvdy.20512

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Lores-Motta L, van Asten F, Muether PS, Smailhodzic D, Groenewoud JM, Omar A, et al. A genetic variant in NRP1 is associated with worse response to ranibizumab treatment in neovascular age-related macular degeneration. Pharmacogenet Genomics (2016) 26:20–7. doi: 10.1097/FPC.0000000000000180

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Sim NL, Kumar P, Hu J, Henikoff S, Schneider G, Ng PC, et al. SIFT web server: predicting effects of amino acid substitutions on proteins. Nucleic Acids Res (2012) 40:W452–457. doi: 10.1093/nar/gks539

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Vimalraj S. Alkaline phosphatase: Structure, expression and its function in bone mineralization. Gene (2020) 754:144855. doi: 10.1016/j.gene.2020.144855

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Komori T. Signaling networks in RUNX2-dependent bone development. J Cell Biochem (2011) 112:750–5. doi: 10.1002/jcb.22994

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Liu Q, Li M, Wang S, Xiao Z, Xiong Y, Wang G, et al. Recent advances of osterix transcription factor in osteoblast differentiation and bone formation. Front Cell Dev Biol (2020) 8:601224. doi: 10.3389/fcell.2020.601224

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Zoch ML, Clemens TL, Riddle RC. New insights into the biology of osteocalcin. Bone (2016) 82:42–9. doi: 10.1016/j.bone.2015.05.046

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Song C, Yang X, Lei Y, Zhang Z, Smith W, Yan J, et al. Evaluation of efficacy on RANKL induced osteoclast from RAW264.7 cells. J Cell Physiol (2019) 234:11969–75. doi: 10.1002/jcp.27852

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Ono T, Nakashima T. Recent advances in osteoclast biology. Histochem Cell Biol (2018) 149:325–41. doi: 10.1007/s00418-018-1636-2

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Appleton BA, Wu P, Maloney J, Yin J, Liang WC, Stawicki S, et al. Structural studies of neuropilin/antibody complexes provide insights into semaphorin and VEGF binding. EMBO J (2007) 26:4902–12. doi: 10.1038/sj.emboj.7601906

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Usuki S, Yasutake Y, Tamura N, Tamura T, Tanji K, Saitoh T, et al. Nrp1 is activated by konjac ceramide binding-induced structural rigidification of the a1a2 domain. Cells (2020) 9:517. doi: 10.3390/cells9020517

CrossRef Full Text | Google Scholar

30. Zhang Y, Su J, Wu S, Teng Y, Yin Z, Guo Y, et al. DDR2 (discoidin domain receptor 2) suppresses osteoclastogenesis and is a potential therapeutic target in osteoporosis. Sci Signal (2015) 8:ra31. doi: 10.1126/scisignal.2005835

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Bielenberg DR, Pettaway CA, Takashima S, Klagsbrun M. Neuropilins in neoplasms: expression, regulation, and function. Exp Cell Res (2006) 312:584–93. doi: 10.1016/j.yexcr.2005.11.024

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Hwang JY, Lee JY, Park MH, Kim KS, Kim KK, Ryu HJ, et al. Association of PLXNA2 polymorphisms with vertebral fracture risk and bone mineral density in postmenopausal Korean population. Osteoporos Int (2006) 17:1592–601. doi: 10.1007/s00198-006-0126-x

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Zhuang J, Li X, Zhang Y, Shi R, Shi C, Yu D, et al. Sema6A-plexin-A2 axis stimulates RANKL-induced osteoclastogenesis through PLCgamma-mediated NFATc1 activation. Life Sci (2019) 222:29–35. doi: 10.1016/j.lfs.2019.01.060

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Howe TE, Shea B, Dawson LJ, Downie F, Murray A, Ross C, et al. Exercise for preventing and treating osteoporosis in postmenopausal women. Cochrane Database Syst Rev (2011), CD000333. doi: 10.1002/14651858.CD000333.pub2

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Liu M, Yao X, Zhu Z. Associations between serum calcium, 25(OH)D level and bone mineral density in older adults. J Orthop Surg Res (2019) 14:458. doi: 10.1186/s13018-019-1517-y

PubMed Abstract | CrossRef Full Text | Google Scholar