# Pharmacogenomics Study for Raloxifene in Postmenopausal Female with Osteoporosis

Authors: Hsing-Fang Lu; Po-Hsin Chou; Gan-Hong Lin; Wan-Hsuan Chou; Shih-Tien Wang; Wirawan Adikusuma; Eko Mugiyanto; Kuo-Sheng Hung; Wei-Chiao Chang
Journal: Disease Markers
 (2020)

---

## Abstract

Osteoporosis is characterized by decreased bone mineral density and increased risk of fracture. Raloxifene is one of the treatments of osteoporosis. However, the responses were variable among patients. Previous studies revealed that the genetic variants are involved in the regulation of treatment outcomes. To date, studies that evaluate the influence of genes across all genome on the raloxifene treatment response are still limited. In this study, a total of 41 postmenopausal osteoporosis patients under regular raloxifene treatment were included. Gene-based analysis using MAGMA was applied to investigate the genetic association with the bone mineral density response to raloxifene at the lumbar spine or femoral neck site. Results from gene-based analysis indicated several genes (GHRHR, ABHD8, and TMPRSS6) related to the responses of raloxifene. Besides, the pathways of iron ion homeostasis, osteoblast differentiation, and platelet morphogenesis were enriched which implies that these pathways might be relatively susceptible to raloxifene treatment outcome. Our study provided a novel insight into the response to raloxifene.


---

## Body

## 1. Introduction

Osteoporosis is defined as the weakened architecture of the bone and elevated risk of fracture [1]. In Taiwan, the prevalence of osteoporosis is about 6.9% in males over the age of 65 and 21.2% in postmenopausal females [2]. Several risk factors of osteoporosis have been identified, including menopause, dairy intake, life style, and genetic factors [3]. Osteoporosis is a polygenic disease that is determined by multiple genetic effects from several genes. Variants near WNT16 and FONG have been reported to be associate with bone mineral density (BMD) [4–7]. The risk of osteoporosis increases dramatically in the elderly population, that bone fractures are known to increase mortality and lead to an enormous healthcare burden on society [8, 9]. Since the aged population is growing rapidly, osteoporosis has become a crucial clinical issue.Nonsteroidal selective estrogen-receptor modulators (SERMs) can inhibit bone resorption. SERMs have been widely used in the prevention and treatment of postmenopausal osteoporosis [10]. Raloxifene is a SERM agent used in the clinic for postmenopausal osteoporosis [11]. The drug is a partial agonist for the estrogen receptor and is known to inhibit proliferation of breast epithelium, while preserving bone mineral density (BMD). This drug increases BMD and reduces the occurrences of new vertebral fractures [12]. However, the therapeutic response of raloxifene varies widely among patients. Previous studies revealed that two polymorphisms of ESR1 (estrogen receptor 1) (rs9340799 and rs2234693) are associated with the BMD response to raloxifene in osteoporosis patients [13, 14]. UGT1A8 (UDP glucuronosyltransferase family 1 member A8) encodes an enzyme involved in extrahepatic glucuronidation of raloxifene [15], and UGT1A8∗2 was found to be associated with more raloxifene glucuronide metabolite formation in human jejunum tissue [16]. Labad et al. also identified a missense variant of UGT1A8, rs1042597, that associated with the treatment responses to raloxifene in postmenopausal schizophrenic women [17].Previous studies focused on the association between therapeutic responses and variants of a single gene. This type of study design called a candidate gene approach, which is a hypothesis-driven method that tests a small number of variants [18]. Limited genetic variants are unable to discover novel mechanisms. On the other hand, genome-wide association study (GWAS) is a hypothesis-free or hypothesis-generating method that provides opportunities to identify novel insights into the phenotype. In addition, GWAS is able to explore new genetic factors of complex human traits [19]. Since the mechanism of how genes influence response of raloxifene remains unclear, the aim of the current study is using gene-based analysis to identify the main network correlated with the responses of raloxifene.


## 2. Materials and Methods

### 2.1. Sample Collection

