Introduction

The human spinal column is a dynamic, segmented, bony and cartilaginous structure that is essential for integrating the brain and nervous system with the axial skeleton while simultaneously providing flexibility in three dimensions1. Idiopathic scoliosis is the most common developmental disorder of the spine, typically appearing during the adolescent growth spurt. Adolescent idiopathic scoliosis (AIS) is reported in all major ancestral groups, with a population prevalence of 1.5-3%2,3. Children with AIS usually present with a characteristic right-thoracic major curve pattern and a compensatory lumbar curve. Major thoracolumbar and lumbar curves are less frequent4. The three-dimensional nature of the deformity results in torsion in the spine that is most significant at the apex of the major curve, and changes in the structures of the vertebrae and ribs may develop as the curve worsens or progresses4. Children with thoracic curves, with larger curves at first presentation, and/or with greater remaining growth potential are at increased risk of progression, but this risk decreases sharply after skeletal maturity4. Sex is a recognized risk factor for AIS, with girls having at least a five-fold greater risk of progressive deformity requiring treatment compared to boys5. This well-documented sexual dimorphism has prompted speculation that levels of circulating endocrine hormones, particularly estrogen, are important exposures in AIS susceptibility6.

The genetic architecture of human AIS is complex, and underlying disease mechanisms remain uncertain. Heritability studies of Northern European7,8, North American9,10, and South Asian11 ancestral groups suggest that disease risk is multifactorial, caused by genetic and environmental contributions2,12. Accordingly, population-based genome-wide association studies (GWAS) in multiple ancestral groups have identified several AIS-associated susceptibility loci, mostly within non-coding genomic regions12. In particular, multiple GWAS have implicated noncoding regions near the LBX113, ADGRG6 (also known as GPR126)14, and BNC215 genes. An association with alleles in an enhancer distal to PAX1, encoding the transcription factor paired box 1, was primarily driven by females, suggesting that it contributes to the sexual dimorphism observed in AIS16. Subsequent meta-analysis of combined AIS GWAS identified additional susceptibility loci. These included variants in an intron of SOX6, a transcription factor, that along with PAX1, is important in early spinal column formation17. Furthermore, gene enrichment analyses found significant correlation of AIS-associated loci with biological pathways involving cartilage and connective tissue development18. A more recent GWAS in a Japanese population identified fourteen additional AIS loci that are candidates for further evaluation19. In separate studies, genome sequencing in AIS cases and families identified enrichment of rare mutations in the COL11A220 and HSPG221 genes, encoding components of the cartilage extracellular matrix (ECM). Hence variation affecting cartilage and connective tissue ECM genes is an emerging theme in the heterogeneous genetic architecture of AIS.

Pre-clinical animal models are essential tools for accelerating mechanistic understanding of AIS and for therapeutic testing12. In zebrafish, several genetic mutants with larval or later-onset spinal deformity have been described, including ptk722,23, c21orf5924, ccdc4025, ccdc15126, dyx1c1, and kif627. In mouse, conditional deletion of Adgrg6 in skeletal cartilage (using Col2a1-Cre) produces a progressive scoliosis of the thoracic spine during postnatal development that is marked by herniations within the cartilaginous endplates of involved vertebrae. Progressive scoliosis, albeit to a lesser extent, was also observed when Adgrg6 was deleted from committed chondrocytes (using ATC-Cre)2830. These studies demonstrate that cartilage and possibly other osteochondroprogenitor cells contribute to the scoliosis phenotype in these models29. Taken together, genetic and functional studies in mouse, although limited, support the hypothesis that deficiencies in biogenesis and/or homeostasis of cartilage, intervertebral disc (IVD), and dense connective tissues, undermine the maintenance of proper spinal alignment during the adolescent growth spurt12.

The combined contribution of reported AIS-associated variants is broadly estimated to account for less than 10% of the overall genetic risk of the disease19. To address this knowledge gap, we sought to define novel loci associated with AIS susceptibility in genes encoding proteins of the ECM (i.e. the “matrisome”31,32). Here we identify new genetic associations with AIS. Further, our functional assessments support a new disease model wherein AIS-associated genetic variation and estrogen signaling perturb a Pax1-Col11a1-Mmp3 axis in cartilage GPCs.

Results

Nonsynonymous variants in matrisome genes are associated with increased risk of AIS

The “matrisome” has been defined as “all genes encoding structural ECM components and those encoding proteins that may interact with or remodel the ECM”33. Proteins comprising the global ECM as currently defined have been identified by both experimental and bio-informatic methods31. We assembled 1027 matrisome genes as previously identified34, including 274 core-matrisome and 753 matrisome-associated genes (N=1,027 total). For the genes encoding these 1,027 proteins, we identified all nonsynonymous common variants (MAF > 0.01) queried by the Illumina HumanCoreExome-24v1.0 beadchip and determined their genotypes in a discovery cohort of 1,358 cases and 12,507 controls (Table 1). After applying multiple quality control measures (see Methods section), we retained 2,008 variants in 597 matrisome genes for association testing (Supplemental Table 1). For controls, we utilized 12,507 non-AMD GRU (non-age-related macular degeneration general research use) subjects of European ancestry downloaded from dbGaP (https://www.ncbi.nlm.nih.gov/gap/). This sample size was estimated to provide at least 80% power to detect significant associations at the matrisome-wide level (α≤2.5e-05), for alleles with population frequency ≥.05 and OR ≥1.5 (Figure 1 – figure supplement 1a). Two nonsynonymous variants, in COL11A1 (rs3753841; NM_080629_c.4004C>T; p.(Pro1335Leu); odds ratio (OR)=1.236 [95% CI=1.134-1.347], P=1.17E−06) and MMP14 (rs1042704; NM_004995_c.817G>A; p.(Asp273Asn); OR=1.239 [95% CI=1.125-1.363], P=1.89E−05) were significantly associated with AIS (Figure 1a). Given the sexual dimorphism in AIS and our prior observation of a female-predominant disease locus16, we tested the 2,008 variants separately in females (N=1,157 cases and 7,138 controls). In females, the association with rs3753841 remained statistically significant, whereas rs1042704, near MMP14, was not associated with AIS in females (Figure 1 – figure supplement 1b). Our study was not sufficiently powered to test males separately.

Matrisome-wide association study.

a. Manhattan plot showing −log10 P-values (y axis) versus chromosomal position (x axis) for the 2,008 common coding variants tested in the discovery study USA (TX). The horizontal line represents the threshold for significance level (P-value<2.5×10-5) after Bonferroni multiple testing correction. b. Tests of association for SNPs rs3753841 and rs1042704 in discovery and independent replication cohorts. RAF reference allele frequency; OR odds ratio; CI-confidence interval.

Participating studies.

NHW: Non-Hispanic White. EAS: East Asian.

To validate these results, we sought to replicate the associations of rs3753841 and rs1042704 in four independent AIS case-control cohorts, from North America, Europe, and eastern Asia, representing multiple ethnicities (total N = 9,161 AIS cases, 80,731 healthy controls, Table 1). Genotypes for both variants were extracted from these datasets and tested for association by meta analysis together with the discovery cohort (see Methods). Meta-analysis of all cohorts together increased the evidence for association of both variants with AIS risk (Figure 1b). While a similar effect size was noted for rs1042704 in Japanese and Han Chinese cohorts, the results were less significant, likely due to lower minor allele frequencies (East Asian MAF = 0.02 compared to total non-Asian cohort MAF = 0.20) in these populations (Figure 1 – figure supplement 1c). Plotting recombination across both regions suggested that these signals were likely confined to blocks of linkage disequilibrium within the COL11A1 and MMP14 genes, respectively (Figure 1 – figure supplement 1d, e).

Although low-resolution structures currently available for collagen triple helices are not useful for modeling the stability change of individual mutations, we noted that the AIS-associated variant P1335L occurs at the third position of a Gly-X-Y repeat and consequently could be structurally important in promoting stability of the triple helix, particularly if it is hydroxylated. Rare dominant mutations in COL11A1, often disrupting a Gly-X-Y sequence, can cause Marshall (MRSHS) (OMIM# 154780) or Stickler syndromes (STL2) (OMIM# 604841) marked variously by facial anomalies, sensineural hearing loss, short stature, spondyloepiphyseal dysplasia, eye anomalies, ectodermal features, and scoliosis. Notably, our AIS cohort and particularly individuals carrying the rs3753841 risk allele were negative for co-morbidities or obvious features of Marshall or Stickler syndromes. Thus, AIS is independently associated with the COL11A1P1335L predicted variant.

COL11A1 is expressed in adolescent spinal tissues

We elected to further characterize COL11A1 in postnatal spine development. COL11A1 encodes one of three alpha chains of type XI collagen, a member of the fibrillar collagen subgroup and regulator of nucleation and initial fibril assembly, particularly in cartilage35. Spinal deformity is well-described in Col11a1-deficient (cho/cho) embryos36. In mouse tendon, Col11a1 mRNA is abundant during development but barely detectable at three months of age37. We analyzed RNA seq datasets derived from adolescent human spinal tissues38, finding that COL11A1 was upregulated in cartilage relative to bone and muscle. In cartilage, PAX1 and COL11A2, showed the strongest expression levels relative to other published human AIS-associated genes1416,18,20,21,39 (Figure 2). In all, the majority of AIS-associated genes showed the strongest expression levels in cartilage relative to other adolescent spinal tissues.

Col11a1 and Mmp14 expression in spine.

a) A heatmap of transcript per million (TPM) values of COL11A1, MMP14, and other published genes associated with AIS. The average TPM value of matrisome genes is represented as MATRISOME.

b) Localization of collagen type XI alpha chain 1 in P0.5 mouse IVD. Immunohistochemistry (IHC) and immunofluorescence (IF) staining in IVD at day P0.5. Negative antibody controls in IVD are shown at left; antibody-positive IHC and IF are shown at right.

