Significant Associations between AXIN1 rs1805105, rs12921862, rs370681 Haplotypes and Variant Genotypes of AXIN2 rs2240308 with Risk of Congenital Heart Defects

This study aimed to investigate possible associations of the susceptibility to congenital heart defects (CHDs) with AXIN1 rs1805105, rs12921862 and rs370681 gene variants and haplotypes, and AXIN2 rs2240308 gene variant. Significant associations were identified for AXIN1 rs370681 and AXIN2 rs2240308 variants. AXIN1 rs370681 variant was significantly associated with decreased odds of CHDs (adjusted OR varying from 0.13 to 0.28 in codominant, dominant and recessive gene models), while the AXIN2 rs2240308 variant was associated with increased odds of CHD in the dominant model. The haplotype-based generalized linear model regression of AXIN1 rs1805105, rs12921862 and rs370681 variants revealed that C-C-C and C-C-T haplotypes significantly increased the risk of CHDs (p < 0.05). No significant second order epistatic interactions were found between investigated variants (AXIN1 rs1805105, rs12921862, rs370681, and AXIN2 rs2240308). Our conclusion is that AXIN1 rs1805105, rs12921862, and rs370681 (C-C-C and C-C-T) haplotypes and AXIN2 rs2240308 contribute to CHDs susceptibility.


Introduction
Congenital heart defects (CHDs) are the most common birth anomaly and occur in approximately 1% of live births [1]. Epidemiological studies revealed that genetic or environmental factors were identified in 25% of cases with CHDs [2]. Recently, one of the largest genetic studies in congenital heart diseases established that 8% of cases had "de novo" autosomal dominant variants, and 2% were due to autosomal recessive inheritance [3]. From a clinical point of view, genetic causes in CHDs are important to be identified for the following reasons: to anticipate the patients' evolution and their

Materials and Methods
This case-control study included 214 Caucasian subjects divided into two groups. The control group consisting of 111 healthy children while the patient's group consisted of 103 CHDs patients from different regions of Romania admitted in the Clinic of Pediatric Cardiology, Emergency Institute for Cardiovascular Diseases and Transplantation Targu Mures, Romania. In the patients group, children with Atrial septal defect (ASD) (51.5%), Ventricular septal defect (VSD) (12%), Tetralogy of Fallot (TOF) (15%), or Double-outlet right ventricle (DORV) (3%) and patients with the combined ASD and VSD (18.5%) were included. The family history was negative in all cases, and the questionnaire completed by the mothers did not highlight prenatal exposure to environmental risk factors related to CHDs (such as maternal medication, alcohol or smoking exposure, parental consanguinity). Patients with a clinical phenotype suggestive for a syndrome (such as Marfan syndrome, Turner syndrome, etc.) or with other congenital malformations were not included in this study. The control group included voluntary non-hospitalized children from different region of Romania, with normal echocardiogram examination, and hospitalized children referred to Pediatrics department for different acute conditions without any congenital structural malformations.
Clinical and demographical data of children included in the study were collected using clinical examinations and medical records (including transthoracic echocardiographic examination) and a questionnaire for the parents or legal guardian of the child.
Ethical approvals were obtained from the Ethical Committee of George Emil Palade University of Medicine, Pharmacy, Science, and Technology of Targu Mures, Romania (47/23.02.2018), and from the Emergency Institute for Cardiovascular Diseases and Transplant Targu Mures. The applied ethical principles were according to the Declaration of Helsinki. The subjects were included in the research activity after a consent, signed by their parents, was obtained.

