Quantitative trait and transcriptome analysis of genetic complexity underpinning cardiac interatrial septation in mice using an advanced intercross line

Abstract

Unlike single-gene mutations leading to Mendelian conditions, common human diseases are likely to be emergent phenomena arising from multilayer, multiscale, and highly interconnected interactions. Atrial and ventricular septal defects are the most common forms of cardiac congenital anomalies in humans. Atrial septal defects (ASD) show an open communication between the left and right atria postnatally, potentially resulting in serious hemodynamic consequences if untreated. A milder form of atrial septal defect, patent foramen ovale (PFO), exists in about one-quarter of the human population, strongly associated with ischaemic stroke and migraine. The anatomic liabilities and genetic and molecular basis of atrial septal defects remain unclear. Here, we advance our previous analysis of atrial septal variation through quantitative trait locus (QTL) mapping of an advanced intercross line (AIL) established between the inbred QSi5 and 129T2/SvEms mouse strains, that show extremes of septal phenotypes. Analysis resolved 37 unique septal QTL with high overlap between QTL for distinct septal traits and PFO as a binary trait. Whole genome sequencing of parental strains and filtering identified predicted functional variants, including in known human congenital heart disease genes. Transcriptome analysis of developing septa revealed downregulation of networks involving ribosome, nucleosome, mitochondrial, and extracellular matrix biosynthesis in the 129T2/SvEms strain, potentially reflecting an essential role for growth and cellular maturation in septal development. Analysis of variant architecture across different gene features, including enhancers and promoters, provided evidence for the involvement of non-coding as well as protein-coding variants. Our study provides the first high-resolution picture of genetic complexity and network liability underlying common congenital heart disease, with relevance to human ASD and PFO.

Keywords: cardiac congenital defects; developmental biology; genetics; genomics; mouse; quantitative trait loci.

Figure 1.. Atrial septal morphology in normal adult and congenital heart disease.

(A) Schematic of human fetal heart highlighting inter-atrial septal structures. In the fetal heart, the septum primum has not yet fused with the septum secundum, as in adult stages; thus, the foramen ovale is patent. The foramen ovale and ostium secundum are offset, creating a flap valve. Arrow shows the direction of the blood shunt, allowing blood to bypass the pulmonary circuit. (B) Schematic of atrioventricular septal complex in the adult heart, and dysmorphology associated with the more common secundum type of atrial septal defect (ASD), and patent foramen ovale (PFO). Red arrow in the normal heart indicates the viewpoint for panel C. Double-headed black arrows indicate blood flow in pathological settings. (C) Light micrograph of the interatrial septum of an adult mouse heart as seen after removal of the left atrial appendage. The annulus of the mitral valve is towards the lower left of this panel. Note that the foramen ovale is covered by the membranous atrial septum primum, which is foremost in the left viewpoint (see red arrow in panel B). Septal landmarks are shown along with quantitative atrial septal traits measured in F2 and F14 studies. The crescent corresponds to the upper edge (see A) of the atrial septum primum. Scale bar, 100 µm. AV: atrioventricular; CRW: crescent width; DMP: dorsal mesenchymal protrusion; FO: foramen ovale; FOW: foramen ovale width; FVL: flap valve length; IVS: interventricular septum; LA: left atrium; LV: left ventricle; MC: mesenchymal cap; RA: right atrium; RV: right ventricle; SP: septum primum; SS: septum secundum.

Figure 1—figure supplement 1.. Histograms of quantitative traits in F14 mice with and without patent foramen ovale (PFO).

Figure 2.. Linkage analysis for heart weight (HW).

Suggestive and significant quantitative trait loci (QTL) identified on MMU7 (A) and MMU11 (B) by retrospective analysis of HW in the F2 study, the latter being fine-mapped by the advance intercross line (AIL) study. Y-axes represent logarithm of odds (LOD) scores for HW adjusted for age, sex, and body weight (BW).