We collected samples from the neurosurgery outpatient clinic in Wang-Fang Hospital, Taipei, Taiwan. Our inclusion criteria were as follows: (1) 60 to 85-year-old postmenopausal female, (2)T‐score≦−2.5, and (3) under regular raloxifene treatment (60 mg/day). Patients with continuous steroid use (of over 6 months) or long-term inflammatory disease were excluded. In total, 41 osteoporosis patients were recruited. Bone mineral density (BMD) was measured by dual-energy radiograph absorptiometry (Lunar Prodigy, version 9.1, GE Healthcare, Madison, WI, USA) with standard protocols at the lumbar spine (LS; L2-4 or L1-4) and femoral neck (FN). Age, gender, body mass index (BMI), medication profile, and follow-up duration of all individual were collected. We also collected the BMD profiles after patients receiving treatment for 1-2 years. The study was approved by the Joint Institutional Review Board of Taipei Medical University. All subjects provided written informed consent.


### 2.2. DNA Extraction and Genotyping Array

Peripheral whole blood was collected from the patients. Gentra Puregene Blood Kit (Qiagen, Valencia, CA, USA) was used to extract the genomic DNA. The extraction process was according to the user manual of the product. The quality of DNA was confirmed by 1% agarose gel electrophoresis. Sample with no obvious degradation can be submitted to further analysis. Affymetrix Axiom TWB SNP Array (Thermo Fisher Scientific, USA) was used to perform GWAS. The process was done by the National Center for Genome Medicine, Taiwan (NCGM).


### 2.3. Quality Control

PLINK [20] was applied to conduct the quality control of genotype data. We excluded SNPs with low genotyping call rate (<0.95), minor allele frequency less than 0.05, and P value of Hardy-Weinberg equilibrium less than 1×10-6. Individuals with a sample call rate less than 0.95 were excluded. The imputation task was done by Michigan Imputation Server (Minimac3). Data from the Haplotype Reference Consortium (HRC) was used as reference panels for imputation in order to increase statistical power. After imputation, only variants which MAF≧0.05 and R2≧0.3 were included in the association analysis. Principle components (PCs) were generated by PLINK with samples from HapMap3. A linear regression model was used for statistically testing the associations between SNPs and phenotypes by PLINK2. We defined the percentage change of bone mineral density from baseline by the following equation:
(1)x=2ndBMD−1stBMD1stBMD.We further normalized phenotype before association analysis. We included age, BMI, follow-up months, and PC1 to PC10 as covariates in the analysis model. The Manhattan plot and qq plot were visualized by the R package “qqman,” and the suggestive significant threshold was1×10−5 as the default setting of this package. The regional plots were generated by LocusZoom [21].


### 2.4. Variant Annotation

We annotated the rs number of candidate variates by ANNOVAR [22] with avsnp150 database. The annotation of nearby genes was according to the NCBI RefSeq hg19 reference panel. In order to compare our results with previous studies, the data of GWAS catalog v1.0.2 (20190322) was downloaded from GWAS catalog website.


### 2.5. Gene-Based Analysis and Gene-Set Analysis

Since our sample size was very limited, the effects of individual markers may be too weak to detect. Therefore, we applied MAGMA to conduct gene-based analysis to increase the power [23]. This method can analyze multiple genetic markers together to determine their joint effect. Input SNPs were mapped to 18306 protein coding genes; therefore, the Bonferroni significant threshold for gene analysis was using P<2.73×10−6. To understand the biological enriched pathways related with the raloxifene treatment response, we used gene-set analysis by MAGMA. The curated gene sets and GO terms were obtained from MsigDB v6.1. Curated gene sets were collected from 9 data resources including KEGG, Reactome, BioCarta, and biomedical literatures. Variates located within transcription start region and stop region of the gene region will be mapped to that gene. We used 1000 genomes phase 3 as reference genome.


## 2.1. Sample Collection

We collected samples from the neurosurgery outpatient clinic in Wang-Fang Hospital, Taipei, Taiwan. Our inclusion criteria were as follows: (1) 60 to 85-year-old postmenopausal female, (2)T‐score≦−2.5, and (3) under regular raloxifene treatment (60 mg/day). Patients with continuous steroid use (of over 6 months) or long-term inflammatory disease were excluded. In total, 41 osteoporosis patients were recruited. Bone mineral density (BMD) was measured by dual-energy radiograph absorptiometry (Lunar Prodigy, version 9.1, GE Healthcare, Madison, WI, USA) with standard protocols at the lumbar spine (LS; L2-4 or L1-4) and femoral neck (FN). Age, gender, body mass index (BMI), medication profile, and follow-up duration of all individual were collected. We also collected the BMD profiles after patients receiving treatment for 1-2 years. The study was approved by the Joint Institutional Review Board of Taipei Medical University. All subjects provided written informed consent.