Genotype Investigation
The DNA was obtained from buccal cells collected with oral swabs (Isohelix, SK-2S, Cell Projects, Kent, UK). For DNA extraction, we used the PureLink Genomic DNA kit (Invitrogen, Carlsbad, CA, USA) and the recommendations of the manufacturer. DNA concentration and absorbance were measured by spectrophotometry (Eppendorf BioSpectrometer, Eppendorf AG, Hamburg, Germany).
For genotyping AXIN1 rs370681 and AXIN2 rs2240308, we used Real-Time Polymerase Chain Reaction (PCR) on 7500 Fast Dx system (Applied Biosystem, Foster City, CA, USA) with TaqMan assays (Thermo Fisher Scientific, Waltham, MA, USA), namely C___2221642_20 and C___2577354_1_. Genotyping of AXIN1 rs1805105 and rs12921862 was performed using polymerase chain reaction-restriction fragment length polymorphism (PCR-RFLP) using specific primers and restriction enzymes, as previously reported [18].
Ensemble genome browser and the Variant Effect Prediction tool of this database were used to annotate and describe the variants.

Data Analysis
Quantitative baseline patients' characteristics are presented as median with interquartile range (IQR) while qualitative characteristics were represented as absolute and relative frequencies.
Mann-Whitney and Chi-square tests were performed for nonparametric data.
Statistical analysis was made in R software v.4.0.0 [19]. Each variant frequency genotype distribution was assessed for departure from Hardy-Weinberg Equilibrium (HWE) using an exact Chi-square test among CHD cases and controls. The exact Chi-square test from "genetics" R package was used [20].
The associations between the studied gene variants and odds of CHD were tested by unconditional logistic regression. According to our knowledge, the age and gender of children were not considered as potential factors for CHD risk in the pediatric population and as a consequence, we did not perform adjustment for age and gender, and for each variant, the odds ratio (OR) and 95% confidence interval (CI) without adjustment for age and gender were estimated. We used four inheritance gene models (codominant, dominant, recessive, and overdominant) to assess the associations between variants and odds of CHD using "SNPassoc" R package [21]. The same package was used to test the potential pairwise interaction between the investigated variants.
Haplotype estimation was performed using R package haplo.stats. The most common haplotype with a significant negative association with the disease risk served as the reference/baseline category. A generalized linear model (GLM) of the haplo.stats package without an interaction term was used.
All statistical tests were two-sided, the threshold of statistical significance being set at α = 0.05. In the case of four inheritance gene models used in binomial logistic regression, the crude p-values were adjusted for multiple testing via Benjamini-Hochberg method.