c) Localization of collagen type XI alpha chain 1 in P28 mouse IVD. Negative antibody IHC control shown at left; antibody positive IHC shown at right. Enlarged, rotated view of white boxed area shows a biphasic staining pattern. CEP – cartilage endplate; GP – growth plate.

We next sought to characterize Col11a1 expression in spines of postnatal mice. To detect Col11a1 protein (collagen α1(XI)), we performed immunohistochemistry (IHC) and immunofluorescence (IF) microscopy using a collagen α1(XI) reactive antibody40 in newborn (P0.5) and adolescent (P28) mice. In spines of P0.5 mice, strong staining was observed in the nucleus pulposus (NP) and in surrounding annulus fibrosus (AF) (Figure 2b). In P28 mice, the compartments of the IVD were more distinct, and strong collagen α1(XI) staining was observed in each (Figure 2c). In regions of the cartilage endplate (CEP)-vertebral bone transition, collagen α1(XI) was detected in columnar chondrocytes, particularly in the hypertrophic zone adjacent to condensing bone (Figure 2c). We also examined collagen α1(XI) expression in ribs, as these structures are also involved in the scoliotic deformity4. In P28 rib growth plates, as in spine, a biphasic pattern was observed in which collagen α1(XI) reactivity was most pronounced around cells of the presumed resting and pre-hypertrophic/hypertrophic zones (Figure 2 – figure supplement a,b). These data show that in mouse, collagen α1(XI) is detectable in all compartments of young postnatal IVD and is particularly abundant in the chondro-osseous junction region of IVD and vertebral growth plate.

Col11a1 is positively regulated by Pax1 in mouse spine and tail

We previously identified AIS-associated variants within a putative enhancer of PAX1 encoding the transcription factor Paired Box 116,18. Pax1 is a well-described marker of condensing sclerotomal cells as they form segments that will eventually become the IVD and vertebrae of the spine4143. We generated Pax1 knock-out mice (Pax1−/−) using CRISPR-Cas9 mutagenesis and validated them using sequencing and Southern blot (Figure 3-figure supplement 3a-d). Homozygous Pax1−/− mice were viable and developed kinks in the tail, similar to other Pax1 deficient mice44. We next compared the expression of collagen α1(XI) protein in IVD and condensing bone of wildtype and Pax1−/− mice by performing IF staining in P28 spines (Figure 3a). In wildtype IVD, strong overlapping expression of collagen α1(XI) and Pax1 cells was observed, mostly within the CEP and chondro-osseous interface (Figure 3a). The IVD in Pax1−/−mice was highly disorganized, without discernable NP, AF, and CEP structures as has been reported45. Cells within the IVD were negative for Pax1 staining as expected, and collagen α1(XI) staining was also clearly reduced in CEP and the chondro-osseous vertebral border. To test the effect of Pax1 on expression of Col11a1 and other AIS-associated genes during embryonic development, RNA was isolated from vertebral tissue dissected from the tails of E12.5 wild-type and Pax1−/− mice and subjected to bulk RNA-seq and quantitative real-time PCR (qRT-PCR). Gene-set enrichment analysis of RNA-seq was most significant for the gene ontology term “extracellular matrix” (Figure 3b). Interestingly, the AIS-associated gene Adgrg6 was amongst the most significantly dysregulated genes in the RNA-seq analysis (Figure 3c). By qRT-PCR analysis, expression of Col11a1, Adgrg6, and Sox6 were significantly reduced in female and male Pax1−/− mice compared to wild-type mice (Figure 3d-g). These data support a mechanistic link between Pax1 and Col11a1, and the AIS-associated genes Gpr126 and Sox6, in affected tissue of the developing tail.

Assessing Pax1 regulation of Col11a1.

a) IF staining of P28 IVD using Pax1 (green) and Col11a1 specific (red) antibodies with DAPI nuclear counterstain. Arrow denotes cells at the CEP-condensing bone interface of IVD that show positive staining for both Pax1 and Col11a1.

b) Gene ontology (GO) analysis of differentially expressed genes in E12.5 tail WT and Pax1-null mice.

c) Heatmap of differentially expressed genes (P-value<0.0001) in E12.5 tails of WT and Pax1-null mice

d-g) Gene expression levels dissected from E12.5 mouse tail from wild type (WT) and Pax1-null (KO) mice as determined by qRT-PCR. Each value represents the ratio of each gene expression to that of β-actin, and values are mean ± standard deviation. The expression value of WT female group was arbitrarily set at 1.0. Each dot represents one embryo and statistical differences were determined using a two-sided unpaired t-test (P*<0.05, P**<0.01, P***<0.001).

Col11a1 regulates Mmp3 expression growth plate chondrocytes

