Whole gene sequencing identifies deep-intronic variants with potential functional impact in patients with hypertrophic cardiomyopathy

Background High throughput sequencing technologies have revolutionized the identification of mutations responsible for genetic diseases such as hypertrophic cardiomyopathy (HCM). However, approximately 50% of individuals with a clinical diagnosis of HCM have no causal mutation identified. This may be due to the presence of pathogenic mutations located deep within the introns, which are not detected by conventional sequencing analysis restricted to exons and exon-intron boundaries. Objective The aim of this study was to develop a whole-gene sequencing strategy to prioritize deep intronic variants that may play a role in HCM pathogenesis. Methods and results The full genomic DNA sequence of 26 genes previously associated with HCM was analysed in 16 unrelated patients. We identified likely pathogenic deep intronic variants in VCL, PRKAG2 and TTN genes. These variants, which are predicted to act through disruption of either splicing or transcription factor binding sites, are 3-fold more frequent in our cohort of probands than in normal European populations. Moreover, we found a patient that is compound heterozygous for a splice site mutation in MYBPC3 and the deep intronic VCL variant. Analysis of family members revealed that carriers of the MYBPC3 mutation alone do not manifest the disease, while family members that are compound heterozygous are clinically affected. Conclusion This study provides a framework for scrutinizing variation along the complete intronic sequence of HCM-associated genes and prioritizing candidates for mechanistic and functional analysis. Our data suggest that deep intronic variation contributes to HCM phenotype.


Introduction
Hypertrophic cardiomyopathy (HCM) is a genetic heart disease associated with sudden cardiac death and progressive heart failure.HCM is considered one of the most common genetic disorders, with an estimated prevalence of 1 in 500 people throughout the world [1].Recognition of the disease is critical for providing treatment and prevention strategies as well as triggering clinical and genetic surveillance of family members [2], [3].Mutation carriers may benefit from lifestyle and medical interventions that improve prognosis, whereas a negative genetic test can reassure individuals that are not at risk [3], [4].
Some rare inherited diseases may mimic the phenotypic and clinical features of sarcomere HCM, as defined by the presence of unexplained left ventricular hypertrophy.These conditions are referred to as HCM phenocopies and represent distinct disease entities with respect to inheritance, pathophysiology, natural history, extra-cardiac features, and management [3], [18].These disorders are not caused by sarcomeric mutations.The most prominent HCM phenocopies in adults include [3] Fabry disease, caused by mutations in the galactosidase-α gene (GLA); Danon disease, a lysosomal storage disease caused by mutations in the lysosomal-associated membrane protein 2 gene (LAMP2); and LVH associated with Wolff-Parkinson-White syndrome, caused by mutations in the regulatory subunit of adenosine monophosphate-activated protein kinase gene (PRKAG2) (Table 1).
In recent years, widespread availability of genetic testing has proved crucial not only to identify the sarcomeric mutations that cause HCM but also to distinguish disorders that can mimic HCM [2,4].However, despite the revolutionary increase in genetic testing capability introduced by next-generation sequencing [19], approximately 50% of individuals with a clinical diagnosis of HCM have no causal mutation identified [4,20], [21].
One possibility to explain why many individuals fail to be genetically diagnosed is the presence of deep-intronic mutations undetected by current clinical genetic testing approaches, which provide information restricted to the exons and exon-intron boundaries.Indeed, recent Genome Wide Association Studies have identified many single nucleotide variants located deep within introns with significant association to diseases [22,23].To date, over 180 deep intronic pathogenic variants located at least 100 bp from the nearest canonical splice site have been reported across 77 different disease genes [24].Most frequently, these mutations lead to pseudo-exon inclusion due to activation of non-canonical splice sites or changes in splicing regulatory elements [25,26].The more common mechanism involves a mutation that creates a novel donor splice site and activates a pre-existing non-canonical acceptor splice site, whereas more rarely the mutation creates a novel acceptor splice site.Alternatively, inclusion of cryptic exons can be induced by mutations that either create an enhancer sequence element or inactivate a repressive sequence [24].For example, a deep intronic mutation (c.639+919G>A) in the galactosidase alpha (GLA) gene, responsible for Fabry disease, disrupts an hnRNP A1 and hnRNP A2/B1-binding splicing silencer motif, thus allowing binding of U1 snRNP to an overlapping cryptic 5'ss that results in pseudo-exon inclusion [27].The appearance of a pseudoexon generally disrupts the reading frame introducing a premature termination codon that targets the mutant mRNA for degradation by nonsense mediated decay (NMD) [28], making mutations that result in abnormal splicing functionally equivalent to null or hypomorphic alleles.In other cases of genetic diseases caused by deep intronic variation, the mutation disrupts transcription regulatory motifs leading to decreased expression of the affected gene [24].In this study, we used targeted high throughput sequencing and computational approaches to identify deep intronic variants that may contribute to HCM phenotype.