Figure 2—figure supplement 1.. Comparison of advanced intercross line (AIL) results for heart weight (HW) (red line) and quantitative traits of atrial septum.

Figure 3.. Comparison of linkage results from F2 (red line) and advanced intercross line (AIL) (blue line) populations.

Y-axis represent logarithm of odds (LOD) scores and x-axis represent genetic map positions. A 1-LOD drop-off for each quantitative trait locus (QTL) is shown on the x-axis representing the confidence interval of the QTL.

Figure 3—figure supplement 1.. Differentiation of overlapping QTL peaks identified by the advanced intercross line (AIL) data.

To determine whether the peaks reflected the same underlying QTL or identified the presence of independent QTL, for each chromosome, the linkage analysis was performed using two different models in each, one peak marker was included as a fixed term. Vertical blue lines correspond to the positions of the AIL markers. The bold arrows indicate the positions of the fixed markers.

Figure 4.. Comparison of advanced intercross line (AIL) results for patent foramen ovale (PFO), flap valve length (FVL), foramen ovale (FOW), and crescent width (CRW).

Figure 5.. Segmentation of the quantitative trait loci (QTL) according to high variant density regions.

The QTL regions for each chromosome are illustrated as black blocks. The red histograms and y-axeis represent the variant density defined as the number of 129T vs QSi5 variants per 600 Kb genomics intervals. High variant density regions (density ≥ 1000) are represented as pink blocks. The x-axis for each chromosome represents the genomic position (Mb).

Figure 5—figure supplement 1.. Combined genome-wide variant zygosity of advanced intercross line (AIL) parental strains.

The majority of single nucleotide polymorphisms (SNPs) are either unique to 129T (1/1:0/0; 30%; dark blue) or QSi5 (0/0:1/1; 29.5%; yellow) or common to both lines (1/1:1/1; 27.5%; light blue) with respect to the C56Bl/6 reference. ‘Other’ includes, for example./.=one call in one line. The low number of heterozygous calls may be annotation artifacts due to the presence of repetitive sequences.

Figure 5—figure supplement 2.. Assessment of 129T2/SvEms vs QSi5 variant density across the genome using different bin sizes.

A bin of 600 Kb provided the clearest distinction between high and low variant density peaks defining a cutoff of 1000 variants/600 Kb.

Figure 6.. Identification of predicted-pathogenic variants.

(A, B) Number of variants per gene for predicted pathogenic genes under quantitative trait loci (QTL) according to classification (high-impact, non-frameshift deletion, or predicted-deleterious missense variant) for (A) 129T2 and (B) QSi5. (C, D) Overlap of genes containing predicted-pathogenic variants with (1) genes in a high confidence/emerging congenital heart disease (CHD) list, (2) genes with a heart development cardiovascular function from mouse genome informatics (MGI), and (3) all genes under QTL intervals that contain protein-modifying variants. Overlaps are shown for predicted pathogenic variants in (C) 129T2 and (D) QSi5.

Figure 7.. Functional analysis of SMAD6 variant R282P.