## 2.2. DNA Extraction and Genotyping Array

Peripheral whole blood was collected from the patients. Gentra Puregene Blood Kit (Qiagen, Valencia, CA, USA) was used to extract the genomic DNA. The extraction process was according to the user manual of the product. The quality of DNA was confirmed by 1% agarose gel electrophoresis. Sample with no obvious degradation can be submitted to further analysis. Affymetrix Axiom TWB SNP Array (Thermo Fisher Scientific, USA) was used to perform GWAS. The process was done by the National Center for Genome Medicine, Taiwan (NCGM).


## 2.3. Quality Control

PLINK [20] was applied to conduct the quality control of genotype data. We excluded SNPs with low genotyping call rate (<0.95), minor allele frequency less than 0.05, and P value of Hardy-Weinberg equilibrium less than 1×10-6. Individuals with a sample call rate less than 0.95 were excluded. The imputation task was done by Michigan Imputation Server (Minimac3). Data from the Haplotype Reference Consortium (HRC) was used as reference panels for imputation in order to increase statistical power. After imputation, only variants which MAF≧0.05 and R2≧0.3 were included in the association analysis. Principle components (PCs) were generated by PLINK with samples from HapMap3. A linear regression model was used for statistically testing the associations between SNPs and phenotypes by PLINK2. We defined the percentage change of bone mineral density from baseline by the following equation:
(1)x=2ndBMD−1stBMD1stBMD.We further normalized phenotype before association analysis. We included age, BMI, follow-up months, and PC1 to PC10 as covariates in the analysis model. The Manhattan plot and qq plot were visualized by the R package “qqman,” and the suggestive significant threshold was1×10−5 as the default setting of this package. The regional plots were generated by LocusZoom [21].


## 2.4. Variant Annotation

We annotated the rs number of candidate variates by ANNOVAR [22] with avsnp150 database. The annotation of nearby genes was according to the NCBI RefSeq hg19 reference panel. In order to compare our results with previous studies, the data of GWAS catalog v1.0.2 (20190322) was downloaded from GWAS catalog website.


## 2.5. Gene-Based Analysis and Gene-Set Analysis

Since our sample size was very limited, the effects of individual markers may be too weak to detect. Therefore, we applied MAGMA to conduct gene-based analysis to increase the power [23]. This method can analyze multiple genetic markers together to determine their joint effect. Input SNPs were mapped to 18306 protein coding genes; therefore, the Bonferroni significant threshold for gene analysis was using P<2.73×10−6. To understand the biological enriched pathways related with the raloxifene treatment response, we used gene-set analysis by MAGMA. The curated gene sets and GO terms were obtained from MsigDB v6.1. Curated gene sets were collected from 9 data resources including KEGG, Reactome, BioCarta, and biomedical literatures. Variates located within transcription start region and stop region of the gene region will be mapped to that gene. We used 1000 genomes phase 3 as reference genome.


## 3. Results

### 3.1. Baseline Characteristics

The distribution of the BMD change percent is shown in Figure1. After treatment, the BMD response to raloxifene was diverse at both lumbar spine and femoral neck sites in patients. In total, 41 patients were assessed for femoral neck BMD (S Table 1), and the average percentage of BMD change was −2.5±15.8%. The mean follow-up time was 35.61±23.54 months among all subjects. The mean age of subjects was 72.4±7.0 years, and the average BMI was 23.96±3.49 kg/m2. Among the 41 patients assessed for FN BMD, 34 were also assessed for lumbar spine BMD. The average age and BMI of this subset were similar to the complete group.Figure 1
The distribution of the BMD change % at the lumbar spine and femoral neck site.


### 3.2. Association Analysis

