Exome sequencing analysis identifies frequent oligogenic involvement and FLNB variants in adolescent idiopathic scoliosis.

BACKGROUND: Adolescent idiopathic scoliosis (AIS) is a genetically heterogeneous disease characterised by three-dimensional deformity of the spine in the absence of a congenital spinal anomaly or neurological musculoskeletal disorder. The clinical variability and incomplete penetrance of some genes linked with AIS indicate that this disease constitutes an oligogenic trait.
OBJECTIVE: We aimed to explore the oligogenic nature of this disease and identify novel AIS genes.
METHODS: We analysed rare damaging variants within AIS-associated genes by using exome sequencing in 40 AIS trios and 183 sporadic patients.
RESULTS: Multiple variants within AIS-associated genes were identified in eight AIS trios, and five individuals harboured rare damaging variants in the FLNB gene. The patients showed more frequent oligogenicity than the controls. In the gene-based burden test, the top signal resided in FLNB. In functional studies, we found that the AIS-associated FLNB variants altered the protein's conformation and subcellular localisation and its interaction with other proteins (TTC26 and OFD1) involved in AIS. The most compelling evidence of an oligogenic basis was that the number of rare damaging variants was recognised as an independent prognostic factor for curve progression in Cox regression analysis.
CONCLUSION: Our data indicate that AIS is an oligogenic disease and identify FLNB as a susceptibility gene for AIS. © Author(s) (or their employer(s)) 2020. Re-use permitted under CC BY-NC. No commercial re-use. See rights and permissions. Published by BMJ.

Introduction

Adolescent idiopathic scoliosis (AIS) is characterised by three-dimensional deformities of the spine. It affects 1%–3% of all children and is the most common type of spine deformity. Although the aetiology and pathogenesis of AIS remain unclear, twin and family studies have indicated the important role of genetic factors in spinal curve formation and progression. Both genetic and environmental factors are known to contribute to sporadic AIS. Recent genome-wide association studies have identified several common variants, including susceptibility loci near LBX1,1 GPR126 2 and PAX1,3 which are associated with AIS. Exome sequencing has also demonstrated the role of rare variants in genes, including FBN1,4 POC5 5 and AKAP2 6 in AIS onset and progression. However, the loci identified thus far only account for a small amount of AIS heritability, and the genetic background of the majority of AIS cases remains poorly explained.

The loci mentioned previously can be inherited in an autosomal-dominant, autosomal-recessive or X-linked manner in familial AIS. The incomplete penetrance of some genes (FBN1 4 and HSPG2 5) with autosomal-dominant inheritance patterns has also been described. Notably, a heterozygous PTK7 missense variant (which provides a faithful developmental zebrafish model of idiopathic scoliosis (IS))7 was identified in one patient with IS and his father without spinal deformity.8 In another study, pathway burden analysis of exome sequence data indicated that patients with AIS harboured multiple rare variants within extracellular matrix genes and that the burden of variants influenced the clinical features of the patients,9 which supports an oligogenic or polygenic inheritance model of AIS.10 11 In addition, in a study of extended families in Utah,10 variable recurrence risk and scoliosis phenotypes (curve severity and type) were observed within families in which multiple individuals were affected, indicating polygenic inheritance of AIS. Another study of multiplex families demonstrated that AIS is a disease with multigenic, multifactorial inheritance in which a greater genetic load is required for men to be affected.11

To identify novel AIS genes and explore the oligogenic nature of the disease, we carried out an exome sequencing study of both AIS trios and patients with sporadic AIS and tested for an association of rare damaging variants with AIS. By studying the gene–gene interactions identified in some AIS trios, we further extended the genetic architecture of AIS.

Methods

Cohort description

Patient with AIS were consecutively recruited from Shanghai Changzheng Hospital. A total of 40 AIS trios (two parents without AIS and one child) and 183 patients with sporadic AIS were recruited (the demographics and clinical characteristics of the patients are summarised in online supplementary table 1). The diagnostic criteria for AIS were as follows: (1) spinal curve of >10° at first presentation and (2) no congenital spinal anomalies or scoliosis secondary to other disorders, including Marfan’s syndrome or neurological disorders. All the patients diagnosed with AIS underwent a physical examination of the spine, including a bending test with a scoliometer. Neurological examination (abdominal reflex test and MRI) was only performed on those patients who were suspected of having an underlying disorder (signs of pyramidal irritation and signs of cerebellar disorder). For the AIS trios, the additional inclusion criteria included the following: (1) the patient and both parents were living, and their DNA was available; (2) both parents showed a normal spine on X-ray examination (mean curve of 4.4°, range 0–8.3°); and (3) no other hereditary diseases were identified. A total of 153 age-matched, sex-matched and ethnicity-matched control subjects were also included as in-house controls. All patients were followed up regularly every 3 months and underwent whole-spine standing anteroposterior and lateral X-ray examination until skeletal maturity (18 years old or Risser sign=5). X-ray examination was also performed for the controls to rule out scoliosis. Blood samples were collected from both patients and in-house controls. In addition, we used the 222 exome data of 222 Han Chinese individuals from the 1000 Genomes Project as controls, including data from the Han Chinese in South China, and Han Chinese in Beijing, China, groups.

