Congenital heart diseases (CHDs) are the most common type of congenital anomaly and the leading noninfectious cause of infant death (Liu et al., Reference Liu, Chen, Zuhlke, Black, Choy, Li and Keavney2019). The mean prevalence of CHD 1970–2017 globally was 8.224 (7.817, 8.641) per thousand live births (Liu et al., Reference Liu, Chen, Zuhlke, Black, Choy, Li and Keavney2019). CHDs are much more common in twin pregnancies, with a reported prevalence of 20 in 1000 live births (Balasubramanian et al., Reference Balasubramanian, Vuppalapati, Avanthika, Jhaveri, Peddi, Ahmed, Reddy and Kaur2021).
Double outlet right ventricle (DORV), a condition when both great arteries arise predominantly from the right ventricle, is a rare form of CHDs. In China, the birth prevalence of total CHDs has increased continuously over time, from 0.201% in 1980–1984 to 4.905% in 2015–2019; the birth prevalence of DORV is 0.059% (L. Zhao et al., Reference Zhao, Chang, Liu, Wang, Zhang, Lu, Zhang, Zhang and Wang2020). Consequently, CHDs cause considerable suffering to patients and their families, and have become a sizable public health concern.
The etiology of the majority of CHDs remains largely unknown (van der Bom et al., Reference van der Bom, Zomer, Zwinderman, Meijboom, Bouma and Mulder2011). It is widely believed that most CHDs arise from the interplay between genetic and environmental factors (Krauss & Hong, Reference Krauss and Hong2016; Tiret, Reference Tiret2002). Epigenetics has been proposed to be one of the main mediators of this interaction, and extensive studies have identified aberrant epigenetic patterns in several polygenic diseases, including CHDs (Cao et al., Reference Cao, Wu, Huang, Wang, Su and Ye2021; Feinberg, Reference Feinberg2018; Gilsbach et al., Reference Gilsbach, Preissl, Gruning, Schnick, Burger, Benes, Wurch, Bonisch, Gunther, Backofen, Fleischmann, Schubeler and Hein2014; Grunert et al., Reference Grunert, Appelt, Grossfeld and Sperling2020; Lyu et al., Reference Lyu, Zhang, Ling, Liu, Zong, Guan, Huang, Sun, Zhang, Li, Nie and Tao2018). Epigenetics involves the control of transcription and translation without changing the nucleotide sequence, and includes DNA methylation and histone modifications, which control gene accessibility, and noncoding RNAs, which primarily control mRNA translation.
DNA methylation, defined as cytosine residues methylation in CpG dinucleotides of DNA molecules, is the most widely explored epigenetic modification. In CHD patients, not only abnormal genomewide methylation has been detected (Grunert et al., Reference Grunert, Dorn, Cui, Dunkel, Schulz, Schoenhals, Sun, Berger, Chen and Sperling2016; Zhou et al., Reference Zhou, Xiong, Dong, Wang, Qian, Ma and Li2021), but also aberrations of methylation at specific genes have been reported, including imprinted genes (Chang et al., Reference Chang, Wang, Xin, Wang, Luo, Wang, Zhang and Li2021; X. Zhao et al., Reference Zhao, Chang, Liu, Wang, Zhang, Lu, Zhang, Zhang and Wang2020), histone-modifying genes (Yi et al., Reference Yi, Jiang, Li and Jiang2017), and heart development-related genes, such as CITED2 (M. Xu et al., Reference Xu, Wu, Li, Yang, Hu, Zheng and Tian2014), NKX2-5, GATA4, HAND1 (Sheng et al., Reference Sheng, Wang, Ma, Qian, Zhang, Wu, Zheng, Chen, Huang and Ma2012), TBX20 (Gong et al., Reference Gong, Sheng, Ma, Huang and Liu2019), RXRA4 (Zhang et al., Reference Zhang, Ma, Wang, Ma and Huang2014), ZFPM2 (Sheng et al., Reference Sheng, Chen, Wang, Ma, Ma and Huang2016). Long noncoding RNAs (lncRNAs) are transcripts >200 nucleotides long that have no known protein-coding function. Animal experimental studies have suggested that several lncRNAs play important roles in cardiac development (Anderson et al., Reference Anderson, Anderson, McAnally, Shelton, Bassel-Duby and Olson2016; Kim et al., Reference Kim, Lee, Son, Nam, Moon, Pyun, Park, Kang, Lee and Cho2021; Ponnusamy et al., Reference Ponnusamy, Liu, Zhang, Li, Zhai, Liu, Zhou, Liu, Yan, Dong, Wang, Qian, Shan, Xu, Wang, Zhang, Li, Zhang and Wang2019). However, only one study have shown the different expression of up to 1500 lncRNAs between normal hearts and hearts from fetuses with ventricular septal defect (Song et al., Reference Song, Shen, Zhu, Liu, Liu, Shen, Zhu, Kong, Yu and Qian2013).
Monozygotic (MZ) twins share almost all of their genetic information stored in DNA, but are often discordant for common complex diseases, such as diabetes, autism, schizophrenia, different types of cancer, as well as CHDs (Grunert et al., Reference Grunert, Appelt, Grossfeld and Sperling2020; Lyu et al., Reference Lyu, Zhang, Ling, Liu, Zong, Guan, Huang, Sun, Zhang, Li, Nie and Tao2018). MZ twins are advantageous material for studying epigenetics by controlling for many potential confounders (Castillo-Fernandez et al., Reference Castillo-Fernandez, Spector and Bell2014). One study conducting genomewide sequencing and methylation analysis in MZ twins discordant for DORV identified few genomic differences but 1566 differentially methylated regions (DMRs) showed a high correlation between hypermethylated promoters at ZIC3 and NR2F2 and downregulated gene expression levels of these two genes in patients with DORV compared to normal controls (Lyu et al., Reference Lyu, Zhang, Ling, Liu, Zong, Guan, Huang, Sun, Zhang, Li, Nie and Tao2018). Another study performing genomewide high-throughput sequencing in MZ twins discordant for tetralogy of Fallot (TOF) observed DNA methylation changes in regulatory regions of known cardiac-relevant genes (Grunert et al., Reference Grunert, Appelt, Grossfeld and Sperling2020).
In the present study, we obtained cardiac tissue samples from one MZ twin pair discordant for DORV, and comprehensively investigated genomewide genetic, DNA methylation, mRNA, lncRNA and protein-level differences, providing further insights into the etiology of CHDs
Materials and Methods
Study Subjects
The discordant monochorionic diamniotic twins were urgently delivered at 35+6 weeks gestation due to severe maternal preeclampsia. However, the twins died within 5 minutes after birth because of fetal status after resuscitation. Twin A had multiple congenital heart defects, including a DORV, with a ventricular septal defect, interrupted aortic arch, hypoplasia of the ascending aorta, thickening of the tricuspid valve leaflets, and right atrial dilation, while twin B had normal cardiovascular structures, but a diagnosis of selective intrauterine growth restriction (sIUGR; Figure S1). There was no evidence that twin–twin transfusion syndrome was applicable, based on the previously published Quintero criteria (Benoit & Baschat, Reference Benoit and Baschat2014). The characteristics of the twins are detailed in Table S1. Monozygosity was confirmed by zygosity identification analysis. Cardiac and muscle tissue samples were collected from the same region of the conusarteriosus and the deltoid muscle of the right arm of the twins. Microsatellite analysis was used to confirm the zygosity of the twins (methodological details are given in Supplementary Additional file 3). Sporadic abortion cases were collected from Shenzhen Maternity and Child Healthcare Hospital and Xijing Hospital. All cases were all diagnosed as DORV by a professional pathologist and a cardiac surgeon.
Whole Genome Sequencing (WGS)
Genomic DNA from conusarteriosus tissue of the MZ twins was extracted using DNeasy Blood and Tissue Kit (Qiagen). A standard CG two-adapter library was constructed. Using short-sequence and paired-end sequencing strategies, genomes were sequenced with Illumina HiseqTM 2000 (Illumina Inc, San Diego, CA, USA). Data processing included generating base sequences, alignment, local assembly of candidate sites, identification of variants, and QC of the results. CGA tools calldiff (http://www.completegenomics.com/analysis-tools/cgatools/) were used to compare the differences between the twins. The statistical results of WGS data output are shown in Table S2.
Whole Genome Bisulfite Sequencing (WGBS)
Genomic DNA from conusarteriosus and skeletal muscle tissue of the MZ twins were processed as follows: fragmentation, end repair, adapters ligation, bisulfite treatment, PCR amplification, library test. The processed libraries were then subjected to paired-end 90bp and single-end 50bp sequencing respectively by Illumina HiseqTM 2000 sequencing system according to the manufacturer’s instructions. The details regarding the sequencing and bioinformatics analysis are available in Supplementary Additional file 5.
The overall data output of WGBS, and the sequencing depth and coverage of CG sites in genomic elements and chromosomes are shown in Tables S3 and S4 respectively.
DMRs were screened with DMR different level ≥ 0.2 and Fisher exact test p value ≤ .005 as filter conditions. Due to the relatively extensive research on the methylation function of the promoter regions, the DMRs that were significant in conus-case and conus-control pairs but not significant in muscle-case and muscle-control pairs were intersected with the gene promoter region, and the gene was outputted when the intersection was 100%.
RNA Sequencing
Total RNA from conusarteriosus and skeletal muscle tissue of MZ twin were processed as follows: purification, integrity assessment, removing ribosomal RNA (rRNA), random fragmentation, and reverse transcription, then synthesized into a double-stranded cDNA, terminal repair (blunt-ending), 3’ terminal addition of A-residues, adaptor ligation, PCR amplification, purification and quality inspection of the library, cluster generation. The libraries were then subjected to paired-end sequencing by Illumina HiseqTM 2000 sequencing system. The quality control standard of sequencing was as follows: the amount of data was about 15G per sample, and the ratio of base quality greater than 20 (Q20) in each direction was not less than 90%. The details regarding the sequencing and bioinformatics analysis of mRNA and lncRNA are available in Supplementary Additional file 8. Differences were compared using edgeR (Robinson et al., Reference Robinson, McCarthy and Smyth2010).
Proteomic Expression Analysis
Four protein samples from the conusarteriosus and skeletal muscle of both twins were processed as follows: enzymolysis, desaltation, iTRAQ labeling, separation of components, and LC-MS/MS analysis. The details regarding the protein expression bioinformatics were available in Supplementary Additional file 9.
Three technical replica were carried out. A summary of the three replicate identifications of proteins are shown in Table S5. The number of peptides and proteins identified three times were 5820 and 1221, accounting for 48.66% and 55.58% of the total peptides (N 11961) and proteins (N 2197) respectively.
Gene Ontology Enrichment and Pathway Analysis
Gene ontology (GO) enrichment (http://www.geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses (http://www.genome.jp/kegg/pathway.html) were applied to determine the functional roles and related pathways of differentially methylated regions, expressed genes or transcripts, or proteins. Whole genome sequencing results were used as enrichment background. Gene networks and canonical pathways representing key genes were identified using the Database for Annotation, Visualization and Integrated Discovery (DAVID; Sherman et al., Reference Sherman, Hao, Qiu, Jiao, Baseler, Lane, Imamichi and Chang2022). Fisher’s exact test and a chi-square test were used to select the significant pathway, and the false discovery rate (FDR) was applied to account for multiple comparisons. The threshold of significance was defined by the p value (p ≤ .005) and FDR (q ≤ 0.20).
Validation of Experiments
Sanger sequencing, liquid chromatography-tandem mass spectrometry (LC-MS/MS; Zhang et al., Reference Zhang, Zhang, Zhou, Ye, Liu, Kang and Cai2011), bisulfite sequencing PCR (BSP), real-time quantitative PCR (qPCR), and multiple reaction monitoring (MRM) were used to validate genomic variants, global DNA methylation, DNA methylation at specific loci, mRNA expression, target proteins, respectively.
Five abortion DORV cases were selected to validate gene expression. Using conus as lesion case and left ventricular as control, qPCR was performed to assess mRNA expression at specific genes. On the other hand, nine sporadic DORV cases and five normal controls were used to validate variation in the methylation levels of specific genes. DNA was extracted from conusarteriosus tissues to measure methylation levels by BSP. Differences in methylation levels between cases and controls were compared, to determine whether these followed the trend for differences between the twins.
Human umbilical vein endothelial cells (HUVECs) were cultured and treated with 4 µM 5-Aza-dC (Sigma–Aldrich, St Louis, MO, USA) for 48 h. Cells were harvested, then BSP and qPCR were performed to assess DNA methylation at specific loci and to validate mRNA expression.
Results
Genomic Variations in MZ Twin Discordant for DORV
The number of differential SNVs between MZ twins was 18,824. Of these, the number of point mutation in coding sequence (CDS), introns, upstream of transcription start sites (TSS-upstream) and untranslated regions (UTRs) were 128, 6302, 685 and 96 respectively. Of the 128 differential SNVs in the CDS region, there were 69 nonsynonymous mutations, of which 29 mutations had a somatic score above -15 (Table S6). The mRNA expression analysis of these 29 SNVs showed that only the NRAP gene had differential mRNA expression. However, the result of Sanger validation showed that the twins had same genotype in NRAP gene; both were homozygous for TT.
Based on the WGS data, genomewide copy number variations (CNVs) were also screened. The number of specific CNVs for case and control was 9 and 13 respectively (Table S7). Among the differential CNVs between twins, only a few overlapped with gene regions, but have not been reported to be associated with heart disease.
DNA Methylation Differences Between Discordant Twin
WGBS revealed that the global methylation level of the affected twin was lower than the unaffected twin in both the conusarteriosus and muscle tissues (Figure 1A). Regarding the CG methylation sites, the global methylation level was 69.10% and 70.38% in the conusarteriosus of the affected and unaffected twin respectively, while that in the skeletal muscle was 69.36% and 70.05% in the affected and unaffected twin respectively. The resulting differences in the global CpG methylation levels between the twins that were relative larger in arterial conus and smaller in muscle tissue were confirmed by LC-MS/MS validation experiments (Figure 1B). The average methylation levels in genomic elements and chromosome are shown in Table S8 and Table S9 respectively. In terms of genomic regions, discordance in CpG methylation was significantly different between the two tissue types in promoter regions, CDSs, introns, and 3′-UTRs, but not in the 5′-UTRs or in the intergenic regions (Figure 1C).
Comparison of DNA methylation levels between the MZ twin and the distribution of DMR. (A) Overview the average genome-scale DNA methylation levels by WGBS. (B) The average genome-scale DNA methylation levels by LC-MS/MS. (C) The differences of methylation levels in genomic elements by WGBS. (D) The distribution of DMR in promoter elements. (E) The hypomethylation of COL1A1 and THBS3 by BSP.
According to the screening criteria, 36,228 DMRs that only occurred in the conusarteriosus tissue were obtained compared to skeletal muscle. A total of 1094 DMRs were listed in promoter region. The distribution of DMRs in promoter elements is shown in Figure 1D. Of these, 3% were located in promoter regions, involving 1039 genes. The GO enrichment in promoter DMR-associated genes in the twin case are shown in Table S10. Functional pathway analysis revealed 44 DMR-associated genes, involved in eight pathways as defined by the KEGG database, as shown in Table 1. Among these, the pathway involving ECM-receptor interaction was the most significant (EASE score = 9.05E-03).The hypomethylation of specific genes, such as COL1A1 and THBS3, were confirmed by BSP validation data (Figure 1E).
Pathway analysis of promoter DMR associated genes
mRNA and lncRNA Expression Levels Between Discordant Twin
RNA-Sequencing revealed that 419 genes that were significantly differentially expressed with an expression fold-change >2 in the conusarteriosus and fold-change <2 in the skeletal muscle between the twin. The expression of these are listed in Table S11. The highly expressed genes by, such as COL1A1 and THBS3, were confirmed by qPCR validation data (Figure 2). The GO enrichment of differentially expressed genes are shown in Table S12. Functional pathway analysis revealed 10 pathways involving 46 genes, as shown in Table 2. Among these pathways, changes in focal adhesion pathways were the most significant (EASE score = 3.29E-08). The expression level of COL1A1and THBS3 by qPCR data were confirmed to be consistent with RNA-Seq results (Figure 2).
The expression level of COL1A1and THBS3 by qPCR data were confirmed to be consistent with RNA-Seq results.
Pathway analysis of differentially expressed genes
Thirteen novel and 33,370 known lncRNA transcripts were obtained. The novel lncRNA transcripts and their corresponding genes are displayed in Table S13. Among these, a total of 1818 lncRNAs were differentially expressed in the conusarteriosus between the twin, including 755 up-regulated and 1063 down-regulated lncRNAs in the case twin. Among the potential target genes of differentially expressed lncRNAs, only 30 mRNAs were differentially expressed; the lncRNA-mRNA pairs are shown in Table S14.
Differential Proteins Between Discordant Twins
Out of the 1676 proteins that were repeatedly detected by LC-MS/MS, 62 proteins that were significantly different in the conusarteriosus tissue (fold-change ≥1.5 and p value ≤ .05) but not in muscle tissue were obtained. These protein-encoding genes function are shown in Table S15. Functional pathway analysis revealed that these genes clustered in eight pathways, as shown in Table 3, among which the extracellular matrix (ECM) –receptor interaction pathway was the most significant (EASE score = 8.89E-10). The MRM results showed that the specific peptide segment of COL1A1 protein was 2.48 times higher in the case twin than in the control twin; this confirmed COL1A1 higher expression by LC-MS/MS.
Pathway analysis of differential protein-encoding genes
The Common Genes and Pathways of Multiple Omics Integrative Analysis
Among the 1039 promoter DMR-associated genes, 42 genes and 9 genes were differently expressed in mRNA and protein levels respectively. Five genes, including BGN, COL1A1, COL3A1, FBLN5 and FLAN, exhibited differences at all three levels. According to the analysis results of DNA promoter methylation, mRNA and protein expression, the pathways that were enriched in all three levels were ECM-receptor interaction, focal adhesion and TGF-β signaling pathway. The three enriched pathways and their corresponding genes that differed at different levels are shown in Figure 3. COL1A1 and COL3A1 participated in ECM-receptor interaction and focal adhesion pathway, and FLNA participated in focal adhesion pathway. Although THBS3 participated in all three pathways, it only showed differences at the methylation and mRNA level.
The three common pathway and genes at three levels.
The Methylation and mRNA Expression Validation
The validation experiment in nine sporadic cases of DORV and five aborted normal controls showed that the average methylation levels at the promoter region of COL1A1 and THBS3 in cases were both lower than controls, following the trend seen in the twins (Figure 4A). The results in the cell model showed that after treatment with 5-Aza-dC, the methylation levels of COL1A1 and THBS3 genes decreased in the HUVECs. Concurrently, the expression of these two genes was increased (Figure 4B). Compared with the control group, 5-Aza-dC treatment led to 25.2-fold and 1.6-fold up-regulated expression of COL1A1 and THBS3 respectively. These results demonstrated that the hypomethylation of COL1A1 and THBS3 may cause up-regulation of their expression.
The methylation and expression validation of COL1A1 and THBS3. A: The methylation level in sporadic cases and controls. B: The methylation and mRNA expression in cell model.
Discussion
In this study, we investigated specific tissues from discordant MZ twins, using multiomics strategies based on high-throughput techniques, which enabled us to perform the first systematic analysis of genome, methylation, transcriptome, proteome-wide molecular changes in CHDs.
We performed genomic sequence variation detection to identify post-twinning SNVs, short indel and CNV differences between the twins that could make one of them susceptible to CHD. However, we did not detect any discordant coding or regulatory SNV or indel variants or CNVs that would have been likely to cause CHD. The results were consistent with the outcomes of previous studies that were conducted in DORV (Lyu et al., Reference Lyu, Zhang, Ling, Liu, Zong, Guan, Huang, Sun, Zhang, Li, Nie and Tao2018), TOF (Grunert et al., Reference Grunert, Appelt, Grossfeld and Sperling2020) or CHDs (Y. Xu et al., Reference Xu, Li, Pu, Cao, Long, Chen, Sun and Xu2017) discordant monozygotic twins. It has been reported that during the twinning process and/or embryonic development, the underlying genetic difference could be involved. For example, one study identified three CNV differences involving copy number gains in the affected member of a MZ twin discordant for CHDs, including 3.7 kb on 12p13.31, 7.5 kb on Xp11.23, and 37.4 kb on Xq28 (Breckpot et al., Reference Breckpot, Thienpont, Gewillig, Allegaert, Vermeesch and Devriendt2012). In addition, studies have reported de novo somatic mutations in CHD-discordant monozygotic twins,(Li et al., Reference Li, Montpetit, Rousseau, Wu, Greenwood, Spector, Pollak, Polychronakos and Richards2014; Zhang et al., Reference Zhang, Thiele, Bartmann, Hilger, Berg, Herberg, Klingmuller, Nurnberg, Ludwig and Reutter2016).
DNA methylation is an epigenetic modification that plays a critical role in heart development (Akerberg & Pu, Reference Akerberg and Pu2020; Gilsbach et al., Reference Gilsbach, Preissl, Gruning, Schnick, Burger, Benes, Wurch, Bonisch, Gunther, Backofen, Fleischmann, Schubeler and Hein2014; Jarrell et al., Reference Jarrell, Lennon and Jacot2019). In the present study, global hypomethylation status was observed in the heart tissue sample of the twin case with DORV. This is consistent with the results of two previous case-control studies: one assessed long interspersed nucleotide elements (LINE-1) methylation, which was representative of genomewide methylation status in right ventricular tissue samples from pediatric patients (Sheng et al., Reference Sheng, Wang, Ma, Qian, Zhang, Wu, Zheng, Chen, Huang and Ma2012), while the other conducted genomewide DNA methylation analysis in myocardium tissue samples from a fetal cases with cardiac defects (Zhou et al., Reference Zhou, Xiong, Dong, Wang, Qian, Ma and Li2021). However, one previous genomewide DNA methylation study found more hypermethylated than hypomethylated DMRs in myocardial samples from pediatric patients with CHDs compared with controls (Grunert et al., Reference Grunert, Dorn, Cui, Dunkel, Schulz, Schoenhals, Sun, Berger, Chen and Sperling2016); another study using heart tissue DNA from fetuses presenting a variety of cardiac defect phenotypes, including DORV, found no significant differences in whole genome methylation pattern between fetuses with isolated cardiac defects and fetuses with normal development (Serra-Juhe et al., Reference Serra-Juhe, Cusco, Homs, Flores, Toran and Perez-Jurado2015).
DNA methylation can change the normal gene expression and regulate crucial aspects of its function. One study revealed that hypermethylation of the NOTCH4 promoter region was associated with decreased NOTCH4 expression in patients with TOF (Zhu et al., Reference Zhu, Ye, Xu, Gu, Ma, Chen, Li, Sheng and Huang2020). Two studies demonstrated that aberrant methylation levels at the promoter CpG island shore of the ZFPM2 and HAND1 genes might be responsible for gene transcription regulation in patients with TOF (Sheng et al., Reference Sheng, Chen, Wang, Ma, Ma and Huang2016; Sheng et al., Reference Sheng, Qian, Wang, Ma, Zhang, Diao, An, Chen, Ma and Huang2013). Another study suggested that demethylation in the TBX20 promoter region may be responsible for overexpression in the cardiac tissue of patients with TOF (Gong et al., Reference Gong, Sheng, Ma, Huang and Liu2019). In this study, hypomethylation and high mRNA expression level of COL1A1 and THBS3 genes were detected in the affected twin and sporadic cases. The cell model experiment demonstrated that decreased methylation levels of COL1A1 and THBS3 genes could lead to increased expression levels. Although there were no reports on the association of THBS3 and COL1A1 genes with cardiac development or CHDs, our results suggested that THBS3, a common gene in all three pathways, and COL1A1, a common gene in the first two pathways, might have potential functions in cardiac development.
In the multiomics analysis, five genes, including BGN, COL1A1, COL3A1, FBLN5 and FLAN, all showed differences in methylation, mRNA expression and protein level. Among them, only COL3A1 has been reported to be associated with the risk of CHDs (Morton et al., Reference Morton, Shimamura, Newburger, Opotowsky, Quiat, Pereira, Jin, Gurvitz, Brueckner, Chung, Shen, Bernstein, Gelb, Giardini, Goldmuntz, Kim, Lifton, Porter, Srivastava and Seidman2021; Tan et al., Reference Tan, Wang, Liu, Peng, Yang, Yu, Xu, Li and Cao2022). In addition, three critical pathways, viz. ECM-receptor interaction, focal adhesion and TGF-β signaling, were dysfunctional in the affected twin and enriched in all three levels. Recent studies have revealed that these pathways may play a role in appropriate cardiac outflow tract morphogenesis and aortic arch remodeling, as they may be involved in differentiation, migration, and survival of cardiac neural crest cells (Restivo et al., Reference Restivo, Piacentini, Placidi, Saffirio and Marino2006) and the second heart field (Dyer & Kirby, Reference Dyer and Kirby2009). Consequently, dysfunction within these pathways may lead to short outflow tracts (Lewandowski et al., Reference Lewandowski, Janardhan and Trivedi2015; Lu et al., Reference Lu, Lin, Kao, Chung, Yeh, Chen and Chen2016), including DORV. According to another theory of left-right (LR) asymmetry, DORV may be caused by the failure of shifting of the conotruncal septum to the left, so that both the aorta and the pulmonary artery would arise from the right ventricle (Rahman et al., Reference Rahman, Zhang, Fan and Wan2020). During this process, TGF-β signaling is important for cardiac looping with appropriate L-R asymmetry (Tadjuidje et al., Reference Tadjuidje, Kofron, Mir, Wylie, Heasman and Cha2016). The ECM and focal adhesion pathways have been reported to be related to heart development and malformation.
This study had several strengths. First, our unique sample consisting of a MZ twin pair discordant for DORV allowed us to detect the epigenetic variation controlling for genotype, age, sex and other potential confounders, and provide credible results. Second, given that epigenetic regulation was tissue-specific, we have produced a tissue-specific multiomics map using tissue samples of the developing heart, which could enhance our capability to pinpoint the underlying mechanism of CHDs. Third, by undertaking a genomewide approach, we were able to uncover phenotype-relevant differentially methylated loci in genomic regions that are both novel and have been previously associated with CHDs. Fourth, compared with the previous studies of MZ twins discordant for CHDs, which only carried out genomewide genetic and DNA methylation analysis (Grunert et al., Reference Grunert, Appelt, Grossfeld and Sperling2020; Lyu et al., Reference Lyu, Zhang, Ling, Liu, Zong, Guan, Huang, Sun, Zhang, Li, Nie and Tao2018), this study systematically performed genomewide genetic and DNA methylation, and transcriptome and protein level analysis, providing a deeper understanding of the etiology of this complex disease.
Our study has several limitations. First, we only validated COL1A1 and THBS3. The other potentially related genes were not validated in the twins and sporadic cases. Second, the sample size of this study was relatively small, and the control fetus in the discordant twins was not completely normal but had iUGR. Although we did the validation experience on sporadic cases, the results on sequencing from only one pair of twins is hard to replicate under other circumstances. Ideally, we would hope to expand the sample size in future studies. Third, due to the relatively stringent screening criteria, genes or pathways that were different in methylation, mRNA expression and protein levels were considered candidates; this potentially reduced the number of candidates that need to be validated. Moreover, the differentially expressed lncRNA-mRNA pairs found in the twin requires further validation in additional cases.
In conclusion, this study is the first time a tissue-specific multiomics profile of DORV was delineated. Global hypomethylation might be associated with the occurrence of DORV. Multiomics molecular analysis suggested that several genes, such as BGN, COL1A1, COL3A1, FBLN5 and FLAN, and some pathways, particularly involving ECM–receptor interaction, focal adhesion and TGF–β signaling, might be involved in the occurrence of CHDs. These results will provide basic data resources for subsequent studies on the pathogenesis of CHDs. However, further studies are still required to investigate their mechanisms in the occurrence and development of CHDs.
Supplementary material
To view supplementary material for this article, please visit https://doi.org/10.1017/thg.2023.51.
Data availability statement
The data involved in this manuscript are listed in the supplements, which can be available for researchers. The raw data also can be obtained by contacting the corresponding author by emails.
Acknowledgments
We are indebted to the pediatrician, pathologists and experimental technicians who were involved in case diagnosis and multiomics experimental detection. We thank all participating families for their cooperation and for providing personal information and biosamples. We also thank the reviewers for their helpful comments.
Author contribution
Dezhi Mu and Jun Zhu contributed to the conception and design of this study; Zhen Liu and Nana Li were responsible for the respondents’ enrollment and specimen collection and drafted the manuscript; Jun Li and Shengli Li performed accurate clinical diagnosis; Xiaoyu Pan, Ping Li, Ying Deng and Ping Yu assisted with the experiments; Qintong Li and Fang Chen performed the statistical analysis; Hui Jiang, Wei Wang, Dezhi Mu, Ping Yu and Jun Zhu critically reviewed the manuscript and supervised the entire study. All authors read and approved the final manuscript.
Financial support
This study was supported by the National Natural Science Foundation of China (grant ID: 81602865, 82103858, 81273086), Major Scientific and Technological Research Projects, Ministry of Education of China (grant ID: 313037), the National Science & Technology basic work Project of the Ministry of Science and Technology of China (grant ID: 2014FY110700), and Sichuan Science and Technology Program (grant numbers 2022NSFSC0661, 2022NSFSC0656).
Competing interests
The authors declare no competing interests.
Ethics approval and consent to participate
Parental permissions and informed consent were obtained from the parents of all patients. This research was approved by the Ethics Committee of Sichuan University (No. 2010004) and followed the tenets of the Declaration of Helsinki.
Open access