We conducted a genome-wide association analysis to reveal the risk loci associated with raloxifene treatment response. A principal component analysis showed that there was no population stratification (S Figure1). Neither site revealed genome-wide significant association signal of BMD change (S Figures 2, 3). However, there were 112 variants with suggestive significance for LS percent BMD change, while 13 variants showed significance for FN percentage BMD change. The most significant variant in the LS site was rs7768089 (P=2.29×10−6), which was an intronic SNP of FUT9 located on 6q16.1. Meanwhile, rs34311394 at 12q21.2 was the most significant SNP (P=4.44×10−6) for the FN site. This SNP is located in the intergenic region between CSRP2 and E2R7. There was no overlap variants among suggestively significant SNPs between the two sites.


### 3.3. Candidate Genes Cause Decreased BMD in Mutant Mice

Next, a gene-based analysis was conducted by MAGMA. However, none of the genes reached the threshold of statistical significance (Figure2). The top ten significant genes at the LS and FN sites are shown in Tables 1 and 2. To understand the potential biological functions of these genes, we queried the phenotypes of mutant mice in the International Mouse Phenotype Consortium (IMPC) database. Among the queried genes, GHRHR, ABHD8, and TMPRSS6 showed decreased BMD in a homozygous mutant mouse in both sexes.Figure 2
Manhattan plot of the gene-based test based on response to raloxifene. The red line represented suggestive significant threshold (P=2.73×10-6): (a) Lumbar spine site and (b) Femoral neck site.
(a)(b)Table 1
Top 10 significant genes of gene-based analysis at the lumbar spine site.
GeneCHRSTARTSTOPNumber of SNPsZSTATp valueKO mouse bone phenotype∗FUT9696463860966634884744.112.0.E−05NAAPBB1IP1026727132268567323043.9623.7.E−05NAMYO15B1773584139736229291043.8087.0.E−05NAUSHBP1191735998517393595543.7289.7.E−05Vertebral fusionDNAH5513690440139446526913.7121.0.E−04NoNRG251392263641394228842383.7091.0.E−04NAKIAA019681260365021261040822113.6911.1.E−04NAPHYH101331979613344412703.6251.4.E−04NoELTD1179355449794724033093.5671.8.E−04NAOR4K114204037672040484273.5452.0.E−04NACHR: chromosome; START: start position of the gene; STOP: stop position of the gene.∗Knockout mouse bone-related phenotypes were queried from International Mouse Phenotype Consortium.Table 2
Top 10 significant genes of gene-based analysis at the femoral neck site.
GeneCHRSTARTSTOPNumber of SNPsZSTATp valueKO mouse bone phenotype∗AC068987.112522037895220663613.9533.90E−05NAGHRHR730978284310328691393.8595.70E−05BMD decreaseF1013113777128113803843353.8366.20E−05NASCN8A1251984050522066482003.6741.20E−04NATMEM128442372694249950843.5621.80E−04NAMRPL34191740341817417652353.512.20E−04NAHSPA95137890571137911133143.5012.30E−04NAABHD8191740294017421045423.4942.40E−04BMD decreaseABLIM24796703981605594883.3943.40E−04NATMPRSS6223746147637505603973.383.60E−04BMD decreaseCHR: chromosome; START: start position of the gene; STOP: stop position of the gene.∗Knockout mouse bone-related phenotypes were queried from International Mouse Phenotype Consortium.


### 3.4. Enrichment of Osteoblast Differentiation Pathway

