Overexpression of SNTG2, TRAF3IP2, and ITGA6 transcripts is associated with osteoporotic vertebral fracture in elderly women from community

Abstract Background Vertebral fractures (VFs) are the most common clinical manifestation of osteoporosis associated with high morbimortality. A personal/familiar history of fractures increases the risk of fractures. The purpose of this study is to identify possible molecular markers associated with osteoporotic VFs in elderly women from community. Methods Transcriptomic analysis using Affymetrix HTA2 microarray was performed using whole blood samples of 240 subjects from a population‐based survey (Sao Paulo Ageing & Health [SPAH] study). Only elderly women with osteoporosis diagnosis by densitometry were analyzed, and divided in two groups: VF: women with osteoporosis and VFs versus no vertebral fracture (NVF): women with osteoporosis and NVFs. They were matched for age, chronic disease, medication use, and bone mineral density (BMD). The logistic regression model adjusted for age was applied for transcriptome data analysis. SYBR green‐based quantitative polymerase chain reaction (qPCR) was used to validate the most significant expression changes obtained in the microarray experiment. Results Microarray analysis identified 142 differentially expressed genes (DEGs, p < .01), 57 upregulated and 85 downregulated, compared VF versus NVF groups. The DEG with the greatest expression difference was the Gamma2‐Syntrophin (SNTG2) (β = 31.88, p = .005). Validation by qPCR confirmed increased expression in VF group of Syntrophin (SNTG2, fold change = 2.79, p = .009), TRAF3 Interacting Protein2 (TRAF3IP2, fold change = 2.79, p = .020), and Integrin Subunit Alpha 6 (ITGA6, fold change = 2.86, p = .038). Conclusion Our data identified and validated the association of SNTG2 (608715), TRAF3IP2 (607043), and ITGA6 (147556) with osteoporotic VF in elderly women, independently of BMD. These results suggest that these transcripts have potential clinical significance and may help to explain the molecular mechanisms and biological functions of vertebral fracture.


| INTRODUCTION
Osteoporosis is a bone disorder characterized by low bone density and microarchitecture deterioration, which may result in fragility and fracture risk, occurring most commonly in postmenopausal women (Brunton et al., 2005).
The Framingham Heart Study showed that VF risk, volumetric bone mineral density (BMD), and geometry at the lumbar spine are complex genetic traits related to gene-environment interactions (Liu et al., 2011). Family history of a VF represents an independent risk factor for incident fractures, (Kanis et al., 2004) showing the significance of heredities in the pathogenesis of bone metabolism deterioration. A genome-wide association study (GWAS) identified 56 BMD loci and reveals 14 loci associated with the risk of all fractures (Estrada et al., 2012).
To date, there are no studies in the literature analyzing the processes associated with VFs on the transcriptomic level. Based on that, the aim of this work was to find differentially expressed genes (DEGs) and pathways associated with osteoporotic VFs.

| Subjects
This study was conducted from June 2015 to July 2016 and subjects were selected from a previous epidemiologic project, which included a population-based survey (São Paulo Ageing & Health [SPAH] study). Two hundred and forty individuals were evaluated and RNA samples were extracted. Of them, 92 men were excluded to avoid comparison with different sexes, 20 women were excluded due to low RNA concentration (<100 ng/dl), 34 women were excluded due to chip hybridization failure, 4 women were excluded in normalization from the data for bioinformatics, 50 women were excluded for BMD T-score >−2.5, and 4 women were excluded for presenting only grade 1 fracture. Therefore, 36 elderly women (≥65-year age) with osteoporosis diagnosis by densitometry were analyzed. They were divided into two age-matching groups: (a) vertebral fracture (VF, n = 24): osteoporosis with VF and (b) no vertebral fracture (NVF, n = 12): osteoporosis with NVF.
All subjects responded to a standardized questionnaire regarding lifestyle and health, including history of hip fracture, fragility and previous fractures during last year, physical activity, alcohol intake, falls in the last 12 months, smoking status, personal history of hypertension, diabetes mellitus cardiovascular events (myocardial infarction and cerebrovascular disease), dietary calcium intake, menopause age, and all were apparently healthy. No individual had malabsorption or chronic diarrhea or hepatic disease or had severe chronic diseases or cancer or using glucocorticoid. BMD and VFs were assessed by dual-energy X-ray absorptiometry (DXA). Race was defined based on the self-reported race of second-generation ancestors, an approach previously used for the Brazilian population. Individuals with four grandparents reported as White were classified as White, while those individuals with Black and White ancestors (mixed race) were classified as non-White. Individuals with Asian ancestors were also classified as non-White. When racial information regarding the grandparents was not available, an individual's race was determined by the race of his or her parents (Lopes et al., 2011).
The Local Ethics in Research Committee of the São Paulo University School of Medicine approved this study, and all participants gave written informed consent.

