Evaluation of genetic susceptibility of common variants in SOX9 in patients with congenital talipes equinovarus in the Han Chinese population

Background Congenital talipes equinovarus (CTEV) is a common birth defect that causes severe deformities of one or both feet. Genetics have been proven to play a key role in the risk of CTEV. Our study aimed to evaluate the genetic susceptibility of common variants in the SOX9 gene to CTEV in a Han Chinese population. Methods In this study, we recruited 2,205 study participants, including 692 CTEV patients and 1513 healthy controls. A total of seven selected single-nucleotide polymorphisms (SNPs) within the SOX9 gene were genotyped, and environmental variables, including maternal smoking and alcoholic drinking habits, were assessed. In addition, bioinformatics analyses were performed to explore the potential biological functions of the associated SNPs. Results The SNP rs73354570 was identified to be significantly associated with the risk of CTEV (OR = 1.53, P = 2.11 × 10−5), and the C allele was associated with an increased risk of CTEV. A dose-dependent pattern could be observed in genotypic analyses. The OR for individuals with AC genotypes was 1.37 (95% CI 1.09–1.71), and the OR for individuals with CC homozygotes was 1.47 (95% CI 1.18–1.82). Further analyses identified that rs73354570 is located within a region of multiple binding proteins, including CEBPB and POLR2A, which suggested that this SNP was also part of genetic motifs that are found within several cell types. Conclusion Our results provide evidence supporting the important role of the SOX9 gene in the contribution to the risk of CTEV.

SOX9 encodes a transcription factor that can recognize and bind to the sequence of CCTTGAG. Early studies have shown that the products of SOX9 can bind to COL9A1 and regulate its transcription [8]. COL9A1 encodes an α chain of type IX collagen. Early studies have shown that mice with a knockout of the Col9a1 gene do not produce any collagen IX polypeptides, which indicates that COL9A1 is essential for the synthesis of type IX collagen. Recent studies have shown that variations in collagen IX-related genes can cause alterations in bone marrow hyperplasia and might therefore be related to articular cartilage-related diseases [8,9]. Zhou et al. investigated the association between common SNPs of COL9A1 and clubfoot in a small sample of individuals with Chinese Han ancestry, and a potentially related SNP was identified [10]. Based on these multiple lines of evidence connecting SOX9 and COL9A1, we hypothesize that SOX9 might be involved in the regulation of articular cartilage development. In a recent study, Wang et al. detected significant differences in the mRNA and protein expression levels of SOX9 between idiopathic CTEV muscle samples and controls [11]. Nevertheless, no evidence has been identified to connect the genetic polymorphisms of SOX9 and the risk of CTEV in the human population.
In this study, we aimed to investigate the potential genetic association between genetic polymorphisms of SOX9 and the risk of CTEV based on thousands of study participants with Chinese Han ancestry. In addition to single marker-based association analyses, we also examined the effects of some environmental factors, including maternal smoking and alcoholic drinking habits. Geneby-environment interactions were also explored.

Study participants
We recruited 692 CTEV patients and 1513 unrelated healthy controls from Honghui Hospital of Xi'an Jiaotong University and Xi'an Children Hospital from January 2012 to April 2017. All participants included in the study were randomly chosen genetically unrelated Han Chinese individuals. All patient diagnoses were confirmed by X-ray examinations. All controls without foot deformities were matched with patients by age and gender. The demographic and clinical characteristics of all participants were obtained from questionnaires and medical records. This study was performed in accordance with the ethical guidelines of the Helsinki Declaration of 1975 (revised in 2008) and was approved through the Ethics Committee of Honghui Hospital of Xi'an Jiaotong University. Written informed consent was obtained from all participants.