Patients
The study population comprised 16 unrelated consecutively evaluated patients (8 males, 8 females) referred to the Cardiology Department at University Hospital Santa Maria.For all probands, the personal and family history, physical examination, ECG and echocardiography were consistent with a diagnosis of HCM according to international criteria [3].Patients were genetically tested at a mean age of 49 years.In addition, family members of two selected probands were clinically and genetically tested.Before blood collection, all patients and relatives provided written informed consent for DNA analysis and received genetic counselling in accordance with guidelines [3].DNA samples used in this study were residual after conventional diagnostic screening by targeted exome and Sanger sequencing.The project was approved by the Lisbon Academic Medical Center Ethics Committee.

Targeted gene enrichment and sequencing
Blood samples (5-8 mL) were collected into EDTA tubes at routine clinic visits, and DNA was isolated from peripheral blood lymphocytes using standard methods.The study was designed to screen the full genomic DNA sequence of 26 genes indicated in Table 1.These genes are included in many commercially available testing panels.A capture library was designed using SureSelect (Agilent) and target regions were sequenced (paired-end) on an Illumina HiSeq platform with 30-97 base read length.Highly repetitive sequences were excluded.Sample preparation was carried out as recommended by the manufacturer.Relatives were genotyped for selected variants by Sanger sequencing.