| Blood collection, RNA extraction, and microarray analysis
Twenty milliliter of peripheral whole blood samples was collected under fasting conditions (between 8:00 a.m. and 10:00 a.m.) and stored at −80°C until analysis. These samples were used to perform biochemical tests of bone metabolism (calcium, phosphorus, alkaline phosphatase, 25OHD, iPTH, CTX, P1NP) and RNA extraction. About 10 ml was collected in PAXgene Blood RNA tubes, and total RNA was extracted, isolated, and purified with RNeasy columns (Qiagen) and treated with RNase-Free DNase Set (Qiagen).
Bioanalyzer 2100 System (Agilent) to analyze ribosomal bands 28S and 18S was used to confirm purified RNA. Samples with a minimum of 100 ng of RNA underwent in vitro transcription with the Ambion WT Expression kit (Life Technologies) in Affymetrix GeneChip ® WT Terminal Labeling and Control kit. The RNA quality was evaluated with the NanoDrop spectrophotometer (ND-1000, NanoDrop, Thermo Scientific). Labeled samples were hybridized to the Affymetrix GeneChip ® Human Transcriptome Array HTA 2.0. Array hybridization, washing, staining, and scanning were done as per subsequent Affymetrix instructions. 67,528 transcribed genes covering 44,699 coding genes and 22,829 noncoding regions were analyzed. After quality control and assessment of the expression signals, CEL files have been subjected to further analysis using Affymetrix GCOS ® .

| Evaluation of the serum biochemical bone parameters (calcium, phosphorus, and alkaline phosphatase)
The concentrations of serum calcium, phosphorus, and alkaline phosphatase were measured using standard automated laboratory methods. Serum concentrations of 25-hydroxyvitamin D (25OHD) and iPTH were determined using a radioimmunoassay technique (DiaSorin) and an immunoradiometric assay (ELISA-PTH; CIS Bio International), respectively. C-terminal telopeptide of type I collagen (CTX) and N-terminal propeptide of type I procollagen (P1NP) serum concentrations were measured by automated Electrochemiluminescence (E411, Roche Diagnostics). The serum CTX limit of detection was 10 ng/ml and intra-assay and inter-assay coefficients variation (CVs) were 2.5% and 3.4%, respectively. Serum P1NP intra-assay and inter-assay CVs were 2.2 and 1.8%, respectively, with a limit of detection of 5% µg/L.

| BMD and semi-quantitative measurements for VF
Bone mineral density was determined by DXA, using a Hologic QDR 4500A densitometry equipment (Hologic Inc., Discovery model) at the lumbar spine, femoral neck, and total hip. The precision errors for BMD measurements were measured based on the standard protocols from The International Society for Clinical Densitometry (ISCD) (Bennett et al., 2006). The prevalent VF was assessed from T4 to L4, a Vertebral Fracture Assessment (VFA) scan of the thoracolumbar spine performed using the same DXA device, with subjects in the supine position. All DXA measurements were performed by the same experienced technologist (LT).
On that same week, standard lateral thoracic and lumbar spine radiographs were taken using a 40-inch tube-to-film distance centered at T7 and L2. All VFA and radiographic images were independently evaluated by two experienced rheumatologists (GHM and LHJN). For the identification of VFs, in both methods, the readers assessed each T4-L4 image to decide whether a fracture was present. A consensus was agreed between the readers for any difference in interpretation. Interobserver agreement was 96% and the Kappa coefficient was 0.83. Nonvisible vertebrae were excluded. Only adequately visualized vertebrae were analyzed for deformity using the Genant semi-quantitative (SQ) method (Genant, Wu, Vankuijk, & Nevitt, 1993). Each identified fractured vertebra was classified by grade based on the Genant SQ scale, where mild (grade 1) is a reduction of 20%-25% of anterior, middle, and/or posterior height relative to the adjacent vertebral bodies; moderate (grade 2) is a reduction of 26%-40% in any height; and severe (grade 3) is a reduction of 40% in any height. For this study, grade I fractures were excluded; therefore, subjects with grade II or grade III VFs were included in the VF group and subjects without any vertebral deformity were included in the NVF group.

