Genomic structural variations lead to dysregulation of important coding and non-coding RNA species in dilated cardiomyopathy

Study population for multi-omics analysis

For this multi-omics analysis, it was required that high-quality material and biopsies were present of each DCM patient (Fig 1A). Included patients with suspicion for primary DCM underwent extensive clinical phenotyping including coronary angiography with myocardial biopsy (Meder et al, 2017). Exclusion criteria were all secondary causes of DCM. A total of n = 50 consecutive patients fulfilling the requirements were included, and detailed clinical baseline characteristics are summarized in Table 1.

NGS was performed on an Illumina platform, and libraries were generated from biopsy and blood using the TruSeq technology. For WGS, the library inserts are small spread, averaged on 296 bp (224–327 bp), which is necessary to fulfill technical requirements to reliably call large deletions, duplications, or inversions with high confidence (Fig 1B), especially in genomic regions of low sequence complexity and regions with difficulties for sequencing and mapping. Deep sequencing gained a coverage ranging from 40× to 66× with an average of 58× ± 6.5. Figure 2A shows the uniformity of the generated sequencing data over all chromosomes (blue line).

High-resolution map of genomic structural variants in DCM

The high-coverage WGS of DCM patients was subsequently subjected to paired-end mapping and split-read detection algorithms for the identification of SVs. The algorithms are implemented in Delly SV calling, an approach also employed in the latest release of the 1000 Genomes Project (Rausch et al, 2012; Sudmant et al, 2015). Inclusion of deletions and duplications with respective read-depth and inversions with support for both breakpoint and complex events resulted in 2,955 high confidence deletions, 797 duplications, and 146 inversions and complex structural variants in the 50 patients. SVs from small to large events in a range of 212 bp to 12.2 Mb (Fig 2B) were detected (deletions: 249 bp to 12.2 Mb, duplications: 212 bp to 1.7 Mb, and inversions and complex structural variants: 216 bp to 44.9 kb). Deletions and duplications were relatively uniformly distributed over all chromosomes and showed an expected modest accumulation in low complex regions, the telomeres, and centromeres (Fig 2A). To further increase confidence on the quality of the SV calls, we utilized raw intensity data from Illumina 450K methylation arrays measurements from the same samples (decreased signal intensities for deletions and increased intensities for duplications) (Feber et al, 2014; Meder et al, 2017). From a total of 633 SVs, in which sufficient methylation sites were available, we could independently validate the majority of deletions and duplications (P-value ≤ 0.05), showing the high quality of the variant calls (Fig EV1A), which was also underlined by PCR validation of randomly selected SV events (Fig EV1B).

Click here to expand this figure.

To our knowledge, we present the most detailed map of genetic variation in a well-phenotyped cardiomyopathy cohort and find that the genome-wide SVs delete important loci directly related to gene expression, such as enhancers (n = 875), TFBSs (66,147), and lncRNAs (3,100) (Fig 2C). The fluctuating numbers of these affected regulatory loci in each patient are also visualized and can be mainly explained by the presence of some very large deletion events in some individuals. Most of the SVs reside within non-coding, intergenic regions. However, we find in total 4,305 exons (1,367 genes) to be entirely deleted by a hetero- or homozygous deletion, together with 630 deletions to be located inside an intron (Fig 2D). Duplications in DCM affect 132 enhancers, 9,674 TFBSs, 658 lncRNAs, and 389 protein-coding genes with 319 exons and 199 duplicated introns (Fig 2C and D).

Deletions of protein-coding exons result in average in a significant downregulation of the transcript (median expression value 0.8 for SV carriers) (Fig 2D, violin plots). Duplications do not show significant changes over all events, most likely due to the incomplete event not carrying the complete promoter or enhancer structure of the host gene. The effect of different SV-allele counts of a deletion on RNA-seq profiles is exemplarily shown in Fig 2E. In this case, the deletion leads to a reduction in exon expression of heterozygous SV carriers (orange line) and a complete absence of RNA expression in the deleted region of homozygous patients (red line, Fig 2E). 23.7% of the protein-coding genes hit by a SV were entirely covered by the genetic variant. In detail, 377 genes were completely deleted at least on one of two chromosomes (126 genes are on the X-chromosome affected by a large deletion that is known from population studies), and 40 genes were completely duplicated. Overall, our results show that each individual carries more than 720 SVs on average.

Single nucleotide and structural variants are linked to myocardial transcript dysregulation in DCM