Bioinformatic data analysis
Raw sequencing paired-end reads (in.fastq format) were aligned using BWA software (version 0.7.12) [29] on the human reference genome (GRCh37) using quality score calibration and Illumina adapter trimming.Following the exclusion of duplicate reads using Picard MarkDuplicates tool (version 1.96) (http://broadinstitute.github.io/picard/),regions around insertiondeletions (indels) were realigned and each base quality score was recalibrated.For variant calling, we used four distinct tools: GATK-UnifiedGenotyper (version 3.4-46) [30] and SAMtools mpileup (version 1.2) [31], which use alignment-based approaches, and GATK-HaplotypeCaller (version 3.4-46) [30] and FreeBayes (version 0.9.21.26) [32], which use haplotype-based approaches.By comparing the performance of each tool against a standard reference (NA12878, published by Genome in a Bottle consortium [33], we observed a concordance of ~85% (S1 Fig) .To take advantage of the strengths of the different tools, we selected variants that were independently called by at least two of them.This strategy showed a better sensitivity (~97%) and precision (~98%) compared to analysis using a single tool (S1 Fig) .Variants that were independently selected by at least two tools, and presented a read depth of 20 or more in the targeted genes, were annotated with ANNOVAR [34] (Fig 1).
As deep intronic mutations may result in altered gene expression through either cryptic splicing or disruption of transcription regulatory motifs [24], we investigated whether the identified intronic variants may disrupt transcription factor binding sites (TFBS).We used available tracks in the UCSC genome browser [49][50][51] (Transcription Factor ChIP-seq Uniform Peaks from ENCODE/Analysis and HMR Conserved Transcription Factor Binding Sites), focusing on TFBS predicted to be targets for transcription factors that have been implicated in pathways related to cardiac regulation, development or pathophysiology.
Variant frequency was determined using the allele frequency estimates from the 1000 genomes project [52] and gnomAD [53] databases (accessed on June 2017).

Quality of sequencing data
Analysis of sequencing data yielded an average of 96.64% confidently mapped reads per gene.For 69% of the targeted genes the average read depth was above 200, and for the remaining genes the average read depth ranged between 130 and 200 (S2A Fig) .The average read depth was slightly lower over noncoding regions (S2B Fig) .The average percentage of covered base pairs was higher than 90 for 85% of the genes, and the lowest coverage was 76% for both coding and noncoding regions (S2 Fig) .Following alignment to the reference genome (GRCh37) and variant calling, we removed variants that were off-target, or had an average read depth below 20 (Fig 1).Single nucleotide substitutions and insertions or deletions of a few bases were identified and considered for further analysis.

Spectrum of exonic and splice site variants
Previously described disease-causing variants in the MYBPC3 gene were detected in 3 patients (Table 2).Rare variants classified in the NCBI ClinVar database and according to the ACMG guidelines as of uncertain significance were additionally detected in the TNNT2, MYBPC3, TTN, TPM1 and MYH6 genes; all these scored as likely pathogenic according to multiple in silico prediction tools (Table 2).Noteworthy, one of the patients harboured, in addition to MYH6 variant rs140596256, a novel variant in the GLA gene that is not listed in online databases but is predicted to be pathogenic by multiple prediction tools (Table 2).

Assessment of deep intronic variants
The noncoding variants were prioritized using GWAVA [44], CADD [45], SPIDEX [22], and Genomiser [46].The genome-wide annotation of variants (GWAVA) is a computational approach that integrates a wide range of available genomic and epigenomic annotations to predict the functional impact of variants.GWAVA results are in the range 0-1, with higher values indicating variants predicted as more likely to be functional.Variants with a GWAVA score above 0.5 were classified as functional, as in previous studies [44].The Combined Annotation-Dependent Depletion (CADD) method provides a metric (C score) for deleteriousness, a property that strongly correlates with functionality and pathogenicity [45].Variants at the top 10% of deleteriousness are assigned a C score of 10, whereas variants at the top 1% are assigned a C score of 20.Variants with C score greater than 15 were selected, as previously described [54].SPIDEX is a computational model that uses the Percentage of Spliced-In (PSI) metric to evaluate whether a certain splicing isoform is more enriched under the presence/absence of a given variant.SPIDEX scores higher than 5 predict that the variant affects RNA splicing [22].The Genomiser framework combines a machine learning method and an integrative algorithm for ranking noncoding variants in whole-genome sequence data [46].Genomiser results are in the range 0-1, with values higher than 0.6-0.9indicating variants more likely to be pathogenic [46].We found that all noncoding variants with higher scores for the different prediction metrics correspond to single nucleotide substitutions, the vast majority of which are located within introns (Fig 1).The position of each variant relative to the nearest canonical splice site ranged between 20 and 2000 nucleotides (Fig 2A).A comparison of variants prioritized as likely pathogenic by the different tools resulted in the identification of 6 variants that scored with high values using GWAVA, CADD and Genomiser metrics (Fig 2B).These include two variants in VCL, two variants in TTN, one variant in ACTC1 and one variant in PRKAG2 (Table 3).Analysis of allele frequency data available for European populations in the 1000 Genomes project [52] and gnomAD [53] databases reveals that two of these variants are more frequent in the patient population than in healthy individuals (Table 3).Namely, the VCL variant c.499+-367T>C (rs113195070) was detected in 6 probands, corresponding to an allele frequency of 19% in the patient population.This contrasts with a frequency of 6-7% in control populations.Similarly, the PRKAG2 variant c.1234-317T>G was present with an allele frequency of 3% in the patient population contrasting with a frequency of 0.1-0.3% in healthy individuals.Such specific enrichment of certain deep intronic variants in the patient population suggests that these may be contributing to the disease phenotype.
We focused on the VCL variant c.499+367T>C.We found that two of the probands were compound heterozygous for this variant and a MYBPC3 mutation previously described as disease-causing.Genotyping of family members of proband #6 showed that the dual presence of the MYBPC3 splice site mutation (c.1227-13G>A) and the VCL variant is associated with the manifestation of the phenotype in the proband (I-2) and his son (II-1), both diagnosed in their 40s (Fig 3).The other children of the proband (II-2 and II-3), while carrying the MYBPC3 mutation, did not develop signs of cardiomyopathy when assessed at a similar age.This suggests a possible modifier effect of the VCL variant, since the presence of the MYBPC3 mutation alone is not sufficient for the phenotype to be manifested.Analysis of this family further suggested that the VCL variant on its own is not sufficient to cause disease.Genotyping of family members of proband #15 (Fig 3) indicated that presence of the MYBPC3 missense mutation (p.Arg495Gln) in the absence of the VCL variant appears sufficient to cause the phenotype, even at a pediatric age (III-1).Clinical characteristics for the two families are detailed in Table 4.
We found that, based on Chip-seq experiments [49,50], the deep intronic VCL variant enriched in the patient population (rs113195070) localizes in a region associated with FOS, JUN and EP300 (Fig 4A).A deep intronic variant in the PRKAG2 gene (rs114394151) prioritized by CADD and Genomiser is enriched in the patient population (Table 3) and localizes in a region associated with FOS and JUN (Fig 4B).FOS and JUN transcription factors are thought to be among the first set of genes to be expressed in the context of pathological cardiac hypertrophy [55], and EP300 has been associated with cardiomyocyte enlargement [56].
Finally, using the SPIDEX tool we identified a variant in the TTN gene (rs72650063) that occurs with a frequency of 2% in control European populations and is present in two probands, corresponding to an allele frequency of 6% in the patient population (Table 3).This variant is predicted by Human Splicing Finder to disrupt binding of splicing factor SRSF1 (Fig 4C) [57].Another variant in the TTN gene (rs142156368) appears highly enriched in our cohort relative to the general population (Table 3).
No potential association of the candidate deep intronic variants with cardiac diseases identified through GWAS was found.

Discussion
Motivated by the clinical heterogeneity of HCM and the lack of a conclusive genetic diagnosis in approximately 50% of the patients [4,20], [21], we hypothesized that genetic variation within deep intronic regions of sarcomere and sarcomere-related genes contributes to the disease mechanism.Using a targeted high-throughput sequencing strategy, we did a comprehensive screening of 26 genes in a cohort of 16 unrelated HCM patients using recently developed computational models to assess variants.

Tools
In agreement with previous studies on larger HCM patient cohorts [61,63], we identified missense and stop-codon mutations, as well as splice site mutations (Table 2).While missense mutations may code for a pathogenic protein that would cause HCM by a gain-of-function mechanism, stop-codon and splice site mutations are more likely to act through a loss-of-function mechanism due to reduced levels of the normal protein.Demonstrating that particular mutations act through a dominant gain-of-function mechanism and others cause loss-of-function would be critical to understand phenotype-genotype correlations.In this regard, we analysed the family members of two probands with either a missense or a splice site mutation in MYBPC3 (Fig 3) and found that all carriers of the missense mutation (family #15) were clinically affected, as expected for a dominant gain-of-function mechanism.In contrast, in the other family, two individuals are carriers of the splice site mutation (family #6) and do not manifest the disease, consistent with a loss-of-function model.Possibly, the MYBPC3 splice site mutation present in family members of proband #6 results in mis-splicing by either decreasing the specificity or fidelity of splice site selection or activating cryptic splice sites that are normally not used.Abnormal splicing often results in a frameshift and consequent introduction of premature termination codons (PTCs), which trigger degradation of the mRNA by nonsense-mediated decay [28].Thus, this mutation can be functionally equivalent to a null or hypomorphic allele associated with loss-of-function of the protein.In contrast, the MYBPC3 missense mutation present in family members of proband #15 presumably leads to an abnormal protein containing an amino acid substitution that may cause a gain-of-function phenotype.
We further show that both probands are compound heterozygous for the missense or splice site MYBPC3 mutation and a deep intronic variant in VCL.This variant (rs113195070) is predicted to be pathogenic based on three independent computational tools, GWAVA, CADD and Genomiser (Table 3).Moreover, it is 3-fold more frequent in our cohort of probands than in normal European populations (Table 3), further suggesting a direct contribution to the   disease phenotype.The variant, which consists of a single nucleotide substitution located at position 367 from the nearest canonical splice site (c.499+367T>C), can potentially disrupt the binding of transcription factors that have been reported as implicated in pathways related to cardiac regulation, development or pathophysiology such as FOS, JUN and EP300 [55], [56].By interfering with the binding of transcription regulatory factors, the variant is expected to alter the transcription rate of the VCL gene.Consistent with this view, sequence elements located within introns of large human genes have been shown to act as transcriptional enhancers [64], and a recent study reported an IRF4 gene variant located in intron 4 that strongly affects IRF4 transcription through disruption of an enhancer element [65].Analysis of family #6 reveals that the presence of the VCL variant or the MYBPC3 mutation in isolation is not sufficient to cause disease phenotype.Indeed, the two clinically affected individuals in this family are compound heterozygous for the VCL variant and the MYBPC3 splice site mutation (Fig 3), suggesting that the combination of the two mutations triggers the disease.A loss-of-function mechanism for the MYBPC3 mutation could explain why in family #6 only compound heterozygous members manifest the disease, whereas the presence of the heterozygous gain-of-function mutation in family #15 would be sufficient to cause disease.Complex genotypes, including individuals that carry 2 or more variants in the same or different sarcomere-related genes, have been reported in 8% of HCM patients [58], and there is evidence indicating that patients with complex genotype and multiple simultaneous mutations may have more severe or early disease expression [66].However, complex genotype-phenotype correlations focusing specifically on carriers of splice site mutations remain to be investigated.
We further identified two single nucleotide substitutions in the titin gene (rs142156368 and rs72650063) that are 3 to 6-fold more frequent in our cohort of probands than in normal European populations (Table 3).These variants are located in the PEVK domain that plays a role in extensibility of the sarcomere and contractility of the titin protein [67,68].Titin is prone to extensive alternative splicing that can change its size and its elastic/stiffness properties; associations have been established between the ratio of expression levels for the main cardiac isoforms (N2BA and N2B) and genetic and non-genetic forms of cardiac diseases [69,70].If these variants do interfere with titin splicing, as predicted by the SPIDEX computational model, they are likely to contribute to HCM phenotype, particularly in combination with other HCM-associated alleles.Supporting this view, titin-truncating splicing isoforms, which are encountered in approximately 1% of the general population, are sufficient to induce molecular and physiological effects on the heart [71].
In conclusion, this study provides a framework for scrutinizing variation along the complete sequence of HCM-associated genes and prioritizing candidates for further analysis.Our data suggest that deep intronic variation contributes to HCM phenotype.Translation of genetic information found in an individual to clinical decision taking requires a precise understanding of the molecular mechanisms underlying the disease phenotype.To date, mechanistic and functional studies have been largely restricted to animal models in part due to difficulties in obtaining human tissue from patients.However, the recent emergence of patient-derived induced pluripotent stem cells (iPSCs) that can be differentiated into functional cardiomyocytes recapitulating HCM-specific characteristics [72,73] holds great promise as an exciting new approach to study how gene mutations relate to clinical outcomes and might be applied to test our hypothesis-generating data.

Fig 1 .
Fig 1. Flowchart of noncoding data analysis.Variants that were independently selected by at least two tools and presented a read depth of 20 or more in the 26 targeted genes were annotated with ANNOVAR.In silico predictions were carried out for noncoding variants that were not classified as either benign/likely benign or pathogenic/likely pathogenic in NCBI ClinVar.All variants with scores above the indicated threshold were single nucleotide substitutions (SNVs).GWAVA: Genome-Wide Annotation of Variants.CADD: Combined Annotation Dependent Depletion.SPIDEX: Splicing Index.R.S.: Region score.PHRED: phred quality score.dPSI: percent of spliced in.https://doi.org/10.1371/journal.pone.0182946.g001

Fig 2 .
Fig 2. Assessment of intronic variants.A) Schematic diagram depicting the position of intronic variants prioritized by each prediction tool.B) Venn diagram illustrating intronic variants that are simultaneously prioritized by multiple prediction tools.https://doi.org/10.1371/journal.pone.0182946.g002

Fig 4 .
Fig 4. Variants located at binding sites for transcription and splicing factors.A) The VCL variant c.499+367T>C (rs113195070) is located at a binding site for transcription factors FOS, JUN and EP300.B) The PRKAG2 variant c.115-30242C>T (rs114394151) is located at a binding site for transcription factors FOS and JUN.C) The TTN variant c.32077 +31C>G (rs72650063) is predicted to disrupt the binding site of splicing factor SRSF1. https://doi.org/10.1371/journal.pone.0182946.g004

Table 1 . Name of the genes analyzed, Ensembl accession number, and chromosomal position. Sarcomere genes most frequently associated with HCM Gene Ensembl ID Chromosome
β-Myosin heavy chain (thick filament)

Table 2 . Putative HCM-causing variants located in exons and exon-intron boundaries.
VUS, variant of uncertain significance.ACMG, American College of Medical Genetics and Genomics.