Methylation of estrogen receptor 2 (ESR2) in deep paravertebral muscles and its association with idiopathic scoliosis

Idiopathic scoliosis (IS) is one of the most common spinal disorders in adolescents. Despite many studies, the etiopathogenesis of IS is still poorly understood. In recent years, the role of epigenetic factors in the etiopathogenesis of IS has been increasingly investigated. It has also been postulated that the development and progression of the disease is related to gender and puberty, and could be associated with estrogen action. Estrogen hormones act via estrogen receptor 1 (ESR1) and estrogen receptor 2 (ESR2). It has been suggested that ESR2 expression is dependent on methylation within its gene promoter. So far, no studies have evaluated local, tissue-specific DNA methylation in patients with IS. Thus, our study aimed to analyze the methylation and expression level of ESR2 in the paraspinal muscles of the convex and concave side of the IS curvature. The methylation level within ESR2 promoter 0N, but not exon 0N, was significantly higher on the concave side of the curvature compared to the convex side. There was no significant correlation between ESR2 expression and methylation level in the promoter 0N on the convexity of thoracic scoliosis, whereas, on the concave side of the curvature, we observed a moderate negative correlation. There was no difference in the methylation levels of the ESR2 promoter and exon 0N between groups of patients with Cobb angle ≤ 70° and > 70° on the concave and convex side of the curvature. We also found no statistically significant correlation between the Cobb angle value and the mean methylation level in either the ESR2 promoter or exon 0N on the convex or concave side of the curvature. Our findings demonstrate that DNA methylation at the ESR2 promoter in deep paravertebral muscle tissue is associated with the occurrence but not with the severity of idiopathic scoliosis.