With the available high-quality mRNA expression data (Fig EV2), we performed SNV-eQTL analyses and found altogether 3,917 transcripts associated with a SNV (FDR significance level ≤ 0.05) (Fig 3A). The five most significant variant-expression associations are linking rs3129888, rs2239802, rs3135390, rs7196, rs2395182 with Histocompatibility Antigen HLA-DR Alpha (HLA-DRA). Other highly significant SNVs link to, for example, the Endoplasmic Reticulum Aminopeptidase 2 (ERAP2) and Transmembrane Protein 117 (TMEM117) (Appendix Table S1).

Click here to expand this figure.

To test whether detected SVs also have a functional impact on the transcriptome as an intermediate cardiac phenotype, we next conducted a SV-eQTL analysis by performing correlation tests on the occurrence of all SV genotypes against all expressed transcript levels. Figure 3A is showing a genome-wide manhattan plot of all SV-eQTLs, with 2,652 reaching significance after correction for multiple testing (FDR P ≤ 0.05). When focusing on cis-regulation within 1 Mb distance (Sudmant et al, 2015), 75 genome-wide significant SV-eQTL links could be established (Appendix Table S2, Fig EV3) for 71 genes. Of those SV-eQTLs, ten expressed genes are altered due to the direct overlap with the SV, seven when considering linkage disequilibrium extension and another 54 in 1 Mb distance.

Click here to expand this figure.

Interestingly, most of the eQTLs are exclusively expressed in heart tissue (left panel) and not in blood (right panel), suggesting a predominant cardiac effect. In 49 out of the 75 SV-eQTLs (65.3%), a SNV-mediated effect is unlikely since no significant SNV-eQTLs are found in this region. Of the overlapping 760 eQTLs, the directional effect sizes of an additive linear model were comparable (Fig 3B). Also of note, we find 37 of the SV-eQTLs to affect non-coding RNAs, emphasizing the importance of the presumable regulative elements of untranslated regions in the human myocardial genome in DCM, which might affect other portions of the transcriptome. Such a complex, interacting system is supported by the differential effects observed for SV-eQTLs lying in enhancers, transcription factor binding sites, or non-coding RNAs (Fig 4A).

The allele frequencies of detected SVs and SV-eQTLs ranged from private (only detected in one individual) to common (frequency ≥ 5% of the cohort). Figure 4B exemplarily shows the effect of common SVs on the cardiac transcription; for example, an upregulation of the gene expression was found in the chromatin remodeler gene Alpha Thalassemia/Mental Retardation Syndrome X-Linked (ATRX), Keratin 1 (KRT1), involved in activation of the immune system and beta catenin signaling, and Serpin family C member 1 (SERPINC1), involved in blood coagulation. The increased amount of mRNA transcript for SERPINC1 detected in SV carriers resulted in higher protein expression when compared to non-SV carrier (Fig 4B), indicating the profound effect to be expected by this class of variation.

Other SVs directly affected the coding transcript either by falling in the exon or introns. Figure EV4 shows altered expression of Amidohydrolase Domain-Containing Protein 1 (AMDHD1) gene that is implicated in amino acid synthesis. As shown, the SV results in an intronic deletion, which directly affects the expression of the mutant allele. To investigate the quantitative effect of the linked SV on the transcription in the corresponding 1-Mb interval, we calculated the determination of coefficient R² per SV event, showing the proportion of differential expression that can be explained by the corresponding SV. For the SV-eQTL AMDHD1, 18% of the total expression variation in the genomic interval can be explained by this approach. Another locus significantly affected by an SV in multiple patients (56%) is harboring five SV-eQTLs (Fig 5A). The 37-kb large deletion, which is spanning the entire Glutathione S-Transferase Theta 2B (GSTT2B) gene and the first two exons of D-Dopachrome Tautomerase-Like (DDTL) gene, accounts for 13% of the total expression alterations in this locus (Fig 5B). The methylation intensities for CpG sites located in this deletion showed decreased signals for heterozygous and homozygous SV carriers, respectively (Fig 5C). Very interestingly, of the five significantly SV-linked genes, two (GSTT1 & 2) are dysregulated in a mouse model of transverse aortic constriction (TAC) induced heart failure (Fig 5D), at least indicating the involvement in adaptive or maladaptive cascades in the failing mammalian heart.

Click here to expand this figure.