| Microarray data preprocessing and analysis
The raw data were preprocessed in R with affy and multtest packages from Bioconductor (Parmigiani, Garrett, Irizarry, & Zeger, 2003). Quality control of individual arrays included various advanced quality metrics, diagnostic plots, and pseudo-images to certify that only excellent quality arrays were used prior to downstream statistical analysis. The robust multiarray average (RMA) algorithm was used for background correction and quantile normalization, and median polish summarization was made using the software expression console. The raw intensity of each probe was processed by R program (version v2.12.1; https://www.r-proje ct.org). Probes that passed the criteria were normalized by 75% median scaling normalization method. Normalized spot intensities were transformed to gene expression log 2 ratios between the control and fracture groups. The probes with log 2 ratio ≥1 or log 2 ratio ≤−1 and p < .01 were defined as differential genes for additional pathway enrichment analysis. In this section, we present the results of the adjustment of the logistic model taking into account the binary variables, osteoporosis and VF .
Statistical analysis of the data was performed in R by using a logistic regression model taking into account age adjustment. In this approach, downregulation reflects an association with no risk of VF, whereas upregulation reflects an association with VF presence. The probes with β ≥1 or β ≤−1 and p ≤ .01 were defined as DEGs and after their identification, additional pathway enrichment analysis was performed.

| Gene ontology and biological pathways enrichment analysis
In further steps of the analysis, DEGs (gene symbols with β ratios) were uploaded into the Ingenuity Pathway Analysis (IPA) software (Ingenuity Systems) to identify relevant differential expression pathways. Gene ontology (GO) analysis and pathway enrichment analysis are commonly used for functional studies of large-scale transcriptomic or genomic data (Hulsegge, Kommadath, & Smits, 2009). We performed these analyses using up-and downregulated DEGs at p ≤ .01 adjusted for multiple comparisons. DEGs were also enriched in the biological process (BP), cellular component (CC), and molecular function (MF) categories of GO database. FUMA GWA was used to search for overrepresentation of the hallmark gene sets, biological pathways. The Database Functional Mapping and Annotation of Genome-Wide Association Studies (FUMA GWA, http://fuma.ctglab.nl) provides a set of data-mining tools and assists in the interpretation of genome-scale datasets by expediting the transition from data collection to biological annotation. GO terms with p < .05 based on a hypergeometric test were defined to be statistically significant. Basically, upregulation reflects an association with risk of developing VFs and downregulation revels decreased risk of VF in women with osteoporosis.

| Interaction network analysis
The interaction network for DEGs was constructed using Cytoscape 3.6 software and stringApp implementing 1.4 String interaction database (Shannon et al., 2003). Interaction confidence cutoff was set as ≥0.4. Tissue expression confidence levels were obtained from stringApp. Age-related expression changes in muscle tissue were obtained from the GTEX study (Yang et al., 2016).

| Quantitative polymerase chain reaction assays of significantly differentially expressed mRNAs
cDNA was synthesized from 100 ng of RNA using cDNA Reverse Transcription kits (Applied Biosystems) and thermocycler PTC-100 (MJ Research) was used for cDNA synthesis. Quantitative PCR was performed using PowerUp ™ SYBR ™ Green Master Mix (Thermo Fisher Scientific) on the StepOne Plus quantitative polymerase chain reaction (qPCR) (Applied Biosystems) cycler according to the manufacturer's instructions. All samples were tested in triplicates.
PCR primers for tested genes were designed based on their DE microarray probe exomic location. Primers were designed using Primer3, OligoCalc, OligoAnalyzer, and Primer-BLAST (Owczarzy, 2018) online tools. Primer sequences are in Table S1. To validate the changes in gene expression, we used the Pfaffl method for relative expression quantification (Pfaffl, 2001). The geometric mean of the two most stable housekeeping genes according to microarrays (RPL27 and GAPDH) was used to compare the gene expression between the groups (Hellemans, Mortier, De Paepe, Speleman, & Vandesompele, 2007). p-Values were calculated using an unpaired t test and 95% confidence intervals were calculated in R assuming unequal variances.

| Subjects' characteristics
From the 36 women with osteoporosis, 24 presented VFs and 12 had NVFs. No difference between the two groups (VF vs. NVF) was observed regarding the lumbar spine, femoral neck, and total hip BMD, frequency of chronic disease, life habits such as smoking, alcoholism, physical activity, and the number of falls. No differences were also observed in biochemical and bone turnover markers (calcium, phosphorus, PTH, alkaline phosphatase, 25-hydroxyvitamin D, P1NP, and CTX). Table 1 presents the participants' characteristics.

| Identification of VFs
Vertebral fractures were detected using Vertebral Fracture Assessment, using the Genant criteria (Genant & Jergas, 2003). From 24 women, 20 presented VFs grade 2 and four women presented VF grade 3. The VFs were found in 72% participants at the thoracic level and 28% at the lumbar level, with biconcave-shape in 49% and wedge-shape in 51%.

| Interaction network analysis
An interaction network constructed from all DEGs showed significant enrichment of interactions between the analyzed genes (PPI enrichment p-value = .00121). Among interacting DEG, we observed especially significant functional enrichment (p FDR corrected ≤.05) within one network containing eight downregulated genes. This enrichment was observed for partially overlapped ubiquitin signaling (UBXN6, CDC34, and UBE2M), adaptive immune system (CDC34, UBE2M, and PRR5), and antigen processing (CDC34 and UBE2M). Other enriched processes associated with interacting DEGs were redoxin functions (PRDX5 and PRDX2) and extracellular matrix binding (ITGA6 and CD248). All of the genes associated with those processes also showed expression in the muscle tissue. PRDX5 and UB2M also showed expression in bone tissue according to StringApp (Figure 2).

| Validation of differential gene expression using qPCR
For qPCR validation, we selected the eight DEGs with the highest β ratios which were also highlighted in enrichment and/or interaction network analysis. We analyzed six upregulated mRNAs (SNTG2, TRAF3IP2, PNO1, CD248, TNXB, and ITGA6) and two downregulated mRNA (PRDX5 and UBXN6). To determine the expression levels of these genes, we performed qPCR (VF: n = 19, Control: n = 12). We chose RPL27 and GAPDH as the endogenous control gene for the real-time qPCR because its expression levels in microarrays were very consistent and stable across all the samples. Six DEGs (SNTG2; CD248; TRAF3IP2; PNO1; ITGA6; UBXN6) showed the direction of the regulation of gene expression consistent with microarray results. Three of them, SNTG2, ITGA6, and TRAF3IP2, showed statistically significant (p < .05) changes in expression in the qPCR analysis. CD248 was almost significant at the 0.05 level (p = .0588). Expression changes of the two DEGs (PRDX5 and TNXB) were not confirmed. The results of qPCR analysis are shown in Figure 3 and in Table 3.

| ROC curve qPCR-RT
We analyzed the sensitivity and specificity of genes for the classification of VF groups. The area under curve (AUC) for three validated genes was significant (ITGA6 = 0.68, p = .05, SNTG2 = 0.71, p = .01, TRAF3IP2 = 0.71, p = .01). These results are shown in Figure 4.

| DISCUSSION
This study demonstrated for the first time the overexpression of transcripts related to VFs in older women with osteoporosis. We found that overexpression of SNTG2 strongly correlates with VF. SNTG2 belongs to the syntrophin family, which is highly expressed in the musculoskeletal system (Moon, Choi, Heo, Kim, & Kim, 2017). Syntrophins are cytoplasmic peripheral membrane proteins that bind to workings of mechanosensitive sodium channels and the extreme carboxyl-terminal domain of dystrophin. SNTG2 is mainly present in the neuromuscular junction (Inoue et al., 2008). Our data suggest a potential interaction between the muscular and skeletal systems in the genesis of osteoporotic VFs.
Moreover, we also found that TRAF3IP2, ITGA6, and CD248 transcripts presented high β ratios and upregulated expression.
Interestingly, TRAF3IP2 interacts with TNF receptor-associated factor 3, and is a central regulator of cytokine F I G U R E 2 Interaction network of the genes differentially expressed (p < .01) between the older women with osteoporosis and vertebral fracture (VF) and osteoporosis with no vertebral fractures (NVFs). On the figure are shown genes with at least one neighbor and differentially expressed interactors selected for qPCR validation. qPCR, quantitative polymerase chain reaction F I G U R E 3 Results for the genes selected for qPCR-RT (quantitative realtime polymerase chain reaction). *Genes Validated. qPCR, quantitative polymerase chain reaction production and type I interferons, and TRAF3 plays a role in various autoimmune and inflammatory diseases (Boyce, 2018;Tseng et al., 2010). Furthermore, TRAF3IP2 is highlighted in regulating RANKL-induced osteoclast development (Xiu et al., 2014) and its expression is altered in muscle tissue through the aging process according to the GTEX database, highlighting the possible role in the pathophysiology of VFs.
Regarding to ITGA6, the network analysis showed that this gene is an interactor with other upregulated gene, TNXB, associated with focal adhesion and integrin binding. Integrins may contribute to osteoclasts/bone adhesion and increase the resorption of the bone matrix, and therefore increase fracture risk. Landowski at al showed that inhibition of ITGA6 cleavage in vivo decreases osteolytic tumor activity and makes a sclerotic reaction in bone lesions (Landowski et al., 2014). The α6 integrin (ITGA6) is also expressed by human hematopoietic stem cell populations with plentiful laminin 511 (LN511) in healthy human bone (Qian et al., 2007). LN511 is a crucial component of the stem cell niche and maintains stem cell survival free of additional cytokines or growth factors (Graf, Heimfeld, & Torok-Strob, 2001). Both genes were shown in the enrichment analysis as a component of the extracellular matrix, which may function by forming structure or by linking connective tissue cells.
We also found that mTOR1 signaling pathway with two genes PNO1 and PDK1 was significantly enriched in patients with VFs. In fact, Glantschniget al. showed that through the mTOR pathway, RANKL, TNF, and M-CSF promote the survival of osteoclasts, inhibiting the apoptosis of osteoclasts (Glantschnig, Fisher, Wesolowski, Rodan, & Reszka, 2003). It is well known, that apoptosis of bone-related cells, including osteoclasts, osteoblasts, and osteocytes, is an important mechanism for bone mass regulation and osteoporosis. Besides, T A B L E 3 Summary of the results for the genes selected for qPCR-RT (quantitative real-time polymerase chain reaction) Gene symbol the bisphosphonates can prevent pathologic bone resorption by making osteoclast apoptosis (Hughes & Boyce, 1997) and experimental studies demonstrated that the rapamycin (an inhibitor of mTOR pathway) increased bone mass by inhibiting osteoclast differentiation (Dai et al., 2017). Furthermore, in mice with senile osteoporosis, rapamycin could restore the biological properties of aged bone marrow-derived mesenchymal stem cells (BMMSCs) by increasing osteogenic differentiation and proliferation (Ma et al., 2018). In fact, the authors found that activation of autophagy reestablished bone loss in aged mice, proposing that autophagy mechanism plays a pivotal role in the aging of BMMSCs, and activation of autophagy could partially converse the aging (Ma et al., 2018). The strength of our work is the comparison focused only on patients with osteoporosis (T-score ≤−2.5), with and without VFs, since VF is the most prevalent clinical outcome in osteoporosis. In fact, the mainstream of genetic studies of bone has been showed studies estimates heritability concentrated on BMD measurements. Those differences may differ in association and magnitude of effect at different skeletal sites. It is known that about 50% of women with osteoporosis with low BMD will not develop a fracture in the future. Some studies showed that not only BMD was hereditary, but also the bone fracture risk (Oei et al., 2014). In our study, we found a result of gene expression regardless of BMD.
Furthermore, the participants of this study were a very similar group selected from the Brazilian community cohort from Sao Paulo (SPAH). Lack of selection bias enables those findings to be extrapolated to the general Brazilian population. We also excluded mild-grade of spine deformities since some authors believe that this grade was not an actual VF, which was a limitation of previous VF study (Grados et al., 2009).
The limitations of this study include: First, the sample size was relatively small, and the data from this study should be replicated in large-scale works and in other phenotypes with different races or from different countries. The consistency of our results might seem limited by the small number of studied individuals, but the high power of the statistical analysis value of 0.9 of the data should give cause for assurance in the accuracy of our results, which were satisfactory to find DEGs (VanIterson, 2009). Second, since we used total blood cells that do not express all bone cell proteins, some protein translation changes may have gone unnoticed during our study. However, research using whole blood mRNA has the advantage of being a less invasive method and can measure the gene expression of systemic diseases, making it more useful and safer to use in clinical practice. Also mRNA reflects the functional state of cells, and integrates responses to both genetic and epigenetic factors of gene regulation, creating it a hopeful way to explore the disease evolution (Mohr & Liew, 2007;Riedmaier & Pfaff, 2013). Other limitations of the study were very low fold changes in gene expression. We explain this effect as associated with the similar phenotype of both analyzed groups. Usage of relatively big cohorts and logistic regression enabled us to identify even very light shifts in expression. Additional validation with qPCR enabled us to select the most promising transcripts associated with the study. qPCR confirmed the direction of expression changes for most of the transcripts supporting the significance of the obtained microarray results.
In conclusion, our data enabled us to identify and validate expression changes of the transcripts SNTG2, TRAF3IP2, and ITGA6 associated with osteoporotic VFs in elderly women, independently of BMD. These findings may serve as the basis for future longitudinal studies in the search for biomarkers of VFs, helping the early diagnosis of this comorbidity. In addition, the study of these molecules SNTG2, TRAF3IP2, and ITGA6 can be tested as therapeutic targets in osteoporosis-related VF.