SNP selection and genotyping
Tagged SNPs located within gene regions with minor allele frequency (MAF) > 0.01 in SOX9 in 1000 Chinese Han genomes were chosen for genotyping. Algorithm Tagger integrated in Haploview [12] was used for SNP tagging, and the r 2 criterion used for tagging was 0.8 for both gene regions. The r 2 criterion was used for tagging. A total of 7 candidate SNPs were selected for genotyping in this study.
Genomic DNA was isolated from peripheral blood leukocytes according to the manufacturer's protocol (Genomic DNA kit, Axygen Scientific, Inc., CA, USA). SNP genotyping was performed using the high-throughput Sequenom MassARRAY platform with iPLEX GOLD chemistry (Sequenom, San Diego, CA, USA) based on the manufacturer's protocols. The results were processed using Sequenom Typer 4.0 software, and genotype data were generated from the samples [13]. The case and control sample results were blinded for quality control during genotyping processes [14], and 5% of samples were randomly processed with a concordance of 100%.

Statistical methods
Genetic association analyses were performed at both genotypic and allelic levels using Plink v1.9 [15]. A linkage disequilibrium (LD) plot was made using Haploview [12]. In addition to ordinary genetic association analyses, we conducted gene-by-environment interaction analyses by logistic models. For gene-by-environment analyses, two environmental factors, i.e., maternal smoking and maternal alcoholic drinking, were included. The smoking and drinking status were recorded as their habit in general but not status during their pregnancy. Each of these seven selected SNPs was examined by pairing with each of the two environmental factors. Both maternal smoking and alcoholic drinking were coded as 0, 1, and 2, indicating never, occasionally, and often, respectively. Age and gender were included in the logistic models described above. Bonferroni corrections were applied to address multiple comparisons. In addition to statistical methods, we conducted bioinformatics analyses to explore the potential biological functions of our targeted SNPs. RegulomeDB, which annotates SNPs using ENCODE data, was used to investigate the potential regulatory roles of significant SNPs [16]. We also examined the potential eQTL patterns of those significant SNPs using GTEx [17]. Moreover, we investigated the gene-gene network of SOX9 using the STRING database, which is a database of known and predicted protein-protein interactions.

Demographic and clinical characteristics of our study participants
Of the patients included in this study, 8.1% had a family history of CTEV compared with 1.4% of the controls.
There was a significant difference in this variable between the study groups (P < 0.001, Table 1). However, no significant differences were found in age, gender, maternal smoking, or maternal drinking between groups (Table 1).

Genetic associations signal
In our study, all seven SNPs selected were in Hardy-Weinberg equilibrium. Basic information, including MAF and the results of Hardy-Weinberg equilibrium tests for these 7 SNPs, is included in Table 2. Weak LD could be identified among these 7 SNPs (Fig. 1). Rs73354570 of SOX9 was the only SNP identified to be significantly associated (odds ratio = 1.53, P = 2.11 × 10 −5 ) with the disease status of CTEV ( Table 3). The C allele was associated with an increased risk of CTEV. A dose-dependent pattern could be observed in genotypic analyses. The ORs for individuals with AC genotypes and CC homozygotes were 1.37 (with 95% confidence interval 1.09-1.71) and 1.47 (with 95% confidence interval 1.18-1.82), respectively. The only interaction signal that achieved nominal significance was identified between SNP rs73354570 and maternal alcoholic drinking habit (P = 0.02). However, this signal was no longer significant if we correct for multiple comparisons (threshold of P value is 0.007). The full results are summarized in Table 4.

Bioinformatics analyses
We explored the potential functional significance of SNP rs73354570 in RegulomeDB and GTEx. The Regulo-meDB score for rs73354570 was 2b. The scoring system of RegulomeDB ranges from 1 to 6, and a lower score often indicates that the SNP has a more significant role in biological functions. A further examination of this SNP identified that rs73354570 is located within a region of multiple binding proteins, including CEBPB and POLR2A. This SNP was also part of genetic motifs found within several cell types. We did not identify significant eQTL signals from GTEx for SNP rs73354570 on SOX9 (Supplemental Table S1). The gene-gene network of SOX9 is shown in Fig. 2. According to the STRI NG database, the protein product of SOX9 interacts experimentally with several proteins encoded by the PRKG2, COL2A1, RUNX2, AMH, NR5A1, FOXL2, and CTNNB1 genes. In addition, SOX9 was also predicted to be connected with other proteins encoded by genes, including COL10A1, ACAN, and SHH.