While common SVs may have a regulatory effect on the transcriptome but theoretically only smaller effects on detrimental phenotypes, rare genetic variations bear the potential to significantly impact on human disease. Using SV-eQTL analysis for rare SVs involving gene regulatory regions, we found regulation of PRELID1P4, where a deletion resulted in an upregulation in two SV carriers (Fig 6A). For ZNF35 a complex event, a combination of an inversion and a duplication in two different DCM patients led to decreased transcript levels. Another complex SV (inversion and deletion) upregulated the expression of CAMK2A (Fig 6B), the gene for which the first knockout mouse was made (Silva et al, 1992). For validation, the genomic region of the complex SV (Fig 6D) was amplified and both single alleles (SV and non-SV) of a SV carrier were subjected to long-read nanopore sequencing, which confirmed the complex SV (Fig 6D). A more detailed analysis of the genetic context revealed that the CAMK2A upregulation on transcription level is exclusively driven by the expression of six C-terminal exons of the 19 exon spanning full-length CAMK2A (Fig 6B). Here, the SV explains 6% of the transcriptional variation of CAMK2A (Fig 6C) and 14% of the depicted genomic region.

Transcriptome-wide effects of SVs in the heart

To examine the global impact of structural variation on myocardial gene expression in humans, we next covered both cis and trans regulatory mechanisms. Of the total of 20,712 genes found to be expressed in the myocardium, 16,449 genes (79%) could be directly or indirectly linked to a SV event at least in one proband (nominal P ≤ 0.01). To determine to what extent those SV-eQTL genes contribute to the overall transcriptional variation found in the 20,712 genes, we applied linear modeling of attributed expression variation. By this approach, we inferred the degree of global SV-associated expressional variation to be 7.5% (Appendix Table S3). A similar picture is seen for the DCM-related genes, where SV-eQTLs were found to explain a fraction of 10.1% of the transcriptional variation.

Patients and study design

The characterization of samples and patient data has been approved by the ethics committee, medical faculty of Heidelberg, participants have given written informed consent, and Care4DCM project was conducted (Meder et al, 2017). Symptomatic DCM patients were consecutively, prospectively enrolled. A prerequisite for enrollment was leftover myocardial tissue from the routine diagnostic workup. For the exclusion of secondary causes of DCM, all patients underwent diagnostic coronary angiography, histopathology of myocardial biopsies, echocardiography, cMRI, comprehensive clinical phenotyping, and biomarker measurements. Patients with valvular or hypertensive heart disease, history of myocarditis, regular alcohol consumption, or cardio-toxic chemotherapy were also excluded.

Biopsy specimens were obtained from the apical part of the free left ventricular wall (LV) from DCM patients undergoing cardiac catheterization using a standardized protocol. Biopsies were rinsed with NaCl (0.9%) and immediately transferred and stored in liquid nitrogen until DNA or RNA was extracted. Total RNA was extracted from biopsies using the RNeasy kit according to the manufacturer's protocol (Qiagen, Germany). RNA purity and concentration were determined using the Bioanalyzer 2100 (Agilent Technologies, Berkshire, UK) with a Eukaryote Total RNA Pico assay chip.

Whole-genome sequencing

1 μg of total gDNA was sheared using the Covaris™ S220 system, applying two treatments of 60 s each (peak power = 140; duty factor = 10) with 200 cycles/burst. 500 ng of sheared gDNA was taken, and whole-genome libraries were prepared using TruSeq DNA sample preparation kit according to manufacturer's protocols (Illumina, San Diego, US). Sequencing was performed on an IlluminaHiSeq 2000, using TruSeq SBS Kit v3 and reading two times 100 bp for paired-end sequencing, on four lanes of a sequencing flowcell.