Results
The ages range from 1 month to 17.5 years in both controls and CHD patients. The median age of CHD patients and control subjects was 6 years (interquartile range [IQR]: 2.0-11.1) and 3 years (interquartile range [IQR]: 0.5-9.0), respectively (p = 0.001).
Concerning prenatal risk factors for congenital heart disease, there was a significant mean difference in maternal age at conception between CHD and control groups (p = 0.012); the mean (SD) of maternal age for CHD patients was 26.9 ± 5.6 years and 25.1 ± 4.9 years for control group. The investigated mothers disclaimed the maternal medication, alcohol or smoking exposure during pregnancy. Similarly, type 1 diabetes, lupus and exposure to nonfertility medications during pregnancy were not observed except for five cases with maternal urinary tract infections that were treated with natural drugs.
The genotypes distributions in CHD and control groups under four inheritance gene models: codominant, dominant, recessive, and overdominant is presented in Table 1. In addition, we explored the difference in genotype distribution in controls and in ASD patient group (with or without VSD). As it may be observed, ASD was the most common CHD found in the present study (72 cases, 69.9%).
In addition, the AXIN2 rs2240308 gene variant was a significant risk factor for ASD development under the allelic model (OR = 1.60, 95% CI: 1.05-2.44), while AXIN1 rs12921862 and AXIN1rs370681 gene variants were protective factors ( Table 2).  Note. a allele frequency in worldwide populations and European (non-Finnish) population from Genome Aggregation Database (gnomAD; https://gnomad.broadinstitute.org); 95% CI = 95% confidence interval; b -unadjusted odds ratio in the allelic model; Significant results were achieved when p-values < 0.05 and the corresponding OR and CI values are highlighted with bold.

AXIN1 rs1805105, rs12921862, and rs370681 Haplotypes and Their Associations with CHD Risk
The haplotype analysis of AXIN1 rs1805105, rs12921862, and rs370681, revealed eight haplotypes. After eliminating one haplotype with an estimated frequency smaller than 0.021, seven haplotypes were used in the haplotype-based association analyses. Haplotypes distribution and their associations with CHDs risk are illustrated in Table 3. The frequency of haplotype C-T-C was lower in CHD cases compared to controls (25.8% versus 33.5%, p = 0.026). To identify possible risk haplotypes and quantify the effect of each haplotype, we performed a haplotype-based GLM regression, in which haplotype C-T-C was chosen as reference. The results of model fitting showed that C-C-C and C-C-T haplotypes remained significantly associated with increased CHD risk after adjusting for maternal age (Table 3). The same haplotypes (C-C-C and C-C-T) had a significant positive association with ASD risk after accounting for maternal age, playing the role of risk factors with an increased risk varying from 4.97 for C-C-C versus C-T-C haplotype and 3.90 for C-C-T versus reference haplotype (Table 4).

Epistatic Pairwise
Interactions between AXIN1 rs1805105, rs12921862, rs370681 and AXIN2 rs2240308 Epistatic pairwise interactions between the investigated variants are illustrated in Table 5. No significant second order epistatic interactions were found. Note. The p-value of epistatic pairwise interactions represented in the upper part of the matrix was obtained using the log-likelihood ratio test (LRT). The p-value from the diagonal of the matrix represent unadjusted (crude) effect on CHD of all investigated variants, obtained from LRT. Also, the p-value from the lower part of the matrix was obtained using LRT comparing the likelihood of the model for variants of two genes and the best model containing one gene variant.

Discussion
AXIN protein and Wnt/β-catenin canonical mechanisms were described as being involved in cardiovascular diseases; disfunction of AXIN protein and Wnt/β-catenin canonical mechanism may cause congenital heart disease [22][23][24]. Wnt pathway plays an important role in the early stage of heart development [25], and inhibition of the β-catenin gene in the endodermal development stage leads to abnormal cardiac development [26].
Similarities and interactions between AXIN1 and AXIN2 genes and between their codified proteins were reported previously [13,27]. Variants of these genes may increase the level of AXIN protein with a normal level of mRNA, suggesting that elevated levels of protein may be explained by the increased stability of AXIN protein [28].
To the best of our knowledge, the present study is the first one that evaluates four variants of the AXIN gene (namely AXIN1 rs1805105, rs12921862, and rs370681 and AXIN2 rs2240308) at the same time and to investigate their interactions.
Our study showed that AXIN1 rs1805105 and rs2921862 variants were not associated with susceptibility to CHDs. A significantly negative association was identified for AXIN1 rs370681 variant, suggesting a potential protective effect of the variant (T) allele. Contrary to our results, Pui Y et al. [18] identified a significant association between variant genotype of the AXIN1 rs12921862 and ASD susceptibility, with an approximately 3-fold increased risk [18]. Moreover, two of the investigated variants in the AXIN1 gene, rs12921862 and rs1805105, were reported by Kai L et al. [15] as being associated with an increased risk for dilated cardiomyopathy, in Chinese Han population [15]. The different results may be explained by the multiple types of CHDs (ASD, VSD, TOF, and DORV) included in our study. Other causes may be the small number of cases and different ethnicity of the subjects, with differences being reported in allele frequency, according with gnomAD database ( Table 2) across all human populations. Moreover, in our investigated population (pediatric patients and controls), only the frequency of homozygous genotype for AXIN1 rs12921862 was similar with the data reported in previous studies [15,18]; the frequency of heterozygous genotype was higher in our investigated population, leading to unsignificant results compared to Chinese Population (23.1% in ASD patients and 8.3% in controls reported by Yan Pu et al. [18], 37.1% in dilated cardiomyopathy patients and 8.6% in controls reported by Kai Li et al. [15], compared to 30.1% in our CHD patients and 40.5% in our controls).
Regarding the AXIN2 rs2240308 variant, our study showed a significant association between the variant genotypes and CHDs risk in the dominant model. According to the latest published data, our study is the first one who analyzed the AXIN2 rs2240308 variant on CHDs patients. Previous studies evaluated the mentioned variant on patients with hypodontia and demonstrated that the presence of a variant allele was associated with frontal agenesis [29]. Also, the mentioned variant was previously associated with cancer risk, such as colorectal cancer [30], breast cancer [31], and Hirschsprung disease [32]. In the case of CHDs patients, other AXIN2 variants [c.28 C > T (p.L10F), c.395 A > G (p.K132R)] were previously reported by Zhu M et al. [33] as pathogenic variants for CHDs.
The haplotype analysis for AXIN1 variants performed in the current study identified significant associations between CHDs and C-C-C and C-C-T haplotypes, the presence of the mentioned haplotypes being risk factors for CHD development. The difference between the mentioned haplotypes is the latest allele, the allele of AXIN1 rs1805105 variant. Moreover, our data showed that C-C-C and C-C-T haplotypes had a significant positive association with ASD risk. The OR in the case of the AXIN1 rs1805105 C allele (OR = 3.49, 95% CI = 1.12-10.87) is higher compared to that for AXIN1 rs1805105 T allele (OR = 2.88, 95% CI = 1.51-5.51). Separately, in the univariate analysis, AXIN1 rs1805105 variant was not associated with CHD risk but our results for haplotype analysis suggest that in the case of C-C-C and C-C-T haplotypes, the risk effect is modulated by AXIN1 rs1805105 variant, with a higher risk for AXIN1 rs1805105 C allele.
Subsequently, we performed an epistatic pairwise interaction among studied variants to evaluate if the investigated variants interact. Even if the interaction between AXIN1 and AXIN2 genes and between their codified proteins was described previously [13], our results did not find epistasis interactions between the investigated variants. No data have been available from other studies until now for variant interactions for AXIN1 rs1805105, rs12921862, and rs370681 and AXIN2 rs2240308.
Regarding the HWE, we observed that AXIN1 rs1805105 was not HWE in the CHDs group. However, the allele frequencies found in our study are similar to those reported in the literature for different populations [15], and according to Ensembl genome browser, all of our investigated variants had similar allele frequencies to those reported by the genome browser.
Briefly, our results showed that the presence of the variant genotypes of AXIN2 rs2240308 and the presence of the C-C-C and C-C-T haplotypes of AXIN1 (rs1805105, rs12921862, rs370681) are associated with a high risk of CHD development.
Concerning prenatal risk factors for congenital heart disease, an increase in maternal age at conception was observed (the mothers of CHD patients 26.9 ± 5.6 years and 25.1 ± 4.9 years for control group). Our observation is in line with that of Fung A et al. [34].
Maternal urinary tract infections were recorded in five cases, and no association was observed with CHD in the study group similarly to the results previously reported by Fung A et al. [34].
Although a family history of CHD, presence of associated congenital anomalies, maternal smoking, and alcohol consumption during pregnancy are reported to be associated with increased odds of CHD according to Fung A et al. [34], our study could not highlight these associations as long as we had cases with no history of CHD and cases with a clinical phenotype suggestive for a syndrome. Moreover, we did not include patients with other congenital malformations.
The strength of our study consists in the simultaneous analysis of four variants of AXIN genes. According to our knowledge, there are no other studies that performed a haplotype analysis and epistatic pairwise interaction between the investigated variants. Also, our study is the first one to evaluate the AXIN2 rs2240308 on CHDs patients. Moreover, a significant association between the variant genotypes of AXIN2 rs2240308 and CHD risk was noticed. On the other hand, our study has some limitations. One of these is the lack of measurement of the Axin protein levels in our subjects. The small number of subjects represents another limitation of our study, therefore, further large-scale studies are required to validate our results. Another limitation was the lack of whole exome sequencing analysis and data regarding socio-economic, educational status and body mass index of the parents.

Conclusions
The present study is an exploratory research that reveals significant positive associations between susceptibility to congenital heart defects and AXIN1 rs1805105, rs12921862 and rs370681 (C-C-C and C-C-T) haplotypes and AXIN2 rs2240308 variant.