We further applied gene-set analysis to reveal the significant enriched gene sets. After FDR adjustment, although none of the pathway achieved the statistical significance (Tables3 and 4), iron homeostasis for the LS site (P=2×10−4) was the most significant gene set. Regarding the femoral neck site, the most significant gene set was platelet morphogenesis (P=1×10−4). In addition, another gene sets of “osteoblast differentiation by phenylamil up” showed nominal significance for the femoral neck site as well (P=4×10−4).Table 3
Gene-set analysis of raloxifene response at lumbar spine site.
DatabasePathwayNumber of genesBetaSEp valueq valueGO bpIron ion homeostasis630.3630.1012.00E−040.884GO bpRegulation of immunoglobulin secretion140.7740.2273.00E−040.884GO bpTransition metal ion homeostasis960.2670.084.00E−040.884Curated gene setsInhibition of replication initiation of damaged DNA by rb1 e2f1100.8790.2665.00E−040.884Curated gene setsApoptosis820.280.0866.00E−040.884GO bpResponse to protozoan200.5690.1776.00E−040.884Curated gene setsCdh1 targets 3 dn580.3310.1069.00E−040.884Curated gene setsTreating iron overload70.8790.2839.00E−040.884GO bpToll like receptor signaling pathway770.30.0991.00E−030.884GO bpCellular iron ion homeostasis400.3740.1241.00E−030.884Curated gene setsTlx targets 60 hr dn2580.1480.0491.00E−030.884GO bpNegative regulation of platelet activation170.6260.2112.00E−030.884Curated gene setsTgfb1 signaling via nfic 10 hr dn290.4640.1582.00E−030.884Curated gene setsResponse to tosedostat 24 hr dn9360.0740.0252.00E−030.884Curated gene setsMetabolism of amino acids and derivatives1860.1660.0572.00E−030.884GO bpCellular defense response540.3230.1112.00E−030.884Curated gene setsAkt phosphorylates targets in the cytosol120.5590.1932.00E−030.884Curated gene setsAsbcell pathway90.7830.2712.00E−030.884Curated gene setsCell cycle s1520.1930.0682.00E−030.884Curated gene setsAr tf pathway470.3520.1242.00E−030.884GO: gene ontology; bp: biological process;q values were from false discovery rate.Table 4
Gene-set analysis of raloxifene response at femoral neck site.
DatabasePathwayNumber of genesBetaSEp valueq valueGO bpPlatelet morphogenesis160.7620.2081.00E−040.859GO bpRegulation of nuclear transcribed mrna catabolic process deadenylation-dependent decay140.6330.1873.00E−040.859GO bpCellular modified amino acid catabolic process170.6160.1854.00E−040.859Curated gene setsOsteoblast differentiation by phenylamil up120.7790.2344.00E−040.859GO bpRegulation of pri mirna transcription from rna polymerase ii promoter160.6910.2095.00E−040.859GO bpRegulation of establishment of protein localization to mitochondrion1180.2340.0737.00E−040.966GO bpRegulation of protein targeting to mitochondrion890.2660.0837.00E−040.966Curated gene setsWerner syndrome and normal aging up910.2510.089.00E−040.991Curated gene setsAml methylation cluster 6 dn320.4180.1371.00E−030.999GO bpEpithelial structure maintenance210.5590.1841.00E−030.999GO bpRegulation of substrate adhesion-dependent cell spreading410.3790.1292.00E−030.999Curated gene setsProteolytic cleavage of snare complex proteins160.7520.262.00E−030.999GO bpCardiac myofibril assembly160.5740.2042.00E−030.999GO bpPositive regulation of pri mirna transcription from rna polymerase ii promoter100.7020.252.00E−030.999Curated gene setsBotulinum neurotoxicity180.6740.2413.00E−030.999GO bpNegative chemotaxis380.3910.1423.00E−030.999GO bpRegulation of endocrine process440.3480.1263.00E−030.999Curated gene setsChondroitin sulfate biosynthesis170.5070.1843.00E−030.999GO bpCellular metabolic compound salvage350.360.1313.00E−030.999GO bpCytokinesis800.2460.093.00E−030.999GO: gene ontology; bp: biological process;q values were from false discovery rate.


### 3.5. Comparison with Previous Candidate Variants

Since some variants have been reported to be associated with response to raloxifene, we next focused two SNPs withinESR1 (rs9340799 and rs2234693) and one SNP in UGT1A8 (rs1042597). As shown in Supplementary Table 2, none of the SNPs demonstrated nominal significance (P>0.05) at the LS or FN site. We further extracted larger regions around these two genes (±500 Mb) of two phenotypes and generated regional plots. However, no variant within these regions reached suggestive significance (S Figure 4-7).


## 3.1. Baseline Characteristics

The distribution of the BMD change percent is shown in Figure1. After treatment, the BMD response to raloxifene was diverse at both lumbar spine and femoral neck sites in patients. In total, 41 patients were assessed for femoral neck BMD (S Table 1), and the average percentage of BMD change was −2.5±15.8%. The mean follow-up time was 35.61±23.54 months among all subjects. The mean age of subjects was 72.4±7.0 years, and the average BMI was 23.96±3.49 kg/m2. Among the 41 patients assessed for FN BMD, 34 were also assessed for lumbar spine BMD. The average age and BMI of this subset were similar to the complete group.Figure 1
The distribution of the BMD change % at the lumbar spine and femoral neck site.