Demultiplexing of the raw sequencing reads and generation of the fastq files was done using CASAVA v.1.82. The raw reads were then mapped to the human reference genome (GRCh37/hg19) with the burrows-wheeler alignment tool (BWA v.0.7.5a) (Li & Durbin, 2009), and duplicate reads were marked (Picard-tools 1.56) (http://picard.sourceforge.net/). Next, we used the Genome-Analysis-Toolkit according to the recommended protocols for variant recalibration (v. 2.8-1-g932cd3a) and variant calling (v.3.3-0-g37228af) as described in the respective best-practices guidelines (https://www.broadinstitute.org/gatk/guide/best-practices) (DePristo et al, 2011). Structural variants were detected using paired-end and split-read algorithms. For this, aligned reads were processed with Delly (v0.6.3) (Rausch et al, 2012), excluding centromeric and telomeric regions of hg19 and employing an insert size cutoff of median + 9 times median absolute deviation for deletions. SVs were included in the analysis, if the median normalized read-depth ratio of heterozygous SV carriers compared to SV non-carriers was < 0.8 or > 1.3 for deletions or duplications, respectively. Only those inversions were kept in the analysis that had inversion-supporting paired ends for both breakpoints. Complex structural variants were also discovered by Delly applying two complementary paired-end clusters and a concurrent read-depth decrease or increase for an inversion and deletion (INVDEL) or inversions and duplication (INVDUP), respectively. Obtained SVs were re-genotyped in the l WGS data of the 1000 Genome Project the phase 3 (http://www.1000genomes.org/) to obtain their allele frequencies.

Expression analysis

From the 50 samples in the cohort, high-quality RNA-seq libraries from left ventricular RNA could be generated for 42 patients using TrueSeq RNA Sample Prep Kit (Illumina). Sequencing was performed 2 × 75 bp on a HiSeq2000 (Illumina) sequencer. Samples were sequenced to a median paired-end read count of 31.5 million (range: 3.1–99.9). Unstranded paired-end raw read files were mapped with STAR v2.4.1c (Dobin & Gingeras, 2015) using GRCh37/hg19 and the Gencode 19 gene model (http://www.gencodegenes.org/). Only uniquely mapped reads were counted into genes using subread's featureCounts program (Liao et al, 2014) (subread version 1.4.6.p1), and mapping percentages were median 87.8 (range: 23.0–90.9). Prior to statistical analyses, genes with very low expression levels (average reads ≤ 1, detected reads in less than 50% of the samples) were removed resulting in 20,712 significantly expressed genes. Count data were normalized by rlog normalization (Love et al, 2014), which is an improved method of the variance stabilization transformation (Anders & Huber, 2010) as recommended for eQTL by the original MatrixEQTL publication (Shabalin, 2012).

SV-eQTL

An eQTL analysis between SVs and gene expressions is performed on the 42 patients with high-quality transcriptome data from biopsy samples. MatrixEQTL (Shabalin, 2012) is used to correlate the 3,897 SV events and the expression profiles of 20,712 genes. To account for LD effect, SV spans are first extended using a twofold extension method as such. The SV is extended to the range from the furthest upstream base pair with LD R² > 0.5 to the furthest downstream base pair with R² > 0.5 to the SV. Then, the range is further extended to the next immediate recombination sites. LD and recombination site data for GRCh36 were first obtained from Hapmap.org and lifted over to GRCh37 using UCSC liftover tool.

For performing an eQTL, it is important to define an interval where a possible link is calculated. To estimate a genuine window between gene and extended SV, we chose to follow methods used in the latest release (phase3) of the 1000 Genomes croject, where an eQTL was considered within a 1 Mb range (Sudmant et al, 2015). Due to the genomic complexity of the human leukocyte antigen locus (HLA), located within the 6p21.3 region on the short arm of human chromosome 6, we excluded these loci from the SV-eQTL analysis.

Coefficient of determination (R²) analysis

From each SV-eQTL with P-value of 1% or less, a residual variance was computed. For a set of SV-linked genes, the residual variance of each gene was summed to generate a residual variance of the set of genes. Together with the total variance of the set, the coefficient of determination of the set of genes was then calculated accordingly, that is, 1 – ∑ (residual variances)/(total variance). The value represents then the proportion of total variance that the model explains. The coefficients of determination are meant to be descriptive, and hence, no associated statistical significances are calculated.

SNV-eQTL

An eQTL analysis between SNVs and gene expressions was performed on the 42 patients with high-quality transcriptome data. MatrixEQTL (Shabalin, 2012) is used to correlate the 14,720,818 SNV events and the autosomal expression profiles of 20,172 genes that are within 1 Mb base pairs of each other. To account for LD effect, SNV spans are first extended using a twofold extension method as such. The SNV is extended to the range from the furthest upstream base pair with LD R² > 0.5 to the furthest downstream base pair with R² > 0.5 to the SNV. Then, the range is further extended to the next immediate recombination sites. LD and recombination site data for GRCh36 were first obtained from Hapmap.org and lifted over to GRCh37 using UCSC liftover tool.

Verification of structural variants using high-density DNA methylation arrays

It has been shown that Illumina 450k methylation assay can be used to profile copy number alterations since the overall signal intensity of the methylated and unmethylated probes reflects the DNA amount and thereby copy number (Feber et al, 2014; Meder et al, 2017). To verify structural variants using MatrixEQTL, we correlated 633 SV events that covered at least one methylation locus measured on the Illumina 450k chip and the overall signal intensity for the respective methylation loci in blood and cardiac tissue using tissue as covariate. In case of duplications, we tested for increased signal and in case of deletions for reduced intensity in the presence of events. For each SV event, an aggregate significance level was obtained using the simes procedure (Rødland, 2006).

Polymerase chain reaction for SV validation

For technical validation of SVs, the PrimeSTAR GXL DNA Polymerase (TaKaRa Bio inc., Tokyo) was used, taking 10 ng of gDNA as template in 20 μl reaction volume using the following primers: XKR9 forward (5′-TTGTGTCCTAGACAGGCGAGTG-3′) and XKR9 reverse (5′-GCCAAATGAGGAGCTTGGCAAT-3′). TNKS2-AS1 forward (5′-TAGTACAGCTGCCCCTTGTGAC-3′) and TNKS2-AS1 reverse (5′-TGGCAGCCTGTTTAGATCCACT-3′). KBTBD11-OT1 forward (5′-ACAAGCGCTTTCAGGGGAAATG-3′) and KBTBD11-OT1 reverse (5′-TTTGGGTGAAGGCGTCTAACCA-3′). KRT1 forward (5′-GGGCGTGGATTCTTGTTCACAG-3′) and KRT1 reverse (5′-GTCTAACTTGGGGGTACGTGCT-3′).

Long-read nanopore sequencing

Genomic intervals were amplified using forward (5′-CCGTAAGTGCAATGCAATCCCT-3′) and reverse (5′-CTCCAGCAGGGTCTGAGGTTAC-3′) primer with Qiagen (Germany) Taq polymerase and separated in 1% agarose gel stained with Midori Green Advance (Nippon Genetics Europe, Germany) and extracted from the gel. 1.25 μg of the variant and wild-type fragment were taken for the library preparation according to the manufacturers protocol (GDE_1002_v1_revB_17Nov2015) with the SQK-MAP006 sequencing kit on a MinION Mk1. Obtained FAST5 sequences were processed with poRetools (Loman & Quinlan, 2014) and aligned with bwa-mem (0.7.10-r789).

Immunoblot analysis

Protein expression analysis of was performed on left ventricular material of independent samples, which were homogenized with ceramic beads (1.4 and 2.8 mm) in 400 μl lysis buffer (50 mM Tris–HCl, 120 mM NaCl, 5 mM EDTA, 0.1% NP-40, 1 mM DTT, 1 mM sodium metavanadate, 1 mM sodium fluoride, 0.2 mM PMSF, protease inhibitor (cOmplete tablet Roche cat# 04693116001) and phosphatase inhibitor (PhosSTOP Roche cat# 04906837001). Equal amounts of protein as tested by the Pan Cadherin control (abcam ab16505; dilution 1:1,000) were diluted in 4× Laemmli sample buffer (Bio-Rad cat# 161-0747) and separated on a 4–20% gradient gel (Bio-Rad, Germany, cat.# 4561094). Primary antibody for SERPINC1 (ThermoFisher Scientific, cat# PA5-13674; dilution 1:750) and secondary HRP-linked anti-rabbit IgG antibody (Cell Signaling Technology, cat#7074; dilution 1:2,500) were used. Pierce™ ECL Western Blotting Substrate (cat# 32209) was used for the detection of HRP.

Transverse aortic constriction

Institutionally available TAC data were used to investigate potential dysregulation of SV-eQTL homologous transcripts in mice with induced heart failure. TAC was performed as previously described (Volkers et al, 2013). Briefly, in 8-week-old male mice, the aorta between the innominate and the left common carotid arteries was ligated with a 7-0 polypropylene suture and a 27-gauge needle, which was removed after ligation. Before extubation and closing the chest, the pneumothorax was reduced. A sham procedure in which the aorta was not bended was also performed.

Statistical analyses and databases

Statistical analyses were carried out in R-3.2.2 (R Development Core Team, 2008). FDR correction of significance levels was performed using the Benjamini–Hochberg procedure (Benjamini & Hochberg, 1995). TFBSs employed in this study were annotated by http://genome.ucsc.edu/cgi-bin/hgTrackUi?db=hg19&g=wgEncodeRegTfbsClusteredV3. Enhancers used for annotation were obtained from FANTOM5 human permissive enhancers phase 1 and phase 2 (FANTOM Consortium and the RIKEN PMI and CLST (DGT) et al, 2014). All other genomic loci used for functional annotation were pulled from GENCODE (v19).

Data availability

The data are freely accessible (accession number CMS-SV-17; https://ccb-web.cs.uni-saarland.de/cms).