COL11A1 has been linked with ECM remodeling and invasiveness in some cancers46. In solid tumors, COL11A1 has been shown to alter ECM remodeling by enhancing MMP3 expression in response to TGFý46. MMP3 encodes matrix metalloproteinase 3, also known as stromolysin, an enzyme implicated in matrix degradation and remodeling in connective tissues47. We confirmed strong MMP3 mRNA expression, relative to COL11A1, in human spinal cartilage and bone, but minimal expression in spinal muscle (Figure 4 – figure supplement 4a). We next generated Col11a1-deficient chondrocytic cell cultures from rib growth plate chondrocytes (GPCs) of P0.5 Col11a1fl/fl mice40 by treating with Cre-expressing adenoviruses. After confirming Col11a1 excision (Figure 4a), we compared Mmp3 expression in these cells to cells treated with GFP expressing adenoviruses lacking Cre activity. We found that Mmp3 expression was significantly increased in GPCs where Col11a1 mRNA expression was downregulated by about 70% (Figure 4b). Furthermore, Western blotting in these cells demonstrated a ∼2-5-fold increase in pro-, secreted, and active forms of Mmp3 protein when collagen α1(XI) was reduced. The proteolytic processing per se of precursor Mmp3 into active forms48 did not appear to be affected by Col11a1 expression (Figure 4c). These results suggest that Mmp3 expression is negatively regulated in Col11a1 in mouse GPCs.

Col11a1 regulation of Mmp3 expression in cartilage.

a) PCR assay of Col11a1 excision in Col11a1fl/fl cultured rib cartilage cells.

b) Gene expression levels from Col11a1fl/fl cultured rib cartilage cells transduced with GFP (Ad5-GFP, left) or Cre expressing adenovirus (Ad5-cre, right) as determined by qRT PCR. Each value represents the ratio of each gene expression to that of GAPDH, and values are mean ± standard deviation. The expression value of control Ad5-GFP results (N≥ 4 experiments) was arbitrarily set at 1.0. Statistical differences were determined using a two-sided unpaired t test (P*<0.05).

c) Representative Western blot of cultured cartilage cells after Ad5-GFP or Ad5-cre transduction. Protein size ladder is shown in lane 1. Quantification of bands detected by Western blotting, where scramble results were set to 1.0 from N=4 independent experiments is shown at right.

d) Gene expression levels from dissected Col11a1fl/fl ATC rib cartilage, analyzed as described in a). Results shown are from N=3 mice.

To test whether Col11a1 affects Mmp3 levels in vivo, we bred Col11a1fl/fl female mice with Col11a1fl/fl ATC males carrying the Acan enhancer-driven, doxycycline-inducible Cre (ATC) transgene49. ATC has been shown to harbor Cre-mediated recombination activity in most differentiated chondrocytes and in nucleus pulposus within two days of treating pregnant mothers with doxycycline starting at E15.549. ATC activity was confirmed by crossing this line to the R26td[Tomato] reporter that ubiquitously expresses the fluorescent gene Tomato after Cre recombination. Strong Cre activity was seen in P0 pups of mothers treated with doxycycline at E15.5 in the NP, CEP, and annulus fibrosus (AF) of the IVD and in chondrocytes of the growth plates (Figure 4 – figure supplement 4b). Pregnant Col11a1fl/fl females were treated with doxycycline water from E15.5 to induce Cre expression in differentiated chondrocytes. Excision of Col11a1 was confirmed in DNA from costal cartilage of Col11a1fl/flATC cre-positive offspring (Figure 4 – figure supplement 4c). Consistent with results obtained by in vitro excision of Col11a1, cartilage from mice deficient in Col11a1 showed ∼4-fold upregulation of Mmp3 mRNA expression relative to Col11a1fl/flmice (Figure 4d).

AIS-associated variant in COL11A1 perturbs its regulation of MMP3

COL11A1 missense mutations have been shown to evoke transcriptional changes in ECM genes in cancer cells50. We therefore tested whether the COL11A1P1335L sequence variant alters its regulation of Mmp3 in cartilage chondrocytes. For this, SV40-immortalized cell lines were established from Col11a1fl/fl mouse rib cartilage and transduced with lentiviral vectors expressing green fluorescent protein (GFP) and COL11A1wt, COL11A1P1335L, or vector alone. After transduction, GFP-positive cells were grown to 50% confluence and treated with Cre-expressing adenovirus to remove endogenous mouse Col11a1 (Figure 5a,c). Using a human-specific COL11A1 qRT-PCR assay, we detected overexpression of COL11A1wt and COL11A1P1335L compared to untransduced cells regardless of Cre expression (Figure 5a). Western blotting with an antibody directed against the HA epitope tag confirmed also confirmed overexpression of human COL11A1 (Figure 5c). Endogenous Mmp3 upregulation was evident by both qRT-PCR and Western blotting in untransduced cells treated with Ad5 Cre, as expected, and the reverse was observed in cells overexpressing human wildtype COL11A1, consistent with the negative regulation we previously observed (Figure 5 b,c). However, the COL11A1P1335L mutant failed to downregulate Mmp3 expression despite being overexpressed (Figure 5b,c). Thus, regulation of Mmp3 is perturbed in the presence of the COL11A1P1335L variant in these cells.

Col11a1P1335L regulation of Mmp3 expression in lentiviral transduced mouse GPCs.

a) qRT-PCR of human COL11A1 and endogenous mouse Mmp3. “SV40” refers to immortalized mouse chondrocytes transduced with the lentiviral vector only. Significant fold changes (P ≤0.5) relative to vector-only transfected cells are shown by *.

b) Western blot corresponding to experiments shown in (a) using antibodies to detect Col11a1, HA epitope tag representing human Col11a1 only, endogenous Mmp3, and GAPDH. Quantification at right shows Mmp3 expression in each experiment relative to GAPDH.

Esr2-Col11a1-Mmp3 form a signaling axis in cartilage cells

The expression of Col11a1, and of other ECM genes, is known to be estrogen-responsive in certain tissues, such as ovarian follicular cells51. Because of the suspected role of endocrine hormones in AIS, we investigated whether Col11a1 expression was responsive to estrogen treatment and/or levels of estrogen receptor following siRNA-mediated knockdown experiments in cultured cartilage GPCs. We first validated that Mmp3 levels were significantly increased after Col11a1 knockdown (Figure 6a). Estrogen receptor 2 (Esr2), but not estrogen receptor alpha (Esr1), was detected in mouse GPCs by qRT-PCR (data not shown). We next tested the consequences of Col11a1 and Esr2 siRNA-mediated knockdown on gene expression in cartilage GPCs. After Esr2 knockdown, Col11a1 as well as Pax1 were significantly upregulated compared to scramble control, while Mmp3 expression was significantly downregulated (Figure 6b). Pax1 expression was also upregulated by Col11a1 knockdown alone. These results suggest a negative feedback loop between Pax1 and Col11a1 in these cells. Simultaneous knockdown of Col11a1 and Esr2 expression reduced Mmp3 expression to normal levels, supporting a possible interaction between Col11a1 and estrogen signaling in regulating Mmp3. Treating GPCs with tamoxifen, an estrogen receptor modulator, also upregulated Col11a1 expression to similar levels as observed after Esr2 knockdown, compared to cells treated with DMSO carrier (Figure 6 figure supplement 6a). These results further support estrogen signaling in regulating Col11a1 expression. In cultured rat CEP cells, Esr2 mRNA was downregulated, and Mmp3 mRNA was upregulated after Col11a1 knockdown, as observed in mouse GPCs (Figure 6c, Figure 6 – figure supplement 6b). However, Esr2 knockdown did not significantly impact Col11a1 or Mmp3 expression in these cells (Figure 6c). Hence, we conclude that in mouse GPCs, expression of the ECM-clearing enzyme Mmp3 gene is normally repressed by collagen α1(XI), but this regulation is disrupted by AIS-associated variation in COL11A1. Further, we find that estrogen signaling also dysregulates the Col11a1-Mmp3 interaction.