## 3.2. Association Analysis

We conducted a genome-wide association analysis to reveal the risk loci associated with raloxifene treatment response. A principal component analysis showed that there was no population stratification (S Figure1). Neither site revealed genome-wide significant association signal of BMD change (S Figures 2, 3). However, there were 112 variants with suggestive significance for LS percent BMD change, while 13 variants showed significance for FN percentage BMD change. The most significant variant in the LS site was rs7768089 (P=2.29×10−6), which was an intronic SNP of FUT9 located on 6q16.1. Meanwhile, rs34311394 at 12q21.2 was the most significant SNP (P=4.44×10−6) for the FN site. This SNP is located in the intergenic region between CSRP2 and E2R7. There was no overlap variants among suggestively significant SNPs between the two sites.


## 3.3. Candidate Genes Cause Decreased BMD in Mutant Mice

Next, a gene-based analysis was conducted by MAGMA. However, none of the genes reached the threshold of statistical significance (Figure2). The top ten significant genes at the LS and FN sites are shown in Tables 1 and 2. To understand the potential biological functions of these genes, we queried the phenotypes of mutant mice in the International Mouse Phenotype Consortium (IMPC) database. Among the queried genes, GHRHR, ABHD8, and TMPRSS6 showed decreased BMD in a homozygous mutant mouse in both sexes.Figure 2
Manhattan plot of the gene-based test based on response to raloxifene. The red line represented suggestive significant threshold (P=2.73×10-6): (a) Lumbar spine site and (b) Femoral neck site.
(a)(b)Table 1
Top 10 significant genes of gene-based analysis at the lumbar spine site.
GeneCHRSTARTSTOPNumber of SNPsZSTATp valueKO mouse bone phenotype∗FUT9696463860966634884744.112.0.E−05NAAPBB1IP1026727132268567323043.9623.7.E−05NAMYO15B1773584139736229291043.8087.0.E−05NAUSHBP1191735998517393595543.7289.7.E−05Vertebral fusionDNAH5513690440139446526913.7121.0.E−04NoNRG251392263641394228842383.7091.0.E−04NAKIAA019681260365021261040822113.6911.1.E−04NAPHYH101331979613344412703.6251.4.E−04NoELTD1179355449794724033093.5671.8.E−04NAOR4K114204037672040484273.5452.0.E−04NACHR: chromosome; START: start position of the gene; STOP: stop position of the gene.∗Knockout mouse bone-related phenotypes were queried from International Mouse Phenotype Consortium.Table 2
Top 10 significant genes of gene-based analysis at the femoral neck site.
GeneCHRSTARTSTOPNumber of SNPsZSTATp valueKO mouse bone phenotype∗AC068987.112522037895220663613.9533.90E−05NAGHRHR730978284310328691393.8595.70E−05BMD decreaseF1013113777128113803843353.8366.20E−05NASCN8A1251984050522066482003.6741.20E−04NATMEM128442372694249950843.5621.80E−04NAMRPL34191740341817417652353.512.20E−04NAHSPA95137890571137911133143.5012.30E−04NAABHD8191740294017421045423.4942.40E−04BMD decreaseABLIM24796703981605594883.3943.40E−04NATMPRSS6223746147637505603973.383.60E−04BMD decreaseCHR: chromosome; START: start position of the gene; STOP: stop position of the gene.∗Knockout mouse bone-related phenotypes were queried from International Mouse Phenotype Consortium.


## 3.4. Enrichment of Osteoblast Differentiation Pathway