Idiopathic scoliosis (IS) is one of the most common spinal disorders in adolescents, with a reported prevalence of 1-3%. IS is a three-dimensional deformity of the spinal column, resulting in curvature in the coronal plane, sagittal plane deviation, axial rotation of the vertebrae, rib prominence, and trunk imbalance. IS may remain stable throughout an individual's life or could progress to a very severe form. The risk of severe curve progression is associated with the growth of the spine 1 . However, at an early stage of the disease, it is difficult to predict the final phenotype. Severe idiopathic scoliosis leads to cardiorespiratory impairment due to lung restriction and results in significant morbidity. The degree of impairment has been shown to correlate with the severity of the spinal deformity [2][3][4] . Therefore, identifying patients who are at risk of developing scoliosis and severe curve progression is crucial in order to provide them with early treatment 1,5 . The idea that IS has a strong genetic background is supported by the higher prevalence of scoliosis in families with an affected member compared to the general population 6 . Candidate genes potentially associated with the occurrence and progression of IS have been identified in many studies, including evaluations of family linkage, genome-wide association studies (GWAS), SNPs (single nucleotide polymorphisms) and gene expression profiling 6 .
The occurrence and progression of IS is related to gender and puberty. The prevalence is much higher, and the curvatures are larger in female adolescents compared to males. Thus, it has been postulated that IS progression could be influenced by estrogen action [7][8][9] . Estrogen hormones act on target cells through two types of receptors: estrogen receptor 1 (ESR1) and estrogen receptor 2 (ESR2) 10,11 and through G protein-coupled estrogen receptor 1 (GPER), which is a membrane protein 12 . It was postulated that two SNPs of the ESR1 were associated with IS [13][14][15][16][17] . However, the relationship between ESR1 polymorphisms and IS remains ambiguous [18][19][20] . Meta-analysis of four studies performed by Chen et al. revealed a non-significant association between rs9340799 ESR1 polymorphism and AIS 21 . The initial results regarding the association of the ESR2 polymorphism rs1256120 with a susceptibility to occurrence and progression of IS were not confirmed in replication studies 19,22,23 . Nevertheless, there is a suggestion that other ESR2 polymorphisms may be associated with IS curvature progression to the severe forms 24 . Additionally, the association of GPER polymorphisms with severity of the curve in patients with idiopathic scoliosis was demonstrated. Abnormalities in GPER presence may influence the deterioration of spine deformity 12 . However, it was not confirmed in replication study 25 .
In recent years, an increasing number of reports have described the effect of estrogens on skeletal muscle 26,27 . ESR2 has been detected at both mRNA and protein levels in skeletal muscle 28 . ESR2 expression was also observed in paravertebral skeletal muscle tissue 26 .
However, the polymorphisms associated with IS are commonly present in the DNA of individuals without IS as well 6 . Therefore, it has been postulated that IS arises from both genetic and environmental factors. Based on observations from twin studies, Grauers et al. estimated that the risk of developing scoliosis is the result of an additive genetic effect in the minority, and the result of environmental effects in the majority of affected individuals 29 . Cheng et al. suggested that genetics play a more significant role during the initiation phase of adolescent idiopathic scoliosis (AIS), whereas during curve progression, environmental factors have a greater contribution 30 . The identification of the factors linking the genome and the environment in IS etiopathogenesis is crucial in the prevention and treatment of the disease. One mechanism often cited as a link between genetics and the environment is epigenetics. However, investigating the potential epigenetic mechanisms that underlie IS etiology is a relatively new avenue of research 31,32 . To date, only a few studies describing DNA methylation in IS were published [32][33][34][35][36] . Additionally, none of the studies published so far evaluated local, tissue-specific DNA methylation.
DNA methylation is an epigenetic DNA modification associated with the regulatory regions of some genes. It has been suggested that ESR2 expression in disease contexts is influenced by methylation at 0N and 0K gene promoters and their corresponding exons 37 . Differences in methylation levels between patients with different IS phenotypes could be associated with disease susceptibility or progression. Thus, we sought to analyze the methylation and expression levels of ESR2 in the paraspinal muscles of the convex and concave side of the IS curvature.
DNA methylation at the ESR2 promoter 0N. The mean methylation level of all CpG sites within ESR2 promoter 0N was significantly higher on the concave side of the curvature comparing to the convex side (3.94% ± 0.74% vs. 3.66% ± 0.7%; P = 0.02; Fig. 1A). The methylation level of seven CpG sites within ESR2 promoter 0N (each CpG analyzed separately) was significantly higher on the concave side of the curvature compared to the convex side (P < 0.05, Additional file 1: Table S1; Fig. 2).  Fig. 1B). Additionally, there was no significant difference in the methylation level at each CpG site (each CpG analyzed separately) between deep paravertebral muscles on the convex and concave side of the curvature (P > 0.05, Additional file 1: Table S2).
Correlation between ESR2 methylation levels and relative expression of the ESR2 gene. On the convex side of thoracic scoliosis, the correlation between ESR2 expression and mean methylation level in promoter 0N was weak, negative, and non-significant (Spearman's rank correlation, R = − 0.3; P = 0.2; Fig. 3A). In contrast, on the concave side of the curvature, a moderate negative correlation was observed (R = -0.5; P = 0.02; Fig. 3B). Additionally, a significant negative correlation between ESR2 mRNA expression and methylation level was observed at 3 and 11 CpG sites in promoter 0N on the convex and concave side of the curvature, respectively (R ranged from − 0.4 to − 0.6, P < 0.5; Fig. 4). In the case of the exon 0N region, the significant, negative correlation between ESR2 expression and methylation level was observed at only one CpG site on the convexity of thoracic scoliosis (R = − 0.4; P = 0.02; Fig. 5). The correlation between ESR2 expression and mean methylation level in exon 0N was weak, positive, non-significant on the convex side of the curvature (Spearman's rank correlation: R = 0.1; P = 0.6; Fig. 3C) and weak, negative, non-significant on the concavity of thoracic scoliosis (Spearman's rank correlation: R = -0.1; P = 0.6; Fig. 3D).

Discussion
The etiopathogenesis of idiopathic scoliosis remains poorly understood 6 . One of the biggest challenges of research into IS etiopathogenesis is that affected individuals possess no obvious structural deficiencies or malformations in the vertebral column and associated soft tissues, except the curvature of the spine 38 . www.nature.com/scientificreports/ Current views considered IS as a multifactorial disease with predisposing genetic factors 39 . Many factors contributing to the occurrence of IS were described, such as mechanical, hormonal, metabolic, neuromuscular, growth, and genetic abnormalities. Amongst these, some factors may be epiphenomena rather than the cause itself. Other factors may contribute to curve progression, rather than curve initiation 40 .
By studying two sides of the paravertebral muscles of IS patients, we aimed to identify the molecular background for AIS-associated paravertebral muscle changes.
The development and growth of the skeleton are controlled by parathyroid, thyroid hormones, and growth hormones, but also by estrogens 41 . It has been reported that factors associated with puberty and sexual dimorphism play an important role in the pathogenesis of idiopathic scoliosis. However, results demonstrating this relationship are ambiguous. Kulis et al. observed that blood estrogen levels in female patients with idiopathic scoliosis are significantly lower compared to unaffected individuals 42 . Conversely, Raczkowski and colleagues found that the level of circulating estrogens in the blood did not differ significantly between females in a control group and the group of patients with IS 43 . Although the role of estrogens in IS has not been fully uncovered, it has been suggested that estrogen-mediated signal transduction may be dysregulated in affected cells 27,43 .
In previous studies, we observed asymmetric expression of ESR2 in deep paravertebral muscles on the concave and convex side of the curvature. The relative expression of the ESR2 gene correlated with Cobb angle values and the risk factor of progression 26 . However, the mechanism of estrogen action in IS has yet to be established. In recent years, the role of epigenetic factors in the etiopathogenesis of IS has been increasingly investigated 32,44 . ESR2 expression may be regulated by methylation of CpG islands within the 0N and 0K promoters and their corresponding exons, which generate transcript variants that differ in their 5′UTR (5′ untranslated) regions 37,45 . However, there has yet to be an assessment of DNA methylation at ESR2 regulatory regions in the deep paravertebral muscles of patients with idiopathic scoliosis.
We observed that methylation levels within ESR2 promoter 0N, but not exon 0N, were significantly higher on the concave side of the curvature compared to the convex side. The difference in methylation levels in deep paravertebral muscles may be associated with their heterogeneous histological architecture 46 . Muscle tissue on the concave side was characterized by a higher content of fibrous elements. Moreover, the ratio of slow-twitch to fast-twitch fibers is shifted, with a higher prevalence of fast-twitch fibers on the concave side of the curvature 47 . www.nature.com/scientificreports/ www.nature.com/scientificreports/ What is more, higher electromyographic activity in paravertebral muscles on the convex side than on the concave side was reported. This phenomenon was studied in relation with potential treatment possibilities and predictive value 48,49 . According to Weiss, this asymmetry can be reduced with specific exercises 50 .
However, it is difficult to distinguish the causative factors from their effects. The difference in methylation may be the cause of the asymmetry in muscles, which may have some contribution to the etiology of IS. Conversely, differences in methylation levels could be a consequence of the muscles being exposed to different conditions on either side of the curvature due to asymmetric loading or other unknown conditions. Taking into consideration that the curvature size did not correlate with the methylation level, the asymmetrical loading is likely not the www.nature.com/scientificreports/ reason for the difference in methylation level, because the asymmetry of loading transmission increases with the curvature size. Postural control deficits in patients with IS was described by Sahlstrand et al. and confirmed in more recent publications and seems to be a repeatable observation 51,52 . When discussing the postural imbalance in IS, neurophysiological function and structure of muscles, the theory of a functional tethering of the spinal cord as a possible background factor for IS needs to be mentioned. Chu et al. described the relative shortening and functional tethering of spinal cord in IS patients. They found an increased prevalence of abnormal somatosensory evoked potentials and the low-lying cerebellar tonsil in severe IS and compared with mild to moderate IS. Similar observations concerning the association of IS with the severe curve with tonsillar ectopia and abnormal somatosensory function were described by Cheng et al. 53 . Deng et al. postulated the cord-vertebral length ratio as a significant predictor for curve progression in IS 54 . Taking into consideration the impact of the functional spinal cord tethering on the signal transduction this phenomenon may affect both function and histological structure of the paraspinal muscles.
It is widely accepted that epigenetic modifications can impact gene expression 55 . Consistent with this, we found an association between the expression of ESR2 and promoter 0N methylation. Specifically, ESR2 expression negatively correlated with methylation level, which is consistent with the phenomenon of methylation-dependent gene silencing. These results provide further support for the hypothesis that 0N promoter methylation modulates ESR2 gene expression levels in the deep paravertebral muscle tissues in the case of patients with IS.
The phenotypic heterogeneity of IS makes it challenging to form an efficient strategy to prevent the disease or its progression 56 . The phenotype of IS differs vastly between patients. Our division of IS phenotypes according to curve severity was based on clinical studies [57][58][59] . According to clinical practice, IS between 50º and 70º requires surgeon intervention due to the risk of further progression. The aim of surgery in such cases is to avoid future complications due to possible progression [57][58][59] . On the other hand, in more severe scoliosis, significant health problems associated with lung function, cardiac function, and back pains can occur 2,4,58 .
Understanding the etiopathogenetic background of IS susceptibility as well as factors enhancing curvature progression has crucial clinical implications. Identification of the patients who are most at risk of curve progression allows clinicians to modify the treatment approach. Currently, no defined threshold marks when severe scoliosis produces a significant impact on a patient's health. Most studies concerning surgical treatment of scoliosis classify severe curvature as a Cobb angle exceeding 70°6 0-62 . Thus, we applied this value to categorize study groups.
To our knowledge, this study is the first to evaluate DNA methylation in muscle tissues affected by IS. Previous studies did not consider locally acting factors, instead focusing on global markers [32][33][34][35][36] . We found just one study investigating epigenetic changes in muscle tissue in IS patients. Jiang et al. evaluated the involvement of non-coding RNA in IS-affected tissue 63 . It is difficult to compare our results directly since we investigated different epigenetic mechanisms. Nevertheless, both studies support the importance of evaluating local changes in the affected region of the body to find causative agents of disease.
Several attempts have been made to evaluate the importance of DNA methylation in IS etiopathogenesis. All these studies provide essential information about the role of DNA methylation in IS etiopathogenesis. However, the chosen tissue (peripheral blood) limited these findings to the evaluation of globally acting factors.
A strength of our study is that we performed a comprehensive evaluation of pathology at the muscular origin of the disease. Thus, our data add valuable insights into the DNA methylation and its biological action in the background of IS. Another key feature of the study is the analysis of IS phenotypes. Although we did not identify an association between methylation and IS severity, our results support the theory that factors associated with occurrence and progression of IS differ considerably.
Finally, we acknowledge that our study also had some limitations. Firstly, our investigation was limited by the sample size. However, the number of patients is comparable to other published research investigating DNA methylation in IS [32][33][34][35][36] . The study was also limited by the lack of tissue samples obtained from healthy controls. We decided not to acquire control samples from individuals who underwent surgery due to degenerative spine disease. This decision was taken because such patients are mostly elderly, and chronic degenerative disease of the spine may cause the muscles atrophy, or induce unknown methylation changes. Similarly, it would be beneficial to compare our results in the case of patients with curves between 20 and 30° to evaluate the impact of ESR2 methylation on IS progression. However, such patients do not need surgical treatment, and there is no possibility to obtain muscle samples.

Conclusions
In conclusion, our results demonstrate that the methylation pattern of CpG sites in the regulatory regions of the ESR2 gene in the deep paravertebral muscle tissue is associated with the occurrence but not with the severity of idiopathic scoliosis. The tissue samples obtained from patients who underwent surgery due to IS showed substantial differences in methylation levels at ESR2. Altogether, the presented findings suggest that differences in methylation levels at the concave compared to the convex side of the spinal curvature may be associated with IS etiopathogenesis. The relationship between the promoter and exon 0N methylation levels and ESR2 gene expression indicates that this type of epigenetic modification may affect the tissue-dependent response to estrogens.

Methods
The study was approved by the Institutional Review Board of the Poznan University of Medical Sciences (No 546/17 and 741/19). Informed consent was obtained from all the patients or their parents/legal guardians in the case of under-aged participants. All methods were carried out in accordance with the approved guidelines.
Patients. Twenty-nine girls with severe IS were included in the study. All patients underwent posterior spinal surgery in one hospital in a Central European country (Poland) from January 2017 until December 2019. Patients were subject to a clinical, radiological, and molecular examination. All patients had undergone standing posteroanterior X-rays before surgery. The curve pattern (number and localization of the curvatures), Cobb angle (angle of curvature size) 64 , and Risser sign (the radiological sign of skeletal maturity) 65 were measured by an experienced spine surgeon. The inclusion criteria were as follows: (1) clinically and radiologically confirmed IS diagnosis, (2) no coexisting orthopedic, genetic or neurological disorders, (3) primary thoracic spinal curvature (4) surgical treatment due to IS. The patients were divided into groups according to disease severity. The first group consisted of ten patients with a moderate form of idiopathic scoliosis with curvature ranging from 50° to 70° and a Risser sign of ≥ 3 or age ≥ 15 years old. The second group consisted of nineteen patients with a very progressive form of idiopathic scoliosis with larger curvatures exceeding 70°, regardless of Risser sign or age.
Tissue samples. During the surgery, two muscle tissue fragments were obtained from each patient. These were from the deep paravertebral muscles (m. longissimus) on the (1) convex and (2) concave side of the curvature (1 cm 3 each). Samples were stored in sterile tubes containing 5 ml nucleic acid stabilizing solution (Novazym).
Genomic DNA methylation analysis. Genomic DNA isolation. Genomic DNA was extracted using Quick-DNA Miniprep Plus Kit (Zymo Research) according to the manufacturer's protocol with modifications. Firstly, tissue samples were ground in liquid nitrogen with a mortar and pestle. Then, 25 mg were incubated overnight at 55 °C with proteinase K. Next, the lysate was centrifuged to remove insoluble debris. At this point, isolation was continued following the steps outlined in the protocol. The DNA quantity and purity were assessed using a spectrophotometer (NanoPhotometer NP80, Implen) by measuring absorbance at A = 260, A = 230 nm, and ratios A = 260/230, A = 260/280. All analyzed samples met the criteria for purity (both ratios ranged from 1.9 to 2.0). DNA integrity was evaluated using a standard 1% agarose gel (Lab Empire) electrophoretic separation in the presence of ethidium bromide (50 ng/ml, Merck).
Bisulfite conversion. One microgram of genomic DNA was bisulfite converted using EZ DNA Methylation Kit (Zymo Research) following the manufacturer's protocol and eluted with 10 µl of M-Elution Buffer.
Polymerase chain reaction and pyrosequencing analysis. Pyrosequencing reactions were preceded by polymerase chain reaction (PCR) with bisulfite converted DNA as the template. Specific primers for both reactions were designed using PyroMark Assay Design 2.0 software (Qiagen) and synthesized by Genomed. The input DNA sequences corresponded to the promoter 0N and exon 0N regions of ESR2 gene reference sequences deposited in the Nucleotide Database of National Center for Biotechnology Information (https ://www.ncbi.nlm.nih.gov; GenBank No.: NG_011535.1). Sequencing, forward, and reverse primers are presented in Table 1. PCR was performed in a total volume of 10 µl using ZymoTaq PreMix (Zymo Research) designed for the amplification of bisulfite-treated DNA ( Table 2). 2 µl of the products were analyzed using standard 2% agarose gel (Lab Empire) electrophoretic separation, and compared to Nova 100 mass marker (Novazym) in the presence of ethidium bromide (50 ng/ml, Merck).
Pyrosequencing analysis was performed using the PyroMark Q48 instrument (Qiagen). CpG assays were designed using Pyromark Q48 Autoprep 2.4.2 software (Qiagen). We analyzed 16 and 19 CpG sites for promoter and exon 0N, respectively. Each reaction contained an internal control that allowed to assess whether the sodium bisulfite treatment was successful. Methylation levels were quantified using Pyromark Q48 Autoprep 2.4.2 software (Qiagen) and determined as a percentage ratio of methylated to non-methylated dinucleotides. Analysis of ESR2 mRNA expression level. Total RNA isolation. Total cellular RNA was extracted using RNA Isolation Reagent (GenoPlast Biochemicals) and Direct-zol RNA Miniprep Kit (Zymo Research) according to the manufacturer's protocol with modifications. Firstly, 25 mg of powdered muscle tissue was transferred to RNA Isolation Reagent and incubated at room temperature for 5 min (min). Next, 200 µl of chloroform was added, and samples were shaken vigorously and incubated at room temperature for 3 min. Then, samples were centrifuged (12,000×g for 15 min at 4 °C), and the upper aqueous phase was subsequently transferred to an equal volume of absolute ethanol. At this point, isolation was continued following the steps outlined in the protocol. The RNA sample quantity and purity were assessed using a spectrophotometer (NanoPhotometer NP80, Implen). All analyzed samples met the criteria for purity. RNA integrity was evaluated by 18S and 28S ribosomal RNA bands presence observation, using a standard denaturating 1% agarose gel (Lab Empire) electrophoretic separation in the presence of ethidium bromide (50 ng/ml, Merck).
Reverse transcription reaction. cDNA (complementary DNA) synthesis was performed using Expand Reverse Transcriptase (Roche) according to the manufacturer's protocol with modifications described below. The total volume of reaction was 10 µl. In the first step, the mixture containing 500 ng of total RNA, DNase, RNase, Quantitative polymerase chain reaction. Quantitative analysis of ESR2 mRNA expression was evaluated using a hydrolysis probe PrimePCR (qHsaCEP0052206, BioRad). The hypoxanthine-guanine phosphoribosyltransferase (HPRT) gene was used as a reference gene (RealTime ready HPRT, 102079, Roche). The total volume of 20 µl reaction mixture contained 5 µl of cDNA, 1 × LightCycler FastStart TaqMan Probe Master (Roche), 1 × ESR2 PrimePCR or 1 × RealTime ready HPRT, and nuclease-free water. qPCR reactions were performed using the LightCycler 2.0 carousel glass capillary-based system (Roche). The thermal profile of the qPCR reaction was as follows: pre-incubation step at 95 °C for 10 min, 45 quantification cycles (denaturation at 95 °C for 10 s, annealing/extension step at 60 °C for 30 s, and the final step at 72 °C for 1 s (fluorescence level acquisition mode) and the final cooling to 40 °C for 30 s. Each sample was analyzed in duplicate with independently synthesized cDNA. The quantitative PCR results were assembled using the LightCycler Data Analysis (LCDA) Software version 5.0.0.38, and the obtained fluorescence measurement results were normalized to standard curves. In each sample, ESR2 expression levels were compared to the reference gene expression level to obtain Cr value (concentration ratio), which corresponded to the relative ESR2 expression level.  ACT TAA CCA TAA ACC CCT  TCT TCC TTT  27  58,9  37,0  + 37  − 11,387   SEQ  ATA TTT TTA GGT TTT ATT  TTA GAT  24