Discussion
In this study, we identified an SNP, i.e., rs73354570, that was significantly associated with the disease status of CTEV in a large Chinese Han population-based sample. To the best of our knowledge, our study is the first to report a significant association between SOX9 and CTEV. SOX9 was not reported to be significant in the first and only previous GWAS that focused on clubfoot, which was performed by Zhang et al. [7]. A GWAS with a stringent P value threshold (5 × 10 −8 ) might miss some true positive hits. Our candidate gene-based study design could avoid this limitation by testing fewer markers of preselected genes. To date, the genes SOX9 and SOX9-AS1 (gene encoding anti-sense RNA1 of SOX9) have been identified to be associated with several traits and disorders, including angiotensin-converting enzyme inhibitor intolerance [18], height [19], liver enzyme levels (gamma-glutamyl transferase) [20], lung function [21], nose morphology [22], and thyroid hormone levels [23]. However, no GWAS has established the connection between SOX9 and CTEV-related disorders.
Our bioinformatics analyses have identified that SNP rs73354570 might have significant biological functions in the gene expression of SOX9. This SNP would likely affect the protein binding of its surrounding areas and thus might affect the gene expression of SOX9. However, this information extracted from RegulomeDB could not be replicated by GTEx eQTL data. No significant differences could be found in the gene expression of SOX9 in different genotype groups of SNP rs73354570 in multiple human tissues. Thus, the eQTL data extracted from  [11]. In the present study, we identified the C allele of rs73354570 to be significantly associated with an increased risk of CTEV. In the future, an eQTL study for rs73354570 on the gene expression of SOX9 based on the target tissue of CTEV patients should be conducted to evaluate the allelic effects on gene expression. With the rapid development of sequencing technology, numerous susceptibility loci contributing to complex diseases have been reported, such as schizophrenia [24][25][26]. Considering that analyses of only some SNPs are not sufficient to draw conclusions [27][28][29][30][31], we conducted gene-byenvironment interaction analyses to investigate potential interactions between our selected SNPs and two environmental factors along with genetic association analyses. However, no significant interaction signals were identified. Nevertheless, we believe that it is not necessary to overinterpret this negative result. Since smoking and alcoholic drinking have been proven to be significantly associated with multiple birth defects [32,33], it is probable that multiple environmental factors could combine with genetic factors to play a role in the process of CTEV onset. Notably, the most significant interaction signal was obtained from the significant hit of association analysis (rs73354570). Although this signal failed to achieve genome-wide significance, it is still worth investigating this interaction signal further in the future because it might partly explain its single marker-based association signal. A potential limitation of the present study is that smoking and drinking status were recorded as habits in general but not status during pregnancy. Additional studies with appropriately measured environmental factors should be conducted to further investigate the combined effects of both environmental factors and genetic factors on CTEV.
Our study has several limitations. The most important factor is the lack of replication. Future studies are needed to replicate our findings regarding SOX9 in Chinese Han and other populations. Moreover, we did not perform any procedures to adjust for population stratification in the study. However, we tried our best to at least partially control this potential confounding factor because geographic location is a good indicator for genetic matching in the Han Chinese population [34,35]. Another limitation is that we only included SNPs located within gene regions. However, several important regulatory regions are located at up-or downstream regions that are not within gene   regions, and a significant portion of the GWAS panel markers cannot be mapped to any gene regions (but are several kb away from the targeted genes). It is very difficult to claim that seven preselected SNPs could represent most of the genetic information of SOX9. In addition, as a candidate gene-based study focusing on the effects of common SNPs, we did not examine any rare or lowfrequency variants. However, a recent study showed that low-frequency and rare variants might play an important role in the onset and development of CTEV [36]. Based on Chinese family samples with CTEV, Zhang et al. identified two pathogenic variations from mediator complex subunit 13L (MED13L) and transforming growth factor-β receptor 2 (TGFBR2). Therefore, a sequencing-based study might provide more information about the genetic etiology of CTEV.
In this study, we identified potential links between genetic polymorphisms of SOX9 and the risk of CTEV. Our study could improve our understanding of the genetic architecture of CTEV and provide a basis for novel intervention plans. Replication studies involving Chinese Han and other populations are still needed in the future.