We further applied gene-set analysis to reveal the significant enriched gene sets. After FDR adjustment, although none of the pathway achieved the statistical significance (Tables3 and 4), iron homeostasis for the LS site (P=2×10−4) was the most significant gene set. Regarding the femoral neck site, the most significant gene set was platelet morphogenesis (P=1×10−4). In addition, another gene sets of “osteoblast differentiation by phenylamil up” showed nominal significance for the femoral neck site as well (P=4×10−4).Table 3
Gene-set analysis of raloxifene response at lumbar spine site.
DatabasePathwayNumber of genesBetaSEp valueq valueGO bpIron ion homeostasis630.3630.1012.00E−040.884GO bpRegulation of immunoglobulin secretion140.7740.2273.00E−040.884GO bpTransition metal ion homeostasis960.2670.084.00E−040.884Curated gene setsInhibition of replication initiation of damaged DNA by rb1 e2f1100.8790.2665.00E−040.884Curated gene setsApoptosis820.280.0866.00E−040.884GO bpResponse to protozoan200.5690.1776.00E−040.884Curated gene setsCdh1 targets 3 dn580.3310.1069.00E−040.884Curated gene setsTreating iron overload70.8790.2839.00E−040.884GO bpToll like receptor signaling pathway770.30.0991.00E−030.884GO bpCellular iron ion homeostasis400.3740.1241.00E−030.884Curated gene setsTlx targets 60 hr dn2580.1480.0491.00E−030.884GO bpNegative regulation of platelet activation170.6260.2112.00E−030.884Curated gene setsTgfb1 signaling via nfic 10 hr dn290.4640.1582.00E−030.884Curated gene setsResponse to tosedostat 24 hr dn9360.0740.0252.00E−030.884Curated gene setsMetabolism of amino acids and derivatives1860.1660.0572.00E−030.884GO bpCellular defense response540.3230.1112.00E−030.884Curated gene setsAkt phosphorylates targets in the cytosol120.5590.1932.00E−030.884Curated gene setsAsbcell pathway90.7830.2712.00E−030.884Curated gene setsCell cycle s1520.1930.0682.00E−030.884Curated gene setsAr tf pathway470.3520.1242.00E−030.884GO: gene ontology; bp: biological process;q values were from false discovery rate.Table 4
Gene-set analysis of raloxifene response at femoral neck site.
DatabasePathwayNumber of genesBetaSEp valueq valueGO bpPlatelet morphogenesis160.7620.2081.00E−040.859GO bpRegulation of nuclear transcribed mrna catabolic process deadenylation-dependent decay140.6330.1873.00E−040.859GO bpCellular modified amino acid catabolic process170.6160.1854.00E−040.859Curated gene setsOsteoblast differentiation by phenylamil up120.7790.2344.00E−040.859GO bpRegulation of pri mirna transcription from rna polymerase ii promoter160.6910.2095.00E−040.859GO bpRegulation of establishment of protein localization to mitochondrion1180.2340.0737.00E−040.966GO bpRegulation of protein targeting to mitochondrion890.2660.0837.00E−040.966Curated gene setsWerner syndrome and normal aging up910.2510.089.00E−040.991Curated gene setsAml methylation cluster 6 dn320.4180.1371.00E−030.999GO bpEpithelial structure maintenance210.5590.1841.00E−030.999GO bpRegulation of substrate adhesion-dependent cell spreading410.3790.1292.00E−030.999Curated gene setsProteolytic cleavage of snare complex proteins160.7520.262.00E−030.999GO bpCardiac myofibril assembly160.5740.2042.00E−030.999GO bpPositive regulation of pri mirna transcription from rna polymerase ii promoter100.7020.252.00E−030.999Curated gene setsBotulinum neurotoxicity180.6740.2413.00E−030.999GO bpNegative chemotaxis380.3910.1423.00E−030.999GO bpRegulation of endocrine process440.3480.1263.00E−030.999Curated gene setsChondroitin sulfate biosynthesis170.5070.1843.00E−030.999GO bpCellular metabolic compound salvage350.360.1313.00E−030.999GO bpCytokinesis800.2460.093.00E−030.999GO: gene ontology; bp: biological process;q values were from false discovery rate.


## 3.5. Comparison with Previous Candidate Variants

Since some variants have been reported to be associated with response to raloxifene, we next focused two SNPs withinESR1 (rs9340799 and rs2234693) and one SNP in UGT1A8 (rs1042597). As shown in Supplementary Table 2, none of the SNPs demonstrated nominal significance (P>0.05) at the LS or FN site. We further extracted larger regions around these two genes (±500 Mb) of two phenotypes and generated regional plots. However, no variant within these regions reached suggestive significance (S Figure 4-7).


## 4. Discussion