(A) Schematic of BMP signaling pathway highlighting roles of SMAD6. Prepared with Biorender. (B) Schematic of hSMAD6 showing the three main functional domains; the conserved DNA-binding MH1 domain (172-270aa), regulatory linker region (271-328aa), and the C-terminal transactivation/cofactor interaction MH2 domain (329-496aa). The position of G471 within the MH2 domain is indicated – mutation of this residue (G471S). Insert shows the amino acid sequence of the linker region containing PPXY (SMURF1 binding) and PLDLS (CtBP binding) domains, with the position of hSMAD R281P and Y279A variants indicated. Known serine phosphorylation sites are marked * above the sequence. (C–E). Functional analysis of hSMAD6 variants. Activity of hSMAD6-G471S (C), hSMAD6-Y279A (D), and hSMAD6-R281P (E), were assessed relative to wild-type (WT) hSMAD6 on a BMP-signaling responsive (BRE)-luc promoter in untreated and BMP4-treated cells. n=6–9; differences between WT and variant were assessed by unpaired t-test; data presented as mean ± SD. Asterisks denote significance, ****p<0.0001, ns = 0.3332. (F,G) Ubiquitination of WT hSMAD6 and variants assessed in MG132 or DMSO (control)-treated HEK293T cells. WT FLAG-ShMAD6 and variants were immunoprecipitated with anti-FLAG antibody with lysates subject to western blotting with anti-ubiquitin antibody (F). Total ubiquitin, FLAG-hSMAD6, and protein indicated. (G) Quantification of FLAG-hSMAD6 ubiquitination. n=3; difference between WT and variants assessed by one-way ANOVA with Tukey’s multiple comparisons; data presented as mean ± SD. Asterisks denote significance, ***p=0.0001.

Figure 8.. Transcriptome analysis of atrial septum.

(A) Principal component analysis (PCA) plot of gene expression profiles from atrial septum RNA-seq libraries including two replicates per time point (E12.5, E14.5, and E16.5) per mouse strain (QSi5 and 129T2/SvEms) for protein-coding genes. (B) Venn diagram with the number of differentially expressed genes between QSi5 and 129T2/SvEms at different time points. (C) MA plots of differentially expressed genes for each time point. (D) STRING network of E14.5 differentially expressed genes and genes with predicted-pathogenic variants from QSi5 or 129T2/SvEms. Genes are colored according to whether they are upregulated or downregulated in 129T2/SvEms or contain a variant that is predicted-pathogenic and differentially expressed.

Figure 8—figure supplement 1.. STRING network of E12.5 differentially expressed genes and genes with predicted-pathogenic variants from QSi5 or 129T2/SvEms.

Genes are coloured according to whether they are upregulated or downregulated in 129T2/SvEms, or contain a variant that is predicted-pathogenic and differentially expressed.

Figure 8—figure supplement 2.. STRING network of E16.5 differentially expressed genes and genes with predicted-pathogenic variants from QSi5 or 129T2/SvEms.

Genes are coloured according to whether they are upregulated or downregulated in 129T2/SvEms, or contain a variant that is predicted-pathogenic and differentially expressed.

Figure 9.. Network analysis of protein kinase A (PKA) and AKT1 phosphorylation targets.

Integrated network analysis of PKA and AKT1 STRING connections and phosphorylation substrate targets from PhosphoSitePlus and PhosphoPICK predictions.

Figure 10.. Enrichment of 129T2/SvEms vs QSi5 variants in genomic features of differentially expressed genes (DEGs) under quantitative trait loci (QTL).

The number of variants directly overlapping a feature is shown in blue, and gray bars represent the expected values based on the mean overlap from 1000 randomly generated interval sets. Error bars show the 95% confidence intervals of the mean. The significance of the enrichment is expressed as p-values, calculated by dividing the number of random samples showing equal or greater overlap than the observed by the total number of permutations (*p< 0.05).

Figure 10—figure supplement 1.. Gene counts by number of 129T2/SvEms vs QSi5 variants.

For each genomic feature, a density plot and histogram representing the log-transformed number of variants per gene (x-axis) and frequency of the genes (y-axis). Each plot illustrates the gene counts for differentially expressed genes (DEGs) (red) and non-DEGs (blue) under QTL.

Figure 10—figure supplement 2.. Impact of increasing enhancer window size on over-representation of variants in differentially expressed genes (DEGs) and non-DEGs under quantitative trait locus (QTL), and on number of enhancers detected.

(A) Total number of enhancers detected with increasing window size. (B) Total number of enhancers corresponding to DEGs and non-DEGs with increasing window size. (C) Effect of increasing enhancer window size on over-representation of variants in DEGs and non-DEGs under QTL.