Effects of estrogen receptor beta on Col11a1 Mmp3 signaling axis.

a) Representative Western blot of cultured cartilage cells after scramble or Col11a1-specific siRNA knockdown. Protein size ladder is shown in lane 1. Quantification of bands detected by Western blotting, where scramble results were set to 1.0 from N=4 independent experiments is shown at right.

b) Gene expression levels of Col11a1 (left), Mmp3 (middle), and Esr2 (right) mRNA in cultured rib cartilage cells showing fold change relative to the scramble control. ddk = double Col11a1-Esr2-specific siRNA knockdowns. Results are representative of n≥3 experiments.

c) Gene expression levels from rat CEP cells, as described in (b).

Discussion

Adolescent idiopathic scoliosis has been described in the medical literature for centuries, yet its underlying etiology has remained enigmatic52. Given that AIS originates in children who appear to be otherwise healthy, even its tissue of origin has been difficult to discern, and long debated12. The advent of powerful genotyping and sequencing methods in the last two decades has led to breakthrough discoveries of genetic loci associated with AIS, most in non-coding regions of the genome that are difficult to interpret biologically12. Aggregating these results, however, provided supportive evidence that pathways of cartilage and connective tissue ECM development are relevant in AIS etiology12, 18. Here, in the largest multi-ethnic human cohort studied to date, we elected to test the hypothesis that alterations in ECM proteins themselves contribute to AIS susceptibility. This approach yielded most significant evidence for a common protein-altering variant in the COL11A1 gene encoding collagen α1(XI), a minor yet critical component of cartilaginous ECM. Moreover, our studies define a COL11A1-mediated disease pathway (Figure 7) and point to the chondro-osseous junction of IVD and vertebrae spine as a relevant cellular compartment in AIS etiology.

Cartoon depiction of a collagen XI-mediated signaling axis in cartilage. Collagen XI is held in the pericellular space by integrins and DDR2. COL11A1, under the regulation of ESR2 and PAX1, signals through unknown mechanisms and inhibits MMP3 transcription.