This study used two skeletal positions to assess the relationships between BMD response to raloxifene and gene network. Although statistical power was limited sample size, several candidates includingGHRHR, ABHD8, and TMPRSS6 still addressed significance after gene-based analysis. These three genes were also associated with decreased BMD in mutant mouse models. From gene-set analysis, we found that iron homeostasis and osteoblast differentiation gene sets were enriched. Previous studies indicated that iron overload affects the bone phenotype and increases risk of bone fracture. Iron overload is an important factor that is involved in the activation of osteoclast, the differentiation process, and the reduction of osteoblast formation [24, 25]. Besides, in patients with thalassemia major, which is often characterized by low BMD values and iron overload, the circulating IGF-1 levels and its binding protein are lower than those in the healthy individuals [26]. Thus, IGF-1 is critical in bone remodeling. The lower level of IGF-1 may increase the risk of hip fracture [27]. Taken together, iron overload may increase bone resorption and decrease bone formation through osteoblast and osteoclast activity or through the endocrine growth hormone- (GH-) IGF-1 axis pathway.Heilberg et al. reported that the rs9340799 (Xbal) AA genotype and rs2234693 (PvuII) CC genotype of theESR1 have higher LS bone density after one-year treatment of raloxifene in Hispanic women [14]. Conversely, Mondockova et al. revealed that the AA genotype of rs9340799 and TT genotype of rs2234693 showed worse treatment response to raloxifene at the LS site in Slovak women [13]. Because of the inconsistent findings, we tried to confirm it by using this study. However, none of the variants were associated with the response to raloxifene in the current study. In another study, individuals with the CC genotype of the missense SNP (Ala173Gly) of UGT1A8 revealed more improvement in negative symptoms after raloxifene treatment [17]. However, this variant cannot be replicated in the current study as well. We found that the most significant association with drug response at LS site was an intronic SNP of fucosyltransferase 9 (FUT9), rs7768089. This gene is a member of the glycosyltransferase family, and it participates in the biosynthesis of Lewis antigen. HaploReg v4.1 demonstrated that rs7768089 may alter the transcription factor (Foxa disc2, Nrd-2 3, and TCF11::MafG) binding affinity, and according to GTEx, this SNP is an expression quantitative locus (eQTL) of FUT9 in the pancreas and brain tissue (S Table 3). However, the mechanistic relationship between FUT9 and BMD remains to be clarified.For the FN site, the most significant SNP (rs34311394) was located within the intron betweenCSRP2 and E2F7. CSRP2 (cysteine and glycine rich protein 2) is involved in development and cellular differentiation, and E2F7 (E2F transcription factors) participates in the regulation of cell cycle progression. This SNP is associated with the expression of CSRP2 in the brain tissue (S Table 3). The most significant association revealed by a gene-based study was GHRHR (growth hormone-releasing hormone receptor). Although no traits related to BMD or drug response have been reported, the Ghrhrtm1.1(KOMP)Vlcg mouse model showed decreased bone mineral density according to the IMPC database. Based on the results in GWAS catalog, the genetic variants of GHRHR associated with congenital left sided heart lesions and adolescent idiopathic scoliosis. Further functional studies are needed to determine the plausible impact of these candidate genes on the bone.Due to its estrogenic effects, raloxifene may increase the risk of thromboembolic events, such as deep vein thrombosis and pulmonary embolism. Estrogen increases the concentration of plasma fibrinogen and coagulation factors and enhances the aggregation of platelets [28]. From the gene-set analysis for the femoral neck site, the most significant treatment-associated gene-set was platelet morphogenesis. These results hint that genetic variants could possibly play a role in the occurrence of side effects. Thus, further studies may focus on the pharmacogenomics of thromboembolic events in raloxifene users. The main limitation of the current study is the sample size which is not able to confer sufficient statistical power to reach genome-wide significance. Larger cohorts are needed to replicate our results. Another limitation is the in the study. The majority of subjects are elderly patients with a mean age of 72.4 years. Previous studies have reported that bone loss may accelerate in the individuals who are over 70 years old [29]. In addition, aging influences GH and IGF-1 serum levels [27], which might result in the change of BMD levels and the pharmacological effects. Further studies should also consider the age, the concentration of raloxifene in the blood, and the drug compliance for each patient.


## 5. Conclusions

In summary, we found plausible gene sets in iron ion homeostasis, osteoblast differentiation, and platelet morphogenesis pathways that may influence raloxifene responses. These findings offer new important insights and warrant further investigations into the pharmacogenomics of raloxifene.



---

*Source: 8855423-2020-08-31.xml*