Exome sequencing and variant annotation

Exome sequencing was performed at 100×coverage by commercial providers (iGene TechTM, China). Genomic DNA was isolated from peripheral blood samples using a QIAamp DNA Blood kit (Qiagen, Germany) according to the manufacturer’s protocol, and the exome was enriched with a TruSeq Exome Enrichment kit (Illumina, California, USA) following the manufacturer’s protocol. Then, an Illumina HiSeq 2500 sequencer was employed to sequence the human exome. Exome capture was performed using the TargetSeq Enrichment kit (iGeneTech, China). Exome data were analysed using the GenPipes DNAseq pipeline (https://bitbucket.org/mugqic/genpipes). Paired-end sequencing reads were trimmed using Trimmomatic to obtain a high-quality set of reads for sequence alignment file generation, and the trimmed reads were aligned to the reference genome (hg19) by Burrows-Wheeler Alignment-MEM algorithm. Genotype and annotation data from the 1000 Genomes Project were used to create files of the gene variants from each of the 222 Han Chinese individuals sequenced in the project. Exome samples from both case and control subjects were processed using a consistent alignment and variant calling pipeline. We used samtools bedcov V.1.3.1 to identify the high coverage regions with the SureSelect V.5 all exon plus untranslated regions (UTR) kit manifest Browser Extensible Data (BED) coordinates. In addition, we used a BED coordinate list that contained the high coverage, unique genomic regions of the exome kit to unify the exome region of the in-house and external data (the summary data and technical characters of exome sequencing data are summarised in online supplementary table 2).

After alignment, we achieved an average mapped read depth of 42×–146× (median 98.4×), and the breadth of coverage was highly comparable between the patients and the control subjects. The aligned reads were then processed according to the Genome Analysis Toolkit (GATK). Variant-based quality control was initially carried out using GATK variant quality score recalibration. We selected the truth sensitivity tranche of 99.5% for single-nucleotide variants. Additionally, we applied genotype-level quality control using KGGseq, in which set low-quality genotypes, with a genotype quality of <20 or coverage by fewer than eight reads (read depth of <8), were set as missing to avoid false-positive and false-negative calls.

We used uniform procedures for variant calling in case and control subjects to avoid introduction of technical bias in the study. The variant call set was then annotated according to the RefGene gene annotations, pathogenicity (MutationTaster and SIFT) and population frequencies (1000 Genomes Project phase 3, Genome Aggregation Database (GnomAD) and National Heart, Lung and Blood Institute (NHLBI) Exome Sequencing Project databases). We combined two in silico prediction tools, MutationTaster and SIFT,12 to evaluate the deleterious properties of missense mutations to achieve low false-negative predictions; the combination of these two prediction tools has previously been used to assess the potential impact of mutations on the function of the filamin B (FLNB) protein.13 We defined rare damaging variants according to the following criteria: (1) missense, nonsense, frame-shift or splice-site variants and (2) variants with a minor allele frequency of <1% in any of the following public databases: 1000 Genomes Project phase 3, GnomAD and NHLBI Exome Sequencing Project.

AIS-associated genes and their interactome

We searched PubMed using the search terms ‘(idiopathic scoliosis) and gene’. Linkage and association studies, genome-wide association studies and exome sequencing studies published in English whose full text was available were further analysed. Twenty-eight candidate genes were implicated in human IS were identified. Genes whose inactivation or mutation results in an IS phenotype in mice according to the Mouse Genomics Informatics databases (http://www.informatics.jax.org/) or Zebrafish, according to the Zebrafish Information Network (https://zfin.org/), were also considered AIS-associated genes (online supplementary table 3). We defined the genes encoding proteins showing protein–protein interactions with significant AIS genes in the InWeb databases (https://omictools.com/inweb-im-tool) as the interactome.

Gene-based burden test for rare variant association

We assessed the enrichment of rare damaging variants of AIS-associated genes in patients with AIS, in-house controls and controls from the 1000 Genomes Project and performed gene-based association tests using multivariate and collapsing tests.

In silico protein structure analysis

To analyse the possible conformational changes in the molecular structures of human Filamin-B, mutation simulations were performed based on the wild-type structures retrieved from the Protein Data Bank. The residues were mutated using the mutate module in the Protein Builder application of MOE 2018.01 software. The Amber 10: EHT force field was used to minimise the wild-type and mutant protein structures. Intramolecular interactions and conformational changes between the wild-type and mutant FLNB protein structures were analysed.

Plasmid constructs

Full-length FLNB with a Flag tag and OFD1 and TTC26 with a myc tag were inserted into the pCDH vector. The clones of FLNB, OFD1 and TTC26 were concordant with the reference sequences NM_001457.3, NM_003611.3 and NM_001144920.2, respectively. PCR mutagenesis was performed using forward and reverse oligonucleotide primers containing single-nucleotide alterations, paired with 5′ and 3′ external primers. The products were subcloned into the pcDNA3.1(−) or Flag (or myc)-tagged PCDH vector. All constructs were verified by DNA sequencing.

Cell culture and transfection

HEK293T cells were cultured in Dulbecco's Modified Eagle Medium (DMEM) containing 10% foetal bovine serum supplemented with penicillin and streptomycin. For plasmid transfection, TransIT-X2 (Mirus) was used according to the manufacturer’s instructions.

Immunoprecipitation assay and western blotting

Cell pellets were lysed in ice-cold lysis buffer (25 mM Tris–HCl, pH 7.4, 150 mM NaCl, 0.4% digitonin, 1 mM EDTA and protease inhibitors). The supernatant was collected by centrifugation for 30 min at 12 000 g at 4°C and further precleared using protein-G Sepharose for 1 hour. After the removal of protein-G beads, the precleared supernatant was incubated with protein-G beads and 1 µg of the indicated primary antibodies (myc tag (2278 1:200 dilution for immunoprecipitation)) or the IgG control overnight at 4°C. After washing, the Sepharose beads were boiled in 1× sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE) loading buffer. Proteins were detected by western blotting.

For western blotting, protein samples were subjected to standard SDS-PAGE and transferred to polyvinylidene fluoride membranes. The membranes were blocked in 2.5% non-fat milk for 1 hour and incubated overnight at 4°C with primary antibodies (Flag tag (8146, dilution 1:50 for immunoprecipitation)). After washing with Tris-buffered saline, 0.05% TWEEN (TBST) three times for 10 min each, the membranes were incubated with secondary antibodies for 1 hour at room temperature. After washing with TBST three times for 10 min each, the membranes were developed with chemical luminescence reagents (Bio-Rad). Images were obtained using a ChemiDoc Touch Imaging System (Bio-Rad).

Statistical analyses

Statistical significance was determined with the two-sided unpaired Student’s t-test or two-way analysis of variance using the Statistical Package for Social Sciences V.19.0. Progression-free survival was calculated from the time of curve onset until the development of a curvature of at least 30°, skeletal maturity (18 years old or Risser sign=5) or the date of the last clinical visit. Kaplan-Meier analysis and Cox proportional hazards models were used to estimate the survival distribution.14 Univariate and multivariate Cox regression analyses were performed to identify the prognostic values of sex, age of curve onset, curve magnitude at presentation, Risser sign, body mass index and the number of rare damaging variants. P values of <0.05 were considered statistically significant. All experiments were replicated at least three times, and contiguous data are shown as means with SDs.

Result

A frequent oligogenic involvement in AIS trios

First, we annotated our exome data from 40 AIS trios for variants in AIS-associated genes. We identified 38 rare damaging variants in 17 genes in patients with AIS (online supplementary table 4). In particular, we identified variants in more than one AIS-associated gene in eight trios (figure 1). Subsequently, we performed a binomial test that demonstrated that the frequency of trios with oligogenic bases was higher than expected on the basis of chance (p=5.82×10−4). We also highlighted five individuals harbouring rare damaging variants in FLNB (Mendelian Inheritance in Man: 603381, figure 1).

Figure 1

Pedigrees with multiple rare damaging variants within AIS-associated genes. Filled symbols for men (squares) and women (circles) denote affected individuals, and empty symbols indicate unaffected individuals. Individuals with heterozygous variants are indicated with ‘m/−’, while ‘−/−’ indicates wild type. AIS, adolescent idiopathic scoliosis.

Variants in AIS-associated genes contribute to oligogenic AIS

To assess oligogenicity, we enrolled 183 patients with sporadic AIS and 153 in-house controls for exome sequencing and applied the data of the 222 Chinese participants from the 1000 Genomes Project for comparison. We identified variants of AIS-associated genes in 50.22% (112/223) of patients with AIS. Notably, 34.82% (39/112) of these patients carried more than one variant. Among the 39 patients with AIS harbouring multiple variants, only one carried biallelic variants in the AKAP2 gene (ie, 0.89% homozygosity or compound heterozygosity), and 38 carried at least two variant alleles of different genes (ie, 33.93% oligogenicity). Taken together, oligogenicity was more frequent in the patients with AIS than in the controls from the 1000 Genomes Project (37/186 vs 8/214; OR 5.321, 95% CI 2.417 to 11.714, p=6.00E-06) or the in-house controls (37/186 vs 3/150; OR 9.946, 95% CI 3.008 to 32.893, p=2.00E-06) (table 1 and online supplementary table 5).

Table 1

Number of alleles with mutations within AIS-associated genes in individuals with AIS and controls

	Patients with mutations >1 (n)	OR (95% CI)(AIS vs controls)	P value(AIS vs controls)
AIS (n=223)	37	1
1000 Genomes Project controls (n=222)	8	5.321 (2.417 to 11.714)	6.00E-06
In-house controls (n=153)	3	9.946 (3.008 to 32.893)	2.00E-06
Combined controls (n=375)	11	6.583 (3.282 to 13.201)	7.44E-09

AIS, adolescent idiopathic scoliosis.

A gene-based burden test identified FLNB as a susceptibility gene for AIS

Next, we focused on AIS-associated genes to identify the top AIS candidate genes. Notably, rare damaging variants in FLNB (p=6.30E-05), TTC26 (p=3.07E-04), PTK7 (p=5.37E-03), CNTNAP2 (p=7.12E-03), FBN1 (p=7.12E-03) and TTLL3 (p=2.33E-02) were over-represented in patients with AIS (table 2). We identified rare damaging FLNB variants in 25/223 (11.21%) of patients with AIS, 16/25 (64.00%) of whom exhibited an additional variant of an AIS-associated gene (eg, TTC26, PTK7, TTLL3 or OFD1; table 3).

Table 2

Rare damaging variants gene-based burden tests for AIS-related genes

Rank	Gene	Count of rare damaging variants		Frequency of rare damaging variants in AIS subjects (%)	Frequency of rare damaging variants in controls subjects (%)	OR (95% CI)(AIS vs controls)	P value(AIS vs controls)
		AIS(n=223)	Controls* (n=375)
1	FLNB	25	11	11.21	2.93	4.178 (2.014 to 8.670)	6.30E-05
2	TTC26	16	5	7.17	1.33	5.720 (2.066 to 15.838)	3.07E-04
43	PTK7	11	4	4.93	1.07	4.813 (1.514 to 15.302)	5.37E-03
4	CNTNAP2	8	2	3.59	0.53	6.940 (1.460 to 32.976)	7.12E-03
5	FBN1	8	2	3.59	0.53	5.720 (2.066 to 15.838)	7.12E-03
6	TTLL3	8	3	3.59	0.80	4.614 (1.211 to 17.577)	2.33E-02
7	FBN2	8	4	3.59	1.07	3.451 (1.027 to 11.596)	6.51E-02
8	IFT122	6	4	2.69	1.07	2.565 (0.716 to 9.189)	1.86E-01

P value is calculated using Fisher’s exact test.

*In-house controls (n=153) and 1000 Genomes Project controls (n=222) were combined in gene-based burden tests.

AIS, adolescent idiopathic scoliosis.

Table 3

Rare damaging FLNB variants detected in patients with AIS

Sample ID	Position	reID	Genotype	Nucleotidechange*	Protein change	Max freq†	Additional mutations of AIS-associated genes
98–1	Chr3:58 064 498		Heterozygotes	c.G596A	p.R199Q	N	ESR2, POC5
98–4	Chr3:58 087 993		Heterozygotes	c.G1409A	p.R470Q	2.39E-05	IFT27
98–7	Chr3:58 133 944		Heterozygotes	c.C5740T	p.R1914W	5.00E-04	FBN1, TTLL3
98–13	Chr3:58 134 076		Heterozygotes	c.C5872T	p.P1958S	N	PTK7
98–16	Chr3:58 139 231	rs199939739	Heterozygotes	c.C6497T	p.T2166M,	5.66E-05	CLUAP1
98–32	Chr3:58 148 980	rs186952950	Heterozygotes	c.G7121A	p.R2374H	5.30E-05	PTK7
98–73	Chr3:58 134 496	rs563096120	Heterozygotes	c.G6008A	p.R2003H	4.38E-05	OFD1
98–84	Chr3:58 112 379		Heterozygotes	c.C4112T	p.S1371L	4.37E-05
Trio-22	Chr3:58 090 893		Heterozygotes	c.G1697T	p.R566L	N	TTC26
Trio-25	Chr3:58 090 892	rs778577280	Heterozygotes	c.C1696T	p.R566W	N	PTK7
Trio-27	Chr3:58 141 758		Heterozygotes	c.G6844A	p.A2282T	N	TTC26
Trio-31	Chr3:58 090 941		Heterozygotes	c.T1745C	p.L582P	N	IFT122
Trio-32	Chr3:58 109 123	rs199959926	Heterozygotes	c.G3430C	p.E1144Q	7.21E-03	COL11A2, PTK7
24–3	Chr3:58 129 322	rs200677473	Heterozygotes	c.A5407T	p.M1803L	6.50E-03
24–17	Chr3:58 145 335	rs754457328	Heterozygotes	c.G6943A	p.A2315T	1.06E-04
24–21	Chr3:58 155 406	rs761994878	Heterozygotes	c.A7507G	p.S2503G	1.74E-04
59–1	Chr3:58 108 868		Heterozygotes	c.G3175A	p.A1059T	N	IFT88, CHL1, PTK7
59–2	Chr3:58 109 276	rs200993986	Heterozygotes	c.G3583A	p.V1195M	8.20E-03	TTLL3
59–10	Chr3:58 110 215		Heterozygotes	c.A3881G	p.Y1294C	N	TTC26
59–15	Chr3:58 121 758		Heterozygotes	c.C4724T	p.T1575M	1.00E-04
59–19	Chr3:58 121 848	rs201630300	Heterozygotes	c.G4814A	p.R1605H	4.00E-03
59–32	Chr3:58 129 322	rs200677473	Heterozygotes	c.A5407T	p.M1803L	6.50E-03
59–47	Chr3:58 133 944		Heterozygotes	c.C5740T	p.R1914W	5.00E-04
43–4	Chr3:58 129 322	rs200677473	Heterozygotes	c.A5407T	p.M1803L	6.50E-03
43–7	Chr3:58 095 412	rs756771275	Heterozygotes	c.C2309T	p.T770I	4.00E-04	SOX6, TTC26

*Nucleotide change is based on FLNB isoform with accession of NM_001457.3.

†Maximum frequency across public databases: 1000 Genomes Project phase 3, Genome Aggregation Database and National Heart, Lung and Blood Institute Exome Sequencing Project databases.

AIS, adolescent idiopathic scoliosis; N, absence in any database.

We then tested whether rare damaging variants in genes within the protein–protein interacting network could occur in patients with AIS. The rare variant burden among the interacting partners of the three aforementioned significant AIS genes, which were FLNB, TTC26and PTK7, was analysed. Among the 46 interactome genes tested, FLNA, encoding an interaction partner of FLNB, was also significantly enriched with rare damaging variants in patients with AIS (p=2.58E-03, online supplementary table 6).

Functional characterisation of AIS-associated FLNB variants

According to alignment to the FLNB protein domain, most of the AIS-associated FLNB variants are located within immunoglobulin-like filamin repeat regions, some of which belong to the domain of interaction with FLNA15 (figure 2A). Of note, p.R199Q is located within the actin-binding domain of FLNB. We transfected wild-type or mutant plasmids into HEK293T cells and found that some FLNB variants (including p.M1803L, p.S2503G and p.T2166M; online supplementary figure 1) resulted in cytoplasmic focal accumulation, and some other FLNB variants (including p.R566L, p.A2282T, p.S2503G, p.R199Q and p.R2003H; online supplementary figure 2) altered actin dynamics (online supplementary figures 1 and 2).

Figure 2

FLNB and TTC26 variants in individuals with AIS. (A) Profiles of rare damaging variants in FLNB. Rare variants are represented by lollipops and counts of alleles with variants in cases (top panel) and controls (bottom panel) are shown. (B, C) Local view of in silico structure analysis for the WT and mutant FLNB structures (B, variant R566L; C, variant A2282T). The WT structure of FLNB is shown in purple, and the mutant structure of FLNB is shown in green. The side chains of R/L566 and A/T2282 are shown as sticks, and the other residues are shown as lines. (D, E) A total of 293 T-cells were transfected with Flag-tagged WT or mutant FLNB (p.R566L, p.A2282T) vector plasmids and myc-tagged WT or mutant TTC26 (p.R297C, p.R50C). Then, communoprecipitation assays were conducted. Western blot images are representative of n=3 experiments. AIS, adolescent idiopathic scoliosis; WT, wild type.

To investigate the protein–protein interactions, we focused on AIS trios with multiple variants. We found that patients in two AIS trios (trios 22 and 27) carried variants in both the FLNB and TTC26 genes (figure 1). We first performed in silico analyses to investigate the potential impact of FLNB variants on protein conformation. Figure 2B, C indicates that there is a large conformational change between the wild-type and mutant p.A2282T and p.R566L FLNB protein structures. The side chains of the residues in the loops of T2282 move in the opposite direction to those of A2282. In the secondary protein structure, the side chain of R566 can form hydrogen bonds with S567, A568, D569 and G561; however, the smaller side chain of the L566 mutant cannot form any hydrogen bonds with the other residues. The distance between the two loops is larger than in the wild-type structure. It has been reported that FLNB localises to the basal body and the proximal regions of the cilium, a non-motile microtubule-based organelle that projects from the cell surface. Since TTC26 is an intraflagellar transport (IFT) protein in cilia,16 we aimed to identify potential interactions between FLNB and TTC26. Using coimmunoprecipitation assays, we found that the myc-tagged mutant p.R50C and p.R197C TTC26 proteins pulled down the Flag-tagged mutant p.A2282T and p.R566L FLNB proteins, respectively (figure 2D, E).

We were also interested in case 98–73, whose twin sister was also diagnosed with AIS (figure 3A). The FLNB missense variant p.R2003H is located in a highly conserved region of the FLNB protein (figure 3B). In silico analyses (figure 3C) indicated that the R2003 residue was solvent accessible and was positioned far from the β-sheet secondary structure. In the FLNB mutant protein structure, the side-chain of H2003 formed a strong hydrogen bond with E2078. Since OFD1 localises to the base of the cilium,17 we assumed that FLNB may interact with OFD1. Coimmunoprecipitation analysis indicated an interaction between wild-type OFD1 and wild-type FLNB, which did not exist between p.R2003H FLNB and p.Y437F OFD1 (figure 3D).

Figure 3

FLNB and OFD1 variants in individuals with AIS. (A) Pedigree of AIS twins. Case 98–73 (proband) is indicated with an arrow. (B) Protein sequences around FLNB. p.R2003 in 11 species. (C) Local view of in silico structure analysis of the WT and mutant FLNB structures (variant H2003). The WT structure of FLNB is shown in purple, and the mutant structure of FLNB is shown in green. The side chains of R/H2003 are shown as sticks, and the other residues are shown as lines. (D) A total of 293 T-cells were transfected with Flag-tagged WT or mutant FLNB (p.R2003H) vector plasmids and myc-tagged WT or mutant OFD1 (p.Y437F). Then, coimmunoprecipitation assays were conducted. Western blot images are representative of n=3 experiments. AIS, adolescent idiopathic scoliosis; WT, wild type.

Clinical significance of oligogenic variants of AIS-associated genes

To further investigate the oligogenic nature of AIS, we investigated the clinical characteristics of patients harbouring different numbers of variants of AIS-associated genes. We used the development of a curvature of at least 30° before any treatment as the outcome variable, since a curvature greater than 30° after skeletal maturity was used to define curve progression.18 Data on progression-free survival were available for all 223 patients. The median follow-up period was 143 weeks (range 17–372 weeks). A total of 31.8% (71/223) of the patients exhibited progressed curves. Univariate survival analysis revealed that curve magnitude at first presentation (>23 vs ≤23, p=2.00E-03, HR=2.137, 95% CI 1.324 to 3.448) and the number of rare damaging variants (≥2 vs 0 or 1, p=7.69E-11, HR=5.098, 95% CI 3.121 to 8.328) were of prognostic significance (table 4). In the multivariate analysis, both curve magnitude at presentation (>23 vs ≤23, p=2.50E-02, HR=1.781, 95% CI 1.074 to 2.955) and the number of rare damaging variants (≥2 vs 0 or 1, p=3.29E-07, HR=4.304, 95% CI 2.458 to 7.537) proved to be independent prognostic factors associated with curve progression.

Table 4

Risk factors and their effect on curve progression by univariate Cox proportional hazard regression analyses

Characteristics	Patients (n)	HR (95% CI)	P value
Sex
Male	24	1
Female	199	1.197 (0.593 to 2.415)	6.16E-01
Age of curve onset (years)
≤12	127	1
>12	96	0.661 (0.406 to 1.077)	9.70E-02
Curve magnitude at first presentation (°)
≤23	130	1
>23	93	2.137 (1.324 to 3.448)	2.00E-03
Risser sign
≤2	123	1
>2	100	0.871 (0.545 to 1.393)	5.65E-01
Body mass index (kg/m2)
≤21.96	86	1
>21.96	137	0.680 (0.426 to 1.084)	1.05E-01
Lenke classification
1	23	1
2	110	1.290 (0.503 to 3.308)	5.96E-01
3	28	1.274 (0.403 to 4.025)	6.80E-01
4	27	2.562 (0.906 to 7.244)	7.60E-02
5	29	1.889 (0.640 to 5.572)	2.49E-01
6	6	1.992 (0.384 to 10.372)	4.12E-01
Number of rare damaging variants
≤1	186	1
>1	37	5.098 (3.121 to 8.328)	7.69E-11

Discussion

In this study, we investigated an oligogenic model for AIS using exome sequencing. Although numerous loci and candidate genes have been revealed to be associated with AIS, no single gene has been established to cause AIS.11 The relatively sporadic characteristics of this disease and its variable phenotypes (scoliosis severity and curve type), even within families in which multiple individuals are affected, suggest that AIS is a multigenic and multifactorial disease.10 Here, we screened the exome data for variants of AIS-associated genes and demonstrated more frequent involvement of oligogenicity in patients with AIS. More importantly, the number of rare damaging variants was recognised as an independent prognostic factor for curve progression, indicating that a higher genetic load might lead to a severe phenotype.

The top gene prioritised by the gene-based burden test was FLNB. FLNB, encoded by FLNB, interacts with the cytoskeleton by linking the actin network with cellular membranes and mediates interactions between actin and transmembrane receptors. FLNB is expressed in human growth plate chondrocytes and in the developing vertebral bodies of the mouse.19 Different variants in the human FLNB gene have been found in a series of skeletal disorders,15 19 with null alleles of FLNB resulting in recessive spondylocarpotarsal syndrome and missense variants or small in-frame deletions or insertions in FLNB causing a group of autosomal-dominant diseases, including atelosteogenesis (AO) I and atelosteogenesis III, boomerang dysplasia (BD), and Larsen syndrome (LS).15 It has been reported that AO–BD–LS-associated variants in the FLNB gene can cause cytoplasmic focal accumulation of FLNB, which may affect the mechanotransduction generated through integrin–filamin interactions that ultimately influence chondrocyte and osteoblast proliferation and differentiation.15 20

Here, we report that rare damaging variants in the FLNB gene might contribute to AIS. Notably, most of the FLNB variants identified in our study are located within the immunoglobulin-like filamin repeat domain, while the variants leading to AO–BD–LS are distributed mainly within the actin-binding domain and hinge region 1,15 indicating that variants in different domains of FLNB may cause various phenotypes through different mechanisms, since no skeletal deformations other than scoliosis were found in our patients with AIS. In the AIS families, patients harbouring variants in FLNB also carried variants in other AIS-associated genes, which indicated that FLNB might play a role as a modifier in AIS.

The second-highest ranked gene was TTC26, which encodes the cilia protein TTC26 (IFT56/DYF‐13). Knockdown of TTC26 in zebrafish embryos causes ciliary defects, including curly body axis, left-right asymmetry defects and hydrocephalus.21 In our study, we found that two patients from AIS families harboured variants in both the TTC26 and FLNB genes. Surprisingly, the variants found in the AIS families result in a stronger interaction between FLNB and TTC26. TTC26, as an IFT complex B protein, is required for the accumulation of other core IFT-B proteins, tubulin transport and Shh signalling in cilia.22 Our data indicated that the interaction of mutant FLNB and TTC26 proteins might influence the role of TTC26 in maintaining IFT-B complex integrity and cilia function, thus affecting the mechanosensory properties or differentiation of chondrocytes and osteoblasts.23

OFD1, variants of which were identified in an AIS twin family, also encodes a cilia-related protein. Variants of OFD1 lead to X-linked oral–facial–digital syndrome24 and Joubert syndrome.25 Ofd 1 antisense morpholinos injected zebrafish embryos exhibit marked body axis curvature, hydrocephalus and other malformations, which are associated with dysregulated ciliary function and Wnt signalling.26 The OFD1 protein is composed of five helical coiled-coil domains and an N-terminally localised LisH domain. It interacts via its fifth coiled-coil motif with PCM-1 and via its C-terminal non-coiled-coil region with BBS4.27 Our observation that OFD1 interacts via its third coiled-coil motif with FLNB might help to explain the phenotypic variety of OFD1 variants.

Rare variants in FBN1 and FBN2 have been reported to be enriched in severely affected AIS cases in individuals of both European and Han Chinese ancestries.4 We identified rare damaging variants of FBN1 or FBN2 in 16/223 (7.17%) patients with AIS, including the full range of curve severities, which was a similar rate to the reported 7.6% of Han Chinese cases with rare coding variants in FBN1 or FBN2.4 However, none of these variants were shared by the IS cohort reported by Buchan et al.4 FBN1 was ranked as the top signal in Buchan’s exome-wide association study but not in our burden test study. This discordance might be related to two factors. First, in Buchan’s study, the analysis was restricted to patients with AIS with severe scoliosis (spinal curves measuring ≥40° or surgically treated), while in this study, we included patients with broader representation of curve severity. Second, in Buchan’ study, all cases and controls were of European ancestry in the gene-burden analysis. In our study, we focused on patients and controls of Chinese ancestry, which may play an important role in the identification of different candidate genes. Interestingly, 8/16 (50%) of the patients harbouring rare damaging variants of FBN1 or FBN2 exhibited an additional variant of an AIS-associated gene, which may help explain the incomplete penetrance of rare FNB1 and FBN2 variants in some AIS families and supports the oligogenic model of AIS.4

Cilia are microtubule-based organelles that extend from the surface of almost all vertebrate cells. Motile cilia beat in a rhythmic fashion to generate fluid flow across surfaces, whereas immotile primary cilia typically function as specialised sensory structures.28 29 Defects in the cilium-driven cerebrospinal fluid flow have been indicated as a contributory factor in IS.7 30 31 Recently, elongated primary cilia of the osteoblasts of patients with IS were reported.32 Considering the emerging evidence indicating that primary cilia act as mechanosensors in bone structures,33 34 dysregulated ciliogenesis and cilia function-associated cellular-molecular mechanisms might provide another explanation for the IS aetiology.35 In our study, we demonstrate that some ciliary genes might contribute to AIS with oligogenic inheritance.

In a gene-based burden test for rare variant associations, we combined data from in-house controls and 1000 Genomes Project controls and found that rare damaging variants in several AIS-associated genes (FLNB，TTC26 and CNTNAP2) were over-represented in patients with AIS. The results were quite similar when we compared the data between the patients with AIS and our sex-matched in-house controls (online supplementary table 7). AIS is a sexually dimorphic disease in which women exhibit a greater risk of progressive curves.36 One study involving a cohort of multiplex AIS families demonstrated that a greater genetic load is required for men to be affected.11 Several sexually dimorphic IS susceptibility loci have been identified,37 38 which may help to explain the genetic influences on the development of AIS in women. The similar results obtained from gene-based burden tests when we combined the in-house and external controls or used only the in-house controls may suggest that a higher genetic load contributes to more progressive curves in the context of AIS without sexual dimorphism. In addition to genetic factors, the correlation between environmental factors, including hormonal dysregulation, and the severity of spinal curvature has also been investigated, although the role of hormones continues to remain controversial.39 40 Additional larger studies are needed to fully understand the pathophysiological mechanism of curve progression in AIS.

In this study, we focused on rare damaging single-nucleotide variants. Copy number variants (CNVs) can also influence the risk of many diseases and are associated with several skeletal phenotypes.41–43 In one study, more than 6% of patients with AIS were found to harbour a clinically important copy number abnormality, and some of these abnormalities may play a role in AIS pathogenesis.44 In another comprehensive, large-scale study, distal chromosome 16p11.2 duplications were identified in approximately 1% of patents with AIS, who also presented an OR that was much higher than any other OR previously reported for AIS.45 The enrichment of CNVs is intriguing, and further studies examining the combination of single-nucleotide variants and CNVs may help to better elucidate the oligogenic model and genetic basis of AIS.