The results of this study together with the previous observation of COL11A2 rare variant enrichment in AIS supports a role for the collagen α1(XI) heterotrimer itself in its pathogenesis20. Collagen type XI, composed of three chains encoded by the COL11A1, COL11A2, and COL2A1 genes (OMIM #s 120280,120290, 120140, respectively), is a minor component of collagen type II fibrils that are abundant in cartilage. Collagen type XI is also broadly expressed in testis, trachea, tendons, trabecular bone, skeletal muscle, placenta, lung, brain neuroepithelium, the vitreous of the eye, and intervertebral discs53. In the pericellular space, collagen α1(XI) initiates fibrillogenesis with collagen type II fibrils, maintaining regular spacing and diameter of the collagen fibrils, while organizing the pericellular surface by interaction with cartilage proteoglycans54,55. Purified human collagen type XI, when added back to chondrocytes in in vitro culture, stimulates chondrogenesis while inhibiting hypertrophy, as measured by histologic staining, proliferation assays, and relative expression of chondrogenic early marker genes56. In newborn and one-month old mice, we found that collagen α1(XI) was abundant in IVD and at the chondro-osseous junction of IVD and vertebrae, particularly concentrated in pre hypertrophic/hypertrophic chondrocytic cells. In long bone growth plates, Long et al.57 recently identified eight distinct cell clusters after unsupervised analysis of single cell (scRNAseq) of flow sorted hypertrophic chondrocytes from Col10a1Cre;Rosa26fs-tdTomato mice. At embryonic stage 16.5 (E16.5), Col11a1 expression was highest in cells with signatures of pre-hypertrophic to hypertrophic transition, and lowest in cells with osteogenic signatures (M. Hilton, personal communication)57. Taken together, these results suggest that collagen α1(XI) normally participates in maintaining growth plate cells in a hypertrophic, pre-osteogenic state, although little is known about its precise molecular function in that compartment, or in the IVD, during spinal development. Spines of Col11a1-deficient mice (cho/cho) show incompletely formed vertebral bodies, spinal curvatures, and decreased separation between vertebrae, which are themselves less mineralized than in wildtype mice36. Notably, common COL11A1 variants also have been associated with adult lumbar disc herniation (LDH) and degeneration (LDD), as well as DXA measured bone size, spinal stenosis and spondylolisthesis5860.

Pax1 is a well-described marker of early spine development, where it activates a gene expression cascade starting at E12.5-13.5 in mouse development44,61,62. Our data provided evidence that Pax1 regulates expression of Col11a1, as well as the AIS-associated genes Sox6 and Gpr126, in E12.5 tails of both male and female mice. The reduction of Col11a1 is consistent with a prior study of gene expression in flow-sorted GFP-labeled Pax1−/− embryonic IVD cells62. It is likely that Col11a1 expression in developing tail is directly activated by binding Sox transcription factors, as a prior genomic study using chromatin immunoprecipitation and sequencing (ChIP-seq) in rat chondrosarcoma cells identified super enhancers near the Col11a1 gene that were bound multiple times by the Sox9 and Sox663. The Sox5/6/9 trio is known to regulate many of the same genes as Pax162, but whether it mediates the Pax1 regulation of Col11a1 observed in our study is unknown.

In mouse postnatal spines, we observed co-localization of collagen α1(XI) and Pax1 protein specifically within the cartilaginous endplate-vertebral junction region that includes the vertebral growth plate. The endplate, which is important as the site allowing diffusion of nutrients from the circulation into the avascular inner IVD, harbors subpopulations of cells expressing type II collagen presumably organized by collagen type XI41,43. While the endplate is continuous with the vertebral growth plate in mice, it is important to note that in humans the endplate and vertebrae become distinctly separate structures with closure of the growth plates at puberty41. This is also the site of the ring apophyses that form the insertion of the IVD into vertebrae64. Lagging maturity of the ring apophysis, combined with mechanical forces across the IVD in the horizontal plane, has been proposed as an initiating factor leading to rotatory decompensation in the adolescent spine in AIS64,65. Recently, Sun et al. reported the discovery of a vertebral skeletal stem cell (vSSC) residing in the endplate and marked by expression of the genes Zic1 and Pax1, along with other cell surface markers66. These vSSCs also express high levels of Col11a1 (M. Greenblatt, personal communication). It is interesting to consider that AIS-associated genetic variation, perhaps together with mechanical forces, could alter the differentiation trajectory of this cell niche. Altogether, extant data and our results strongly suggest that cell populations at the IVD-vertebral junction region are relevant in AIS pathogenesis. Further investigation is warranted to understand the developmental programmes of cells in this region of the spine.

Matrix metalloproteinase 3, also known as stromolysin, is a secreted enzyme expressed in connective tissues and in regions of endochondral ossification67. MMP3 has degradative activity toward a variety of ECM components, including proteoglycans, fibronectin, laminin, but notably not type I collagen68. Additionally, in chondrocytes MMP3 also has been shown to translocate to the nucleus, where it activates transcription of connective tissue growth factor (CTGF/CCN2) by binding to an element known as transcription enhancer dominant in chondrocytes (TRENDIC)69,70. Our observations of a Col11a1-Mmp3 signaling axis in GPCs and CEP cells raise the possibility that Col11a1 variation may have consequences for both MMP3 enzymatic activity and for MMP3 mediated transcriptional programming in these cells. COL11A1 missense mutations, usually altering glycine or proline in Gly-X-Y repeats in the collagen α1(XI) helical domain as with COL11A1P1335L, are reported to be frequent in cutaneous squamous cell carcinomas and have been linked to transcriptional changes and tumor invasivenes50. The mechanisms by which GPCs or other cells sense such single amino acid changes in COL11A1 and induce subsequent transcriptional responses are unknown but may involve direct interaction with integrins in the pericellular space50.

We found that Col11a1 expression is sensitive to estrogen receptor status in GPCs, suggesting that in these cells it is responsive to estrogen signaling. Collagen α1(XI) is also a key player in organizing the pericellular space, which is critical for transmitting mechanical forces from the ECM to the cell71. Our findings support a new model wherein environmental cues, i.e. mechanical forces and estrogen signaling, intersect with genetic variation in a Pax1-Col11a1-Mmp3 signaling axis in the IVD-vertebral junction and thereby increase AIS susceptibility.

Methods

Discovery study

The cases in the discovery stage (USA TX: n=1,358) were recruited at Scottish Rite for Children as approved by the Institutional Review Board of the University Texas Southwestern Medical Center as previously described72. Subjects were genotyped on the Illumina HumanCoreExome BeadChip (Illumina, San Diego, CA, USA). For controls, we utilized 12,507 non-AMD GRU (non-age related macular degeneration general research use) subjects of the European ancestry downloaded from dbGaP web site (https://www.ncbi.nlm.nih.gov/gap/) from the International Age-Related Macular Degeneration Genomics Consortium study (IAMDGC: phs001039.v1.p1.c1; https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs001039.v1.p1). The subjects from the IAMDGC study were also genotyped on the Illumina HumanCoreExome Beadchip-24v1.0 platform73. We merged cases and controls and applied quality controls to the genotypes for 468,801 overlapping SNPs using PLINK.1.9. 74 as described in 18. In summary, samples with sex inconsistencies or from duplicated or related individuals or ancestral outliers as identified by principal component analysis (PCA) were removed, leaving 13,865 samples in the analysis. Genotypes were corrected for strand direction, and SNPs with call-rate per marker <95%, deviating from Hardy-Weinberg equilibrium (HWE) (cutoff P-value = 10−4), or with significant missingness rate between cases and controls (cutoff P-value = 10−4) were removed, leaving 341,759 SNPs in the analysis. Genotypes for SNPs across autosomal chromosomes were imputed using Minimac3 with the 1000G-Phase3.V.5 reference panel as described in the instructions available from the software website75. Protein coding changes were annotated with ANNOVAR using RefSeq-based transcripts76. External databases included allele frequencies from gnomAD77; variant pathogenicity in Clinvar78; Combined Annotation Dependent Depletion (CADD) scores79; Genomic Evolutionary Rate Profiling (GERP) scores 80 and protein domains in IntroPro 81. Only bi-allelic common (MAF>0.01) protein-altering SNPs with imputation quality Rsq >=0.3 within matrisome genes34 were included for further analysis. Matrisome genes used can be found in the Molecular Signature Database (MsigDB)82,83 (https://www.gsea-msigdb.org/gsea/msigdb/cards/NABA_MATRISOME). Genetic association for the imputed allele dosages in the discovery cohort (USA TX) was performed in Mach2dat 84 using logistic regression with gender and ten principal components as covariates. The genomic regions of the associated loci were visualized with LocusZoom software85 utilizing linkage disequilibrium (LD) information from 1000 Genomes EUR populations.

Meta-analysis study

For the meta-analysis stage we utilized four cohorts [USA MO: n=2,951(1,213 cases and 1,738 controls], Swedish-Danish populations [SW-D: n= 4,627 (1,631 cases and 2,996 controls), as described in 18,86], Japan [JP: n=79,211(5,327 cases and 73,884 controls)] and Hong Kong [HK: n=3,103 (990 cases and 2,113 controls)] to check significant candidates from the discovery study. Summary statistics across the discovery study and the 4 replication cohorts [total N=103,757 (10,519 cases and 93,238 controls)], were combined as previously described18 using METAL87. USA MO cohort: Whole exome sequencing (WES) data from 1,213 unrelated idiopathic scoliosis cases of European ancestry with spinal curvature greater than 10° Cobb angle were derived from the Adolescent Idiopathic Scoliosis 1000 Exomes Study (dbGAP accession number: phs001677), and recruited from St. Louis Children’s Hospital, and St. Louis Shriners Hospital for Children. Patients and/or parents provided study consent and IRB approval was obtained from each contributing institution. For controls, exome data from 1,738 unrelated samples of European ancestry were provided by investigators at Washington University School of Medicine in St. Louis, MO (dbGAP accession numbers: phs000572.v8.p4 and phs000101.v5.p1), and Oregon Health & Science University in Portland, OR (https://gemini.conradlab.org/). Exome data were aligned to the human genome reference (GRCh37) using BWA-MEM (v0.7.15). Variant calling of single nucleotide variants (SNVs) and insertion and deletion variants (INDELs) were generated first for each single sample in cases and controls and then combining all samples with joint genotyping method, described in GATK Best-Practices (Genome Analysis Toolkit (GATK v3.5) https://gatk.broadinstitute.org/hc/en-us/sections/360007226651-Best-Practices-Workflows). All cases and controls samples were classified as unrelated and of European ancestry using relationship inference 88 and principal component analysis 74. Association analysis of variants rs3753841 and rs1042704 were performed using logistic regression adjusting for sex and PCs in PLINK74.

JP cohort

Informed consents were obtained from all the subjects or their parents and the ethics committees of RIKEN and participating institutions approved this study. 5,327 case subjects were recruited from collaborating hospitals (Japanese Scoliosis Clinical Research Group) as previously described19. For controls, 73,884 subjects were randomly selected from the BioBank Japan Project, and subjects were genotyped on Illumina Human BeadChips as previously described 15,19. Imputation and association analyses in JP were performed as previously described89.

HK cohort

3,103 subjects were recruited at The Duchess of Kent Children’s Hospital as approved by the Institutional Review Board of the University of Hong Kong/Hospital Authority Hong Kong West Cluster (IRB approval number: UW 08-158). All 990 cases were characterized by Cobb angles greater than 40 degrees with onset age between 10 and 18 years old. Congenital, neuromuscular and syndromic scoliosis samples were excluded. We used 2,113 controls from the Chinese population with no spinal deformities on MRI scans 90. Cases and controls were genotyped using the Illumina Infinium OmniZhongHua-8 BeadChip and analyzed with GenomeStudio 2.0 software. The quality control approach adopted the GWA tutorial developed by Andries et al 91. The filtered genotyping data of cases and controls was phased and imputed using SHAPEIT 92 and IMPUTE2 93, respectively. Logistic model association analysis was performed using PLINK 1.974.

Stratification-by-sex test

To investigate sex-specificity in the COL11A1 and MMP14 loci, we performed stratification-by-sex analysis in the discovery study (USA_TX). Association for the imputed allele dosages in rs3753841 and rs37538 were computed separately for females (1,157 cases and 7138 controls) using logistic regression with ten principal components as covariates in Mach2dat84.

RNAseq of human tissues

RNAseq was performed as previously described94. Read counting and transcript quantification were performed using HTSeq95. Finally, reads were normalized using DESeq2 tools96 and TPM values were generated using the Kalisto pipeline97.

Animal studies

All mouse and rat work was conducted per IACUC approved protocols and in accordance with AALAC and NIH guidelines.

Generation of Pax1 knockout mice

Mouse work was approved by the UCSF IACUC, protocol number AN181381. Two gRNAs were designed to target the 5′ and 3′ ends of Pax1 gene (gRNA sequence shown in Figure 5-figure supplement 5a) using the gRNA design tool on the Integrated DNA Technologies (IDT, Newark, NJ, USA) website and selected based on low off-target and high on-target scores. The knockout allele was generated using i-GONAD98 as previously described 99.

To validate proper generation of knockout, mice were analyzed by genotyping with primers shown in Appendix A Key Resources Table, Sanger sequencing of PCR-amplified DNA, and Southern blot (Figure 5-figure supplement 5a). For Southern blot analyses, genomic DNA were treated with NcoI (Cat #R0193, New England Biolabs, MA, USA) and fractionated by agarose gel electrophoreses. Following capillary transfer onto nylon membranes, blots were hybridized with Digoxigenin (DIG)-labeled DNA probes (corresponding to chr2:147,202,083-147,202,444; mm9) amplified by the PCR DIG Probe Synthesis Kit (Cat #11636090910, Sigma-Aldrich, MO, USA). The hybridized probe was immunodetected with antidigoxigenin Fab fragments conjugated to alkaline phosphatase (Cat # 11093274910, Sigma-Aldrich, MO, USA) and visualized with a CDP star (Cat #11685627001, Sigma-Aldrich, MO, USA) according to the manufacturer’s protocol. Chemiluminescence was detected using the FluorChem E (Cat #92-14860-00, ProteinSimple, CA, USA).

Col11a1fl/fl mice

The Col11a1fl/fl mouse line was kindly provided by Dr. Lou Soslowsky with permission from Dr. David Birk. All Col11a1fl/fl and wildtype mice were handled per UTSW IACUC, protocol number 2016-101455.

Other mice

Cartilage was harvested from C57B46 wild type mice for siRNA-mediated knockdowns experiments.

Histologic methods

For thin cryostat sections, P0.5 mouse whole body was fixed in 4% paraformaldehyde for 6 hours followed by 10% sucrose for 12 hours, then transferred to 18% sucrose for 24 hours. Tissues were then embedded in optimal cutting temperature compound (OCT) and sectioned using low-profile blades on a Thermo Shandon Cryostar NX70 cryostat and all sections were lifted on APES clean microscope slides. For whole mount images, samples were treated similarly with the addition of 2% polyvinylpyrrolidone (PVP) during the cryoprotection step and frozen in 8% gelatin (porcine) in the presence of 20% sucrose and 2% PVP. Samples were sectioned at a thickness of 10 μm. Slides were stored at −80° C until ready for use. For P28 and older mice, spines were removed then fixed, decalcified, and embedded in OCT. Spines were processed by making 7 um thick lateral cuts the length of the spine.

Collagen α1(XI) was detected by immunohistochemistry staining using affinity-purified antisera against peptide (C) YGTMEPYQTETPRR-amide (Genescript, NJ, USA) as described40 and secondary horseradish peroxidase (HRP) conjugated affinity purified secondary antibody (Cat # AP187P, Millipore-Sigma Aldrich, MO, USA). Briefly, frozen sections were equilibrated to room temperature for 1 hour, then fixed with 4% paraformaldehyde in PBS at 4 degrees Celsius for 20 minutes. Slides were washed, treated with 3% H2O2 in methanol for 10 minutes to block endogenous peroxidase, washed, and transferred to PBS with 0.05% TWEEN 20 (Cat #P3563-10PAK, Sigma-Aldrich, MO, USA) pH 7.4. Slides were blocked with 0.5 % goat serum in PBS mix with 0.2% Triton 100 (Cat #T8787, Sigma-Aldrich, MO, USA) at room temperature for 1.5 hours. The primary collagen α1(XI) affinity purified antibody was applied at 0.40 mg/ml and slides were incubated overnight at 4°C. Afterward slides were washed in PBS Tween 20 for three times and treated with goat anti rabbit-HRP for 1.5 hours, then washed three times in PBS Tween 20. After applying 3,3’-diaminobenzidine (DAB) solution, slides were washed and counterstained with Mayer’s hematoxylin (Cat #MHS80, Sigma-Aldrich, MO, USA), washed, dehydrated, and mounted.

For collagen α1(XI) and Pax1 IF studies, P0.5 mice, P28 spine and ribs sections were fixed in 4% paraformaldehyde (PFA) for 20 minutes then washed with PBS + 0.1% Triton three times, before incubation with 10% normal goat serum in PBS + 0.1% Triton for 30 min to block the background. Slides were incubated with goat anti-mouse collagen α1(XI) antibody at 1:500 dilution and mouse anti-rat Pax1(Cat #MABE1115M, Sigma-Aldrich, MO, USA), in PBS + 0.1% Triton + 1% normal goat serum at 4 °C overnight. Secondary antibodies used were 1:5000 anti-rat Alexa488 and anti mouse Alexa594 conjugated antibodies (Cat #A32740 Invitrogen, CA, USA). The sections were mounted using ProLong Gold with DAPI (Cat #S36964 Invitrogen, CA, USA) for imaging as described100.

#S36964 Invitrogen, CA, USA)). All images were taken with Carl Zeiss Axio Imager.M2 fluorescence microscope (Zeiss,Oberkochen, DE).

Rib cartilage and IVD cell culture

Mouse rib cartilage cells were isolated from the rib cage and sternum of P0.5 mice. Rat IVD was removed intact from 1-month female rats and immediately separated into NP, AF, and CEP isolates. Subsequently, tissues were incubated and shaken with 2 mg/ml Pronase solution (Cat #10165921001 Sigma-Aldrich, Inc., St. Louis, MO, USA) for 1 hour, followed by 1.5 hours digestion with 3 mg/ml Collagenase D solution. (Cat #11088882001 Sigma-Aldrich, Inc., St. Louis, MO, USA), then 5 hours digestion with 0.5 mg/ml Collagenase D solution before 3 times PBS wash. Filtered, dissociated cells were seeded in Dulbecco’s modified Eagle’s medium (DMEM; Cat #MT15017CV Thermo Fisher Scientific, MA, USA) containing 10% fetal bovine serum (FBS), 100 μg/ml streptomycin and 100 IU/ml penicillin. Remaining cartilage tissues underwent further digestion in 86 U/ml type 1 collagenase (Cat # SCR103 Sigma-Aldrich, Inc., St. Louis, MO, USA) overnight and cells collected.

SV40 immortalization of primary chondrocytes

Col11a1fl/fl mouse rib cartilage cells were isolated from the rib cage and sternum of P0.5 mice. The cells were transduced with pRRLsin-sv40 T antigen-IRES-mCherry lenti-virus101 for 48 hours, then sorted for mCherry-positive cells by flow cytometry. mCherry-positive cells were then infected with plv-eGFP, plv-eGFP-COL11A1-HA, plveGFP-COL11A1P1335L-HA constructs. After expansion, GFP-positive cells were sorted by flow cytometry.

Adenovirus treatment

SV40 induced Col11a1fl/fl mouse rib cartilage cells seeded on 6 wells plate and grown to 50% confluency. Afterward, cells were treated with 2ul Ad5-CMV-cre adenovirus (titer 1.8×1011pfu/ml) and Ad5-CMV-eGFP (titer 1.65×1010pfu/ml) as control. Both virus strains were from the Gene Vector Core facility, Baylor College of Medicine. After 48-hours the cells were harvested for mRNA and protein lysate.

RNA-seq and qRT-PCR

For Pax1 knockout studies, total RNA was collected from E12.5 tails using TRIzol (Cat #15596026, Thermo Fisher Scientific, MA, USA) and converted to cDNA using ReverTra Ace qPCR-RT master mix with genomic DNA (gDNA) remover (Cat #FSQ-301, Toyobo, Osaka, Japan). Sequencing was done using an Illumina Novaseq platform and the data were analyzed using Partek Flow (Version 10.0) and Gene ontology102. qPCR was performed using SsoFast EvaGreen supermix (Cat #1725205, Bio-Rad, CA, USA). Primer sequences used for qPCR are shown in Appendix A, Key Resources table.

To quantify the expression level of Col11a1, Mmp3, and marker genes in IVD compartments and rib cartilage, cultured cells were collected in RNeasy (Qiagen, Inc.) for RNA purification. Taqman Reverse Transcription Kit (Cat #4387406 Thermo Fisher Scientific, MA, USA) was used to reverse transcribe mRNA into cDNA. Following this, RT-qPCR was performed using a Power SYBR Green PCR Master Mix Kit (Cat #1725271, Bio-Rad, CA, USA). The primer sequences for the genes used in this study are listed in the Appendix A Key Resources Table. Gene expression was calculated using the ΔΔCT method after normalizing to GAPDH.

siRNA knockdown

Mouse rib cartilage cells seeded in 6 well-plates were 60-80% confluent at transfection. Lipofectamine RNAiMAX reagent (Cat #13778030 Thermo Fisher, Waltham, MA, USA) was diluted (9ul in 500ul) in Opti-MEM Medium (Cat # 31985070 Thermo Fisher MA, USA). 30 pmol siRNA was diluted in 500ul Opti-MEM medium, then added to diluted Lipofectamine RNAiMAX Reagent. siRNA-lipid complex was added to cells after 5-minute incubation at room temperature. Cells were harvested after 72 hours.

Western blotting

For Mmp3 western blotting, a total of 30 µg protein mixed with SDS-PAGE buffer was loaded on 12% SDS-polyacrylamide gel for electrophoresis. For collagen α1(XI) western blotting, 50 µg protein mixed with SDS-PAGE buffer was loaded on 4-20% SDS-polyacrylamide gel. The separated proteins were then transferred to nitrocellulose membranes (Cat #77010 Thermo Fisher Waltham, MA, USA) at 100 volts for 2-3 hours. The membrane was first incubated with blocking buffer containing 5% defatted milk powder, and then exposed to 0.1 mg/mL mouse anti-rabbit Mmp3 (Cat #ab214794 Abcam, Cambridge, MA, USA) or anti-rabbit Col11a1 (Cat #PA5-38888 Thermo Fisher Waltham, MA, USA) overnight. The samples were then washed thoroughly with TBS buffer, followed by incubation with HRP-labeled antirabbit IgG secondary antibodies 1:5000, (Cat #32460 Thermo Fisher Waltham, MA, USA) overnight. The membranes were then washed with TBS buffer. GAPDH was detected by a rabbit anti-mouse antibody (Cat #14C10 Cell Signaling, MA, USA) and used as the internal loading control.

Data availability

Summary data for GWAS3 is deposited in the NHGRI-EBI GWAS catalog (accession #GCST90179478). Otherwise, data generated or analyzed in this study are included in the manuscript.

Acknowledgements

We thank the patients and their families who participated in these studies. We are grateful to the clinical investigators who referred patients into the study from Japan Scoliosis Clinical Research Group (JSCRG), Scottish Rite for Children Clinical Group (SRCCG). We are grateful for the ScoliGeneS study group for the help with patient recruitment and analyses of the Swedish-Danish cohort. The names and affiliations for JSCRG, TSRHCCG, ScoliGeneS are included Appendix B. We also thank Drs. Carlos Cruchaga, Sanjay Jain, Matthew Harms, and Don Conrad for allowing us to use exome data from their cohort studies. The University of Pennsylvania is acknowledged for providing the Col11a1fl/fl mouse line. We thank Drs. D. Burns and B. Evert for their help interpreting histological studies. We thank Dr. M. Hilton for sharing scRNAseq data as described in reference (51). This study was funded by JSPS KAKENHI Grants (22H03207 to SI), the Swedish Research Council (number K-2013-52X-22198-01-3 and 2017-01639 to PG and EE), the regional agreement on medical training and clinical research (ALF) between Stockholm County Council and Karolinska Institutet (to PG), Center for Innovative medicine (CIMED), Karolinska Institutet (to PG), the Department of Research and Development of Vasternorrland County Council (to PG), and Karolinska Institutet research funds (to PG) and Research Grants Council of Hong Kong (No:17114519 to YQS), The National Institutes of Health GM127390 (to NG) and the Eunice Kennedy Shriver National Institute of Child Health & Development of the National Institutes of Health (P01HD084387 to CAW and NA). Content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH. The data used for the analyses described in the USA MO cohort were obtained from the database of Genotypes and Phenotypes (dbGaP) Adolescent Idiopathic Scoliosis 1000 Exomes Study (phs001677.v1.p1) which included investigators at Washington University in St Louis (M. Dobbs), Shriners Hospital for Children, St Louis (M. Dobbs), University of Colorado (N. Miller), University of Iowa (S. Weinstein and J. Morcuende), University of Wisconsin (P. Giampietro), and Hospital for Special Surgery (C. Raggio). Sequencing of these samples was funded by NIH NIAMS 1R01AR067715 (M. Dobbs and CG). Use of the dbGAP dataset phs001039.v1.p1 is gratefully acknowledged. All contributing sites and additional funding information for dbGAP dataset phs001039.v1.p1 are acknowledged in this publication: Fritsche et al. (2016) Nature Genetics 48 134–143, (doi:10.1038/ng.3448); The International AMD Genomics consortium’s web page is: http://eaglep.case.edu/iamdgc_web/, and additional information is available on: http://csg.sph.umich.edu/abecasis/public/amd015/. The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing computing resources that have contributed to the research results reported within this paper. URL: http://www.tacc.utexas.edu

Author contributions

H.Y., A.M.K., A.U., X-Y. L, N.O, Y.K., E.E., Y.F., L.A., Y.H.K., R.C., R.S, Y.Z., J.P. performed genetic, bioinformatic, or proteomic analyses or molecular biology experiments. N.V.G., B.M.E., J.P.Y.C., and J.A.H. provided technical or clinical advice. C.T., Y.-Q.S., C.A.G., P.G., and S.I. supervised genetic experiments and provided summary data. J.J.R., N.A., and C.A.W. supervised the work. C.A.W conceived the study and wrote the manuscript. All co-authors reviewed and approved the final manuscript.

Competing interests

The authors declare no competing interests.

Appendix A. Key Resources Table

siRNA and Primer Sequences

Appendix B. Clinical groups

Texas Scottish Rite Hospital for Children Clinical Group (TSRHCCG)

Lori A. Karol1, Karl E. Rathjen1, Daniel J. Sucato1, John G. Birch1, Charles E. Johnston III1, Benjamin S. Richards1, Brandon Ramo1, Amy L. McIntosh1, John A. Herring1, Todd A. Milbrandt2, Vishwas R. Talwakar3, Henry J. Iwinski42, Ryan D. Muchow3, J. Channing Tassone4, X.-C. Liu4, Richard Shindell5, William Schrader6, Craig Eberson7, Anthony Lapinsky8, Randall Loder9 and Joseph Davey10

  1. Department of Orthopaedic Surgery, Texas Scottish Rite Hospital for Children, Dallas, Texas, USA.

  2. Department of Orthopaedic Surgery, Mayo Clinic, Rochester, Minnesota, USA

  3. Department of Orthopaedic Surgery, Shriners Hospitals for Children, Lexington, Kentucky, USA.

  4. Department of Orthopaedic Surgery, Children’s Hospital of Wisconsin, Milwaukee, Wisconsin, USA.

  5. OrthoArizona, Phoenix, Arizona, USA.

  6. Departments of Orthopedics, Sports Medicine, and Surgical Services, Akron Children’s Hospital, Akron, Ohio, USA.

  7. Pediatric Orthopaedics and Scoliosis, Hasbro Children’s Hospital, Providence, Rhode Island, USA.

  8. University of Massachusetts Memorial Medical Center, Worcester, Massachusetts, USA.

  9. Indiana University-Purdue University Indianapolis, Indianapolis, Indiana, USA.

  10. University of Oklahoma Health Sciences Center, Oklahoma City, Oklahoma, USA.

Japan Scoliosis Clinical Research Group (JSCRG)

Kota Watanabe1, Nao Otomo1, 2, Kazuki Takeda1, 2, Yoshiro Yonezawa1, 2, Yoji Ogura1, 2, Yohei Takahashi1, 2, Noriaki Kawakami3, Taichi Tsuji4, Koki Uno5, Teppei Suzuki5, Manabu Ito6, Shohei Minami7, Toshiaki Kotani7, Tsuyoshi Sakuma7, Haruhisa Yanagida8, Hiroshi Taneichi9, Satoshi Inami9, Ikuho Yonezawa10, Hideki Sudo11, Kazuhiro Chiba12, Naobumi Hosogane12, Kotaro Nishida13, Kenichiro Kakutani13, Tsutomu Akazawa14, Takashi Kaito15, Kei Watanabe16, Katsumi Harimaya17, Yuki Taniguchi18, Hideki Shigemats19, Satoru Demura20, Takahiro Iida21, Ryo Sugawara22, Katsuki Kono23, Masahiko Takahata24, Norimasa Iwasaki24, Eijiro Okada1, Nobuyuki Fujita1, Mitsuru Yagi1, Masaya Nakamura1, Morio Matsumoto1

  1. Department of Orthopaedic Surgery, Keio University School of Medicine, Tokyo, Japan

  2. Laboratory of Bone and Joint Diseases, Center for Integrative Medical Sciences, RIKEN, Tokyo, Japan

  3. Department of Orthopaedic Surgery, Meijo Hospital, Nagoya, Japan

  4. Department of Orthopaedic Surgery, Toyota Kosei Hospital, Nagoya, Japan

  5. Department of Orthopaedic Surgery, National Hospital Organization, Kobe Medical Center, Kobe, Japan

  6. Department of Orthopaedic Surgery, National Hospital Organization, Hokkaido Medical Center, Sapporo, Japan

  7. Department of Orthopaedic Surgery, Seirei Sakura Citizen Hospital, Sakura, Japan

  8. Department of Orthopaedic Surgery, Fukuoka Children’s Hospital, Fukuoka, Japan

  9. Department of Orthopaedic Surgery, Dokkyo Medical University School of Medicine, Mibu, Japan

  10. Department of Orthopaedic Surgery, Juntendo University School of Medicine, Tokyo, Japan

  11. Department of Advanced Medicine for Spine and Spinal Cord Disorders, Hokkaido University Graduate School of Medicine, Sapporo, Japan

  12. Department of Orthopaedic Surgery, National Defense Medical College, Tokorozawa, Japan

  13. Department of Orthopaedic Surgery, Kobe University Graduate School of Medicine, Kobe, Japan

  14. Department of Orthopaedic Surgery, St. Marianna University School of Medicine, Kawasaki, Japan

  15. Department of Orthopaedic Surgery, Osaka University Graduate School of Medicine, Suita, Japan

  16. Department of Orthopaedic Surgery, Niigata University Hospital, Niigata, Japan

  17. Department of Orthopaedic Surgery, Graduate School of Medical Sciences, Kyushu University Beppu Hospital, Fukuoka, Japan

  18. Department of Orthopaedic Surgery, Faculty of Medicine, The University of Tokyo, Tokyo, Japan

  19. Department of Orthopaedic Surgery, Nara Medical University, Nara, Japan

  20. Department of Orthopaedic Surgery, Kanazawa University School of Medicine, Kanazawa, Japan

  21. Department of Orthopaedic Surgery, Dokkyo Medical University Koshigaya Hospital, Koshigaya, Japan

  22. Department of Orthopaedic Surgery, Jichi Medical University, Simotsuke, Japan

  23. Department of Orthopaedic Surgery, Kono Othopaedic Clinic, Tokyo, Japan

  24. Department of Orthopaedic Surgery, Hokkaido University, Sapporo, Japan

Scoliosis and Genetics in Scandinavia (ScoliGeneS) study group

Tian Cheng1, Juha Kere2, Aina Danielsson3, Kristina Åkesson4, Ane Simony5, Mikkel Andersen5, Steen Bach Christensen5, Maria Wikzén6, Luigi Belcastro6

  1. Department of Clinical Sciences and Technology, Karolinska Institutet, Huddinge, Sweden

  2. Department of Biosciences and Nutrition, Karolinska Institutet, Huddinge, Sweden

  3. Department of Orthopedics, Sahlgrenska University Hospital, Gothenburg, Sweden

  4. Department of Orthopedics and Clinical Sciences, Lund University, Skane University Hospital, Malmö, Sweden

  5. Sector for Spine Surgery and Research, Middelfart Hospital, Middelfart, Denmark

  6. Department of Reconstructive Orthopaedics, Karolinska University Hospital, Stockholm, Sweden.