Discovering miRNA Regulatory Networks in Holt–Oram Syndrome Using a Zebrafish Model

1. Introduction

MicroRNAs (miRNAs) are small non-coding RNAs of about 20–23 nucleotides that play an essential role in a variety of biological important pathways from development and physiology to diseases such as cancer (Chen and Rajewsky, 2007; Small and Olson, 2011). miRNAs are mostly known to function by targeting complementary sequences in mRNA transcripts, usually in the 3′ untranslated region (3′ UTR) and so inhibiting the translation and altering the stability of mRNA (Bartel, 2004; Yates et al., 2013). The identification and validation of miRNA–mRNA interactions is fundamental for discerning the role of miRNAs in the complex context of regulatory networks. However, since the miRNA binding is mostly not a perfect one-to-one match with the complementary target sites, it is difficult to predict miRNA targets. Consequently, several computational methods and tools have been developed in the last years (Yue et al., 2009; Peterson et al., 2014). They encompass a range of different computational approaches, from the modeling of physical interactions, exploiting common features like seed match, conservation, free energy, and site accessibility to the incorporation of less common features extracted through machine learning techniques. Computational methods predict hundreds of thousands target mRNAs per miRNA, instead the number of experimentally validated targets is very low. One possibility to reduce the false positive rate is to combine high-throughput experimental data with sequence-based predictions (Huang et al., 2007; Muniategui et al., 2013). Although, this approach does not allow to identify miRNA targets that are repressed exclusively at the translational level. Since many miRNAs cause degradation of their targets (Baek et al., 2008; Hendrickson et al., 2009; Guo et al., 2010; Subtelny et al., 2014), the integration of expression profiles has been proposed to be an effective strategy to discover true miRNA–target interactions (Gennarino et al., 2009; Nazarov et al., 2013; Albert et al., 2014).

In this work, we used expression values of miRNAs and mRNAs obtained with high-throughput technologies to study complex regulatory networks altered in the Holt–Oram syndrome (HOS). HOS is a rare autosomal congenital disease characterized by cardiac and upper limb malformations (Basson et al., 1997). Mutations in the T-box gene Tbx5, which encodes a key transcription factor for vertebrate heart development, are responsible for HOS (Horb and Thomsen, 1999; Goetz et al., 2006). Family members with identical Tbx5 mutations can display large variations in malformation severity and HOS penetrance. This peculiar characteristic of HOS can be explained by the observation that Tbx5 is part of an extremely complex regulatory network. Due to the high number of messenger RNAs that are targeted by one miRNA, miRNAs are the best candidates to orchestrate the downstream regulation of Tbx5 gene expression in embryonic heart development. We have recently shown that miRNAs are crucial components of this network (Chiavacci et al., 2012, 2015). In fact, we proved that in mouse cardiac cells and zebrafish embryos, Tbx5 is able to regulate several miRNAs and, in particular, miR-218 and miR-19 (Chiavacci et al., 2012, 2015). The dysregulation of both miRNAs has a severe impact on heart development, affecting early heart morphogenesis.

As a model system, zebrafish has been extensively used for studying early vertebrate development (Kimmel et al., 1995; Yao et al., 2014) over the last 20 years. In particular, the HOS model called heartstring (hts) mutant has been well established in zebrafish, and it recapitulates almost completely the HOS characteristics. Furthermore, the zebrafish hts mutant can be easily replicated with the injection in zebrafish eggs of a specific Tbx5 morpholino (small antisense ribo-oligonucleotides, which blocks target translation) (Garrity et al., 2002).

Here, we propose an integrative approach, which uses experimental data from zebrafish HOS model system and computational methods for investigating in vivo complex regulatory networks perturbed in this pathology across two different stages of zebrafish development, 24 and 48 hpf. Those two stages were chosen since they mark fundamental steps in heart development. By 24 hpf, the migration phase is concluded, and the heart tube lies along the anteroposterior axis of the embryo with the atrial end to the left of the midline. By 48 hpf, the heart development is substantially completed: the heart terminated the looping phase and functional valves are formed (Kimmel et al., 1995; Yao et al., 2014). We show that it is possible to use data integration methods for studying rare diseases, providing significant insight into biological processes, and identifying new potential markers and drug targets of clinical interest.

2. Materials and Methods

2.1. Embryos Injection

The zebrafish line used in this study is the wild-type AB strain, the animals were raised and maintained under standard laboratory conditions (Westerfield, 1993). To silence the zebrafish gene, Tbx5a we used the antisense morpholino oligonucleotide MO-Tbx5a against the translational start site of the gene, the sequence of MO-Tbx5a was 5′-GAA AGG TGT CTT CAC TGT CCG CCA T-3′ (Garrity et al., 2002). The sequence of the control morpholino, MO-Ct, was 5′-CCT CTT ACC TCA GTT ACA ATT TAT A-3′. All morpholinos were supplied by Gene Tools LLC. Zebrafish morpholinos were injected into the yolk of 1-cell stage embryos with a constant injection volume, ~1 nl, using a microinjector (Tritech Research, Los Angeles, CA, USA). Zebrafish eggs were injected with 1.5 ng of MO-Tbx5a or 1.5 ng of MO-Ct, and embryos were collected at 24 and 48 hpf.

2.2. RNA Extraction, Library Preparation, Sequencing, and Microarray

For high-throughput DNA sequencing, total RNA was extracted from batch of n = 50 zebrafish embryos. The library preparation was done as described in (Baumgart et al., 2012). In detail, 500 ng of total RNA was used as input material. Library preparation was done using the TruSeq Small RNA Sample Prep (Illumina). The purified libraries were quantified on the Agilent DNA 1000 chip, diluted to 10 nM and subjected to sequencing-by-synthesis on Illumina HiSeq 2000 producing single-end 51 bp read length. Two independent batches of embryos were used for MO-Tbx5a and MO-Ct at 24 hpf, one for both condition at 48 hpf.

To measure mRNA expression, the Agilent Low Input Quick Amp labeling kit was used to retrotranscribe into the cDNA (from 200 ng total RNA), amplify, and incorporate the cyanine 3-labeled CTP (cRNA). The method uses the T7 RNA polymerase, which simultaneously amplifies and incorporates cyanine 3-labeled CTP. The fluorescent cRNA was purified and hybridized to the Agilent Zebrafish V3 Gene Expression Microarray 4 × 44, according to the manufacturer protocol. Three independent batches of embryos were assessed for MO-Tbx5a, while two for MO-Ct both at 24 and 48 hpf stages. Resulting images were quantified and text files containing raw values were analyzed.

2.3. Analysis of Sequencing and Microarray Data

Raw sequences were obtained and de-multiplexed using the Illumina pipeline CASAVA v1.8.2 FastQC v0.10.1¹, which was used for quality check, and primary reads were initially trimmed off to remove adapters sequence using Cutadapt v.1.2.1 (Martin, 2011). Employing FASTX_Toolkit (0.0.13.1), the reads with N calls were discarded. Remaining high quality reads, with a minimum length of 17 bp and a maximum 38 bp after clipping, were clustered for unique hits and mapped to zebrafish pre-miRNA sequences present into the mirBase (release 20) employing miRExpress (v2.1.3; Wang et al., 2009). We allowed 95% of sequence identity between read and reference sequence and a length tolerance range of 4 bp for mapping. miRNAs expression profiles were built by calculating the sum of read counts for each miRNA according to the alignment criteria. Differential expression analysis of miRNAs identified by miRExpress was performed using Bioconductor’s package DESeq (Anders and Huber, 2010). The reads count, used as measure of miRNAs quantification, was first normalized by library size factors to a common scale. The analysis was then performed and p-values were estimated using a negative binomial distribution model and local regression to estimate the relationship between the dispersion and the mean of each miRNAs. Raw p-values were finally adjusted for multiple testing using the Benjamini and Hochberg (1995) procedure controlling the false discovery rate (FDR). miRNAs with an adjusted p-value <0.05 were considered to be differentially expressed.

For microarrays, pre-processing of the data included background correction using a normal-exponential convolution model (offset = 16) (Ritchie et al., 2007) and cyclic loess normalization (Ballman et al., 2004) implemented in Limma package v.3.14.4 (Smyth, 2004). Low-expressed probes were filtered out keeping probes that are at least 10% brighter than the 95th percentile of the negative controls on at least 2 arrays. The Agilent Single channel Expression Microarray 4 × 44K for Zebrafish contains 39344 probes, 39162 are unique. For 35073 of them, we were able to retrieve gene accession ids corresponding to 21956 unique gene IDs. The linear modeling approach and empirical Bayes statistics implemented by Limma were used for assessing differentially expression. Finally, p-values were adjusted for multiple testing by means of the Benjamini and Hochberg method to control the false discovery rate. Genes with FDR less than 0.05 and fold change (FC) higher than 1.3 were selected for downstream analysis.

All statistical analyses were conducted using R and available Bioconductor packages.²

2.4. Integrated Analyses of Zebrafish miRNA and mRNA Expression Profiles

In order to discover miRNA-target pairs involved in HOS, we combined inverse correlations between miRNA and mRNA expression for improving in silico microRNA target predictions (see Figure 1). We selected the significant differential expressed miRNAs and mRNAs and performed target prediction analysis. Since miRNAs act at the post-transcriptional level downregulating their targets binding on the 3′-UTR of mRNAs, in this study, we focused our attention on these sequences that we retrieved from the UCSC Table Browser.³ We predicted miRNA target sites in the 3′-UTR using TargetScan Fish 6.2 (Lewis et al., 2005) and Pita (Kertesz et al., 2007) algorithms and then selected the consensus. Finally, we extracted the inversely correlated interactions (to reflect the typical miRNA–mRNA relationship) obtaining the final miRNA-target list.

FIGURE 1

Figure 1. Overview of the analytical workflow used in the study to identify inversely correlated putative target genes and to build altered regulatory networks in HOS.

3. Results

In the following sections, we detail the expression profiles of both miRNAs and annotated genes, which resulted altered by Tbx5a depletion during early zebrafish developmental stages (24 and 48 hpf). Small RNAseq and microarray analysis were performed to generate, respectively, miRNA and mRNA profiles. Moreover, we describe the main results obtained by integrating experimental data with computational methods to investigate in vivo regulatory networks modified by Tbx5 dosage alteration.

3.1. Sequencing and Annotation of miRNAs Modulated by Tbx5a at 24 and 48 hpf

In order to assess miRNAs expression modulation in zebrafish embryos after Tbx5a depletion, we conducted massive parallel sequencing experiments producing between 12.7 and 25.1 million total sequencing reads were obtained for each given library (16.8 mean) and this ranged from 10.5 to 20.2 million reads of 17–38 length after adapter trimming. On average, around 760 thousands of reads mapped to zebrafish miRNAs, annotated in miRBase v20 identifying 367 mature miRNAs on average per sample (see Table S2 in Supplementary Material). Among them, 19 miRNAs resulted to be significantly modulated at 24 hpf and 33 at 48 hpf (Table S3 and S4 in Supplementary Material respectively). We selected the most variable miRNAs, in terms of expression fold-change between Tbx5 and Ct morphants for downstream analysis: miR-34a, miR-10d-5p, miR-30a, miR-210-3p, and miR-5p at 24 hpf, miR-34a, miR-462, miR-146a, miR-21, miR-7b, and miR-190b at 48 hpf (Table 1). Differently from experiments reported in our previous work (Chiavacci et al., 2015), the downregulation of miR-19a at 48 hpf was not significant. However, this miRNA was included in the list of miRNAs modulated by Tbx5 because: (1) Q-RT PCR analysis performed in four different sets of experiments confirmed miR-19a downregulation (Figure 2B), (2) this downregulation was clearly supported by physiological data and by in situ hybridization experiments already presented (Chiavacci et al., 2015). Besides miR-19a-3p, other seven differentially modulated miRNAs were measured by quantitative RT-PCR and fold-changes were compared in Figures 2A,B for 24 and 48 hpf, respectively. All modulation was confirmed except for miR-210-5p, which resulted not significative.

TABLE 1

Table 1. Selected differentially expressed miRNAs at 24 and 48 hpf.

FIGURE 2

Figure 2. Validation of small RNA seq profiling and array-based gene expression profiles by Quantitative RT-PCR. (A) Sequencing and corresponding Q-RT PCR expressions of eight of the miRNAs reported in Table 1 and identified as differentially modulated in MO-Tbx5 vs. MO-Ct embryos at 24 hpf and at 48 hpf are reported. Values are expressed as fold change of MO-Tbx5 relative to MO-Ct. For Q-RT PCR, values are normalized on U6 expression. As pointed out in the results, miR-19a RNA-seq value is from Chiavacci et al. (2015). miR-210 is the 5p isoform. (B) Microarray and corresponding qRT-PCR expressions of eight genes showing differential expression in MO-Tbx5 compared to MO-Ct embryos at 24 and 48 hpf are reported. Values are expressed as fold change of MO-Tbx5 relative to MO-Ct. For Q-RT PCR, values are normalized on EF1, beta actin, and 18S expression. The values reported in the Q-RT PCR analysis are the mean of at least three independent microinjection experiments, t-test was used for statistical analysis: *p < 0.05, **p < 0.01, and ***p < 0.001.

3.2. Tbx5 Sensitive Genes in Early Developmental Stages of Zebrafish

To characterize the gene expression profiles at 24 and 48 hpf of zebrafish development and evaluate how altered Tbx5 dosage influences the genome-wide transcription, we measured mRNAs using expression microarray technology (see Materials and Methods for details). mRNAs were extracted from zebrafish embryos injected with MO-Tbx5a or MO-Ct and collected at 24 and 48 hpf. Using an absolute FC cut-off of 1.3 and an adjusted p-value of 0.05, we identified 7100 differentially modulated genes after TBx5a silencing at 24hpf, while 2276 genes at 48 hpf. The magnitude of differential expression was formally tested to be biologically significant using the t-test relative to a threshold (TREAT) method (McCarthy and Smyth, 2009) implemented in Limma. The complete lists are available as Supplementary Material. Validation by relative Q-RT PCR was performed for some of the genes taking into consideration the relevance for the cardiac context. Q-RT PCR analysis confirmed the microarray data (Figure 2A for 24 hpf and Figure 2B for 48 hpf).

To highlight most relevant functional categories among identified modulated genes, we performed a Functional Annotation Clustering using The Database for Annotation, Visualization and Integrated Discovery (DAVID) tool (Da Wei Huang and Lempicki, 2008). The Functional Annotation Clustering integrates the Kappa statistics to measure common genes between two annotations (e.g., ontological terms), and the fuzzy heuristic clustering to classify the groups of similar annotations according to kappa values. The resulting groups have similar biological meaning due to share similar gene members. We considered KEGG pathways and Gene Ontology terms performing the enrichment analysis of downregulated and upregulated genes (separately) at 24 and 48 hpf. Clusters that resulted significant (p < 0.05) are reported in Table 2. Interestingly, at 24 hpf, the top scoring functional categories contained genes that are involved in cell adhesion and ion binding. It is well known that morphogenesis requires specific cell adhesion molecules that are expressed in a precise developmental time, and the altered gene expression leads to heart defects (Buck et al., 1993; Kwee et al., 1995). In accordance with this observation and with our results, genes annotated with the term homophilic cell adhesion were also identified as significantly upregulated following inhibition of Tbx5a in a microarray-based expression profile performed in 56 hpf Tbx5 morphant zebrafish embryos (Mosimann et al., 2015). Furthermore, we found that the majority of modulated genes consisted of the cation and ion binding categories. In this context, the cation calcium has an important role in heart development, functions, and diseases (Arnolds et al., 2012; Crocini et al., 2014). The Calcium Binding Proteins (CaBPs) share a very similar domain organization with Calmodulin (CaM) and have been shown to have coevolved in vertebrate animals (McCue et al., 2010). The CaBPs have an important role during the development and in several diseases, such as Diastolic dysfunction that is characterized by slow or incomplete relaxation of the ventricles during diastole and is an important player in heart failure pathophysiology (Asp et al., 2013).

TABLE 2

Table 2. Most significant categories from functional annotation clustering analysis of the deregulated transcripts were reported.

3.3. Changes in miRNA Expression and Integration with mRNA Profile Identify Potential miRNA–mRNA Target Pairs Involved in HOS

In this study, we integrated two target prediction algorithms, TargetScan and Pita, with miRNA and gene expression data to refine in silico predictions and reduce the number of false positive interactions. The resulted miRNA-target pairs consisted of 122 potential targets at 24 hpf for upregulated miRNAs and 372 potential targets for downregulated miRNAs (complete lists are in Table S5 and S6 in Supplementary Material). At 48 hpf, we discovered 142 potential targets for upregulated miRNAs and 162 for downregulated miRNAs (see Table S7 and S8 in Supplementary Material). Among them, several miRNA–mRNA interactions involved genes that are known to be connected to heart development or cardiac functions published in previous works (Table S9 and S10 in Supplementary Material). We summarized these finding in Figures 3 and 4 and explore most interesting functional relations in the next section.

FIGURE 3

Figure 3. Regulatory network altered in zebrafish HOS model at 24 hpf. Potential interactions involving Tbx5, transcriptional factors, miRNAs, and their targets are shown together with the functional impact in heart development.

FIGURE 4

Figure 4. Regulatory network altered in zebrafish HOS model at 48 hpf. Potential interactions involving Tbx5, miRNAs, and their targets are shown together with the functional impact in heart development.

4. Discussion

Tbx5 is a crucial transcription factor in heart development. In HOS murine model, it has been shown that even small alterations of this gene cause modulation of hundreds of genes (Mori et al., 2006). It has been suggested that the strong impact that Tbx5 has on gene expression is mainly the result of its ability to modulate other regulators, such as different transcription factors, in a very complex regulatory network. Our previous studies suggest that Tbx5 affects the embryo development by modulating also miRNAs (Chiavacci et al., 2012). Moreover, the fact that miR-19a replacement is able to partially rescue fins and cardiac defects in a zebrafish model of HOS, strongly supports the importance of miRNAs in Tbx5 regulatory circuits (Chiavacci et al., 2015).

In this study, we analyzed miRNA and mRNA expression profiles at two fundamental time points (24 and 48 hpf) of zebrafish embryos development after depletion of Tbx5 and compared them to the wild-type ones. We employed expression data to improve miRNA-target predictions of computational sequence-based methods by means of anticorrelation analysis of miRNA–mRNA expression levels. Repression by animal miRNAs, differently from plant miRNAs, leads to decreased translational efficiency and/or decreased mRNA levels. Although, the relative contributions of these two outcomes is still unknown and increasing experimental evidences show that changes in mRNA levels closely reflects the impact of miRNAs on gene expression suggesting that destabilization of target mRNAs by exonucleolytic activity is the main mechanism to decrease protein output (Baek et al., 2008; Hendrickson et al., 2009; Guo et al., 2010; Subtelny et al., 2014). Therefore, the anticorrelation analysis of miRNA–mRNA expression levels may contribute to elucidate large portion of miRNA–mRNA regulatory networks affected by pathological conditions. This approach allowed us to identify putative miRNA–target interactions, and cardiac transcription factors that are particularly interesting in the context of HOS.

One of the most interesting miRNA identified as upregulated both at 24 and 48 hpf was miR-34a. MiR-34 family members (miR-34a, -34b, and -34c) are upregulated in the heart in response to stress and the silencing of the entire miR-34 family could protect the heart against pathological cardiac remodeling and improve cardiac functions (Bernardo et al., 2012). Moreover, miR-34a is induced in the aging heart and in vivo silencing of miR-34a reduces age-associated cardiomyocyte cell death. The inhibition of miR-34a reduces cell death and fibrosis following acute myocardial infarction and improves recovery of myocardial function (Boon et al., 2013). These recent studies show an emerging role of miR-34a (and the miR-34 family) as potential regulator of heart remodeling. Therapies that inhibit miR-34a could be useful for cardiac pathologies and HOS.

We discovered several potential miR-34a targets with a possible connection with heart development and HOS (Figure 2). The ATPase Na⁺/K⁺transporting, alpha 2a polypeptide (ATP1a2a) transcript is downregulated in 24 hpf Tbx5a morphants. Since the ATP1a2a contributes to the Ca homeostasis by pumping sodium ions (Na⁺) out of cells and potassium ions (K⁺) into cells, a downregulation of this enzyme might have a negative impact in cardiac contractility and control of arrhythmias. This observation is consistent with the crucial role of Tbx5 in the regulation of cardiac contraction in embryos and in adults. Interestingly, HOS patients show diastolic filling abnormalities and downregulation of ATP2a2, which regulates Ca fluxes in the SR (Zhu et al., 2008).

Furthermore, our data suggest that miR-34a might impact cardiac contraction by regulating a member of the ras homolog gene family, rhoad. RhoA, controlling the Rho-kinase pathway, plays an important role in various fundamental cellular functions, including contraction and motility (Satoh et al., 2011). Moreover, in line with the pro apoptotic role exerted by miR-34a (Raver-Shapira et al., 2007), we observed upregulation of dual-specificity tyrosine-(Y)-phospho-regulated kinase 2 (dyrk2), putative miR-30a target, which negatively regulates the cardiomyocyte growth by mediating repressor function of GSK-3 beta on eIF2B (Weiss et al., 2013) and upregulation of caspase 8, putative target of miR-219-5p downregulated at 48 hpf (Figure 3).

As expected, several genes affecting cardiac morphogenesis were identified in our analysis (Figures 2 and 3). Specifically Roundabout Guidance Receptor 1 (Robo1), which is involved in heart tube formation in zebrafish (Fish et al., 2011) and tmem2, whose expression in myocardial and endocardial tissues in zebrafish and mouse is required for regionally restrict atrioventricular canal boundary and endocardial cushion development. Both genes are putative miR-30a targets at 24 hpf.

Recently, a role for Tbx5 in the establishment of correct heart asymmetry in zebrafish embryos has been highlighted (Pi-Roig et al., 2014). Our data suggest that miR-30a and miR-10d might be contribute to this specific Tbx5 function by controlling respectively bmpr2a (Monteiro et al., 2008) and camk2d1 (Francescatto et al., 2010).

Another interesting miRNA that we found upregulated at 48 hpf is miR-21 whose deregulation in heart has been reported to contribute to cardiovascular disease (Jazbutyte and Thum, 2010). More recently, a crucial role of this miRNA in heart valve formation has been also shown in zebrafish (Banjo et al., 2013), and the alteration of cardiac valve morphology is one of the hallmark of zebrafish HOS phenotype (Camarata et al., 2010; Chiavacci et al., 2012). Two predicted targets of miR-21 are NDRG1B and NDRG4, members of the N-myc downstream regulated gene (NDRG) family, which are downregulated in Tbx5 morphants at 48 hpf. Alterations of NDRG4 cause several of the cardiac defects that characterize the heartstring mutants and are significantly decreased in hearts with reduced Tbx5 activities (Qu et al., 2008). Therefore, we hypothesized that Tbx5 might affect NDRG4 expression through miR-21 modulation.

Although in this study, we used whole embryos for our analysis, we discovered important alterations on transcription factors with crucial roles in heart development. In particular, we observed downregulation both at 24and 48 hpf of the bHLH transcription factor Hand2 (Yelon et al., 2000). Mutations in the hands off locus, which encodes for this transcription factor, cause defects in myocardial development in an early stage, produce a reduced number of myocardial precursors, and show delayed differentiation of the pectoral fin mesenchyme (Schindler et al., 2014). All these phenotypic characteristics are in line with the observed Tbx5 morphant phenotype. In HOS mouse hearts, a strong downregulation of Hand1 was observed (Mori et al., 2006). In mouse, Hand1 and Hand2 are members of the Hand subfamily and have partially redundant functions (Yelon et al., 2000; Tamura et al., 2014). However, in zebrafish, only one member of the hand family has been identified, the Hand2 transcription factor, which is able to perform several of the functions that in mammals are Hand1 specific (Togi et al., 2006; Reichenbach et al., 2008). Therefore, modulation of Hand1 in mouse or Hand2 in zebrafish might have similar functional consequences. Indeed, it has been shown that Hand2 is able to downregulate Nppa, a direct target of Tbx5 and Irx4 an other important cardiac transcription factor strongly downregulated in HOS mouse heart (Bruneau et al., 2001; Mori et al., 2006). In our analysis, we were not able to detect a significant modulation of Nppa gene whose expression is restricted to the cardiac tissue. On the contrary, at 48 hpf, we detected a downregulation of Irx4. This gene is not only expressed in the heart tissues but also present in the eye, a district which is relatively large at this time of development and where Tbx5 is functionally active.

Differently from Hand2, MEF2AA, MEF2CA, SRFB, and SRFA resulted upregulated in 24 hpf embryos depleted for Tbx5a. Among them, MEF2CA is a putative target of miR-10d, both miRNAs already mentioned as downregulated at 24 hpf. All of them codify for transcription factors largely expressed in mesodermal tissues and involved in cardiac developmental patterns highly active at 24 hpf. Consequently, alterations in the expression of these factors have important effects on cardiac development. However, it is hard to predict whether dysregulation of these genes might have positive or negative regulatory effects on their targets. For Tbx5 direct interactors, such as MEF2CA, the ratio between interactors seems more important than the absolute level of the specific factor (Takeuchi et al., 2011). For SRFs, it has been shown that a mild increase may pose either positive and/or negative modulatory effects on their target activation depending on the cofactor recruited (Miano, 2003; Zhang et al., 2003). Interestingly, a negative functional cooperator of SRF is SRFBP1 that we identified as upregulated in 24 hpf Tbx5a morphants. SRFBP1 is highly expressed in fetal and adult mouse heart and functionally interacts with SRF and myocardin in repressing the atrial natriuretic factor promoter activity (Zhang et al., 2004). The data suggest that the observed mild increase of SRF and SRFBP1 in zebrafish Tbx5a morphants might contribute in the downregulation of Nppa, which characterizes the HOS disease.

In conclusion, in this study, we proposed an integrative analysis of miRNA and mRNA expression profiles in a zebrafish model to study the impact of the downregulation of Tbx5 responsible of the HOS. We found several deregulated transcripts including important transcription factors for heart development and diseases, and several deregulated miRNAs with a potential role in the pathology. This model uncovered novel miRNAs and protein coding genes not considered before in the HOS such as miR-34a and miR-30 and their targets. Further dissection of these regulatory circuits will shed light on fundamental pathways in heart development that can contribute to the pathogenesis of human heart diseases. Identification of new TBX5 targets might not only help understand the complexity of HOS phenotype but also contribute in finding novel therapeutic strategies to treat congenital disease. Future experiments are needed to test the role of the identified miRNAs regulated by Tbx5 and the effects on their downstream targets.

4.1. Data Accession Codes

The data discussed in this publication have been deposited in NCBI’s Gene Expression Omnibus and are accessible through GEO Super Series accession number GSE64466.⁴

Author Contributions

RD and LP conceived the work and interpreted the data; RD analyzed sequencing and microarray data; FR performed the integrative analysis; EC performed the in vivo experiments; MB and MG carried out sequencing experiments, MD and IA carried out microarray experiments; RD, FR, EC, and LP wrote the manuscript; GR, LP, and MP supervised the work; all authors read and approved the final manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

FR has been supported by a fellowship sponsored by Progetto Istituto Toscano Tumori Grant 2012 Prot.A00GRT.

Funding

The present work is partially supported by the Flagship project InterOmics (PB.P05, CUP B91J12000270001), funded by the Italian MIUR and CNR organizations, and by the joint IIT-IFC Laboratory of Integrative Systems Medicine (LISM).

Supplementary Material

The Supplementary Material for this article can be found online at https://www.frontiersin.org/article/10.3389/fbioe.2016.00060

Footnotes

1. ^http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
2. ^http://www.bioconductor.com
3. ^http://genome.ucsc.edu/
4. ^http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE64466

References

Albert, M. H., Mannert, J., Fleischmann, K. K., Schiemann, M., Pagel, P., Schmid, I., et al. (2014). Mirnome and transcriptome aided pathway analysis in human regulatory t cells. Genes Immun. 15, 303–312. doi:10.1038/gene.2014.20

PubMed Abstract | CrossRef Full Text | Google Scholar

Anders, S., and Huber, W. (2010). Differential expression analysis for sequence count data. Genome Biol. 11, R106. doi:10.1186/gb-2010-11-10-r106

PubMed Abstract | CrossRef Full Text | Google Scholar

Arnolds, D. E., Liu, F., Fahrenbach, J. P., Kim, G. H., Schillinger, K. J., Smemo, S., et al. (2012). Tbx5 drives scn5a expression to regulate cardiac conduction system function. J. Clin. Invest. 122, 2509–2518. doi:10.1172/JCI62617

PubMed Abstract | CrossRef Full Text | Google Scholar

Asp, M. L., Martindale, J. J., Heinis, F. I., Wang, W., and Metzger, J. M. (2013). Calcium mishandling in diastolic dysfunction: mechanisms and potential therapies. Biochim. Biophys. Acta 1833, 895–900. doi:10.1016/j.bbamcr.2012.09.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Baek, D., Villén, J., Shin, C., Camargo, F. D., Gygi, S. P., and Bartel, D. P. (2008). The impact of microRNAs on protein output. Nature 455, 64–71. doi:10.1038/nature07242

PubMed Abstract | CrossRef Full Text | Google Scholar

Ballman, K. V., Grill, D. E., Oberg, A. L., and Therneau, T. M. (2004). Faster cyclic loess: normalizing RNA arrays via linear models. Bioinformatics 20, 2778–2786. doi:10.1093/bioinformatics/bth327

PubMed Abstract | CrossRef Full Text | Google Scholar

Banjo, T., Grajcarek, J., Yoshino, D., Osada, H., Miyasaka, K. Y., Kida, Y. S., et al. (2013). Haemodynamically dependent valvulogenesis of zebrafish heart is mediated by flow-dependent expression of mir-21. Nat. Commun. 4, 1978. doi:10.1038/ncomms2978

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartel, D. P. (2004). MicroRNAs: genomics, biogenesis, mechanism, and function. Cell 116, 281–297. doi:10.1016/S0092-8674(04)00045-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Basson, C. T., Bachinsky, D. R., Lin, R. C., Levi, T., Elkins, J. A., Soults, J., et al. (1997). Mutations in human cause limb and cardiac malformation in holt-oram syndrome. Nat. Genet. 15, 30–35. doi:10.1038/ng0197-30

CrossRef Full Text | Google Scholar

Baumgart, M., Groth, M., Priebe, S., Appelt, J., Guthke, R., Platzer, M., et al. (2012). Age-dependent regulation of tumor-related microRNAs in the brain of the annual fish Nothobranchius furzeri. Mech. Ageing Dev. 133, 226–233. doi:10.1016/j.mad.2012.03.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Statist. Soc. B 57, 289–300.

Google Scholar

Bernardo, B. C., Gao, X.-M., Winbanks, C. E., Boey, E. J., Tham, Y. K., Kiriazis, H., et al. (2012). Therapeutic inhibition of the mir-34 family attenuates pathological cardiac remodeling and improves heart function. Proc. Natl. Acad. Sci. U.S.A. 109, 17615–17620. doi:10.1073/pnas.1206432109

PubMed Abstract | CrossRef Full Text | Google Scholar

Boon, R. A., Iekushi, K., Lechner, S., Seeger, T., Fischer, A., Heydt, S., et al. (2013). MicroRNA-34a regulates cardiac ageing and function. Nature 495, 107–110. doi:10.1038/nature11919

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruneau, B. G., Nemer, G., Schmitt, J. P., Charron, F., Robitaille, L., Caron, S., et al. (2001). A murine model of holt-oram syndrome defines roles of the t-box transcription factor tbx5 in cardiogenesis and disease. Cell 106, 709–721. doi:10.1016/S0092-8674(01)00493-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Buck, C. A., Baldwin, H., DeLisser, H., Mickanin, C., Shen, H., Kennedy, G., et al. (1993). Cell adhesion receptors and early mammalian heart development: an overview. C. R. Acad. Sci. III 316, 838–859.

PubMed Abstract | Google Scholar

Camarata, T., Krcmery, J., Snyder, D., Park, S., Topczewski, J., and Simon, H.-G. (2010). Pdlim7 (LMP4) regulation of tbx5 specifies zebrafish heart atrio-ventricular boundary and valve formation. Dev. Biol. 337, 233–245. doi:10.1016/j.ydbio.2009.10.039

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, K., and Rajewsky, N. (2007). The evolution of gene regulation by transcription factors and microRNAs. Nat. Rev. Genet. 8, 93–103. doi:10.1038/nrg1990

PubMed Abstract | CrossRef Full Text | Google Scholar

Chiavacci, E., D’Aurizio, R., Guzzolino, E., Russo, F., Baumgart, M., Groth, M., et al. (2015). MicroRNA 19a replacement partially rescues fin and cardiac defects in zebrafish model of holt oram syndrome. Sci. Rep. 5, 18240. doi:10.1038/srep18240

PubMed Abstract | CrossRef Full Text | Google Scholar

Chiavacci, E., Dolfi, L., Verduci, L., Meghini, F., Gestri, G., Evangelista, A. M. M., et al. (2012). MicroRNA 218 mediates the effects of tbx5a over-expression on zebrafish heart development. PLoS ONE 7:e50536. doi:10.1371/journal.pone.0050536

PubMed Abstract | CrossRef Full Text | Google Scholar

Crocini, C., Coppini, R., Ferrantini, C., Yan, P., Loew, L. M., Tesi, C., et al. (2014). Defects in t-tubular electrical activity underlie local alterations of calcium release in heart failure. Proc. Natl. Acad. Sci. U.S.A. 111, 15196–15201. doi:10.1073/pnas.1411557111

PubMed Abstract | CrossRef Full Text | Google Scholar

Da Wei Huang, B. T. S., and Lempicki, R. A. (2008). Systematic and integrative analysis of large gene lists using david bioinformatics resources. Nat. Protoc. 4, 44–57. doi:10.1038/nprot.2008.211

CrossRef Full Text | Google Scholar

Fish, J. E., Wythe, J. D., Xiao, T., Bruneau, B. G., Stainier, D. Y. R., Srivastava, D., et al. (2011). A slit/miR-218/robo regulatory loop is required during heart tube formation in zebrafish. Development 138, 1409–1419. doi:10.1242/dev.060046

PubMed Abstract | CrossRef Full Text | Google Scholar

Francescatto, L., Rothschild, S. C., Myers, A. L., and Tombes, R. M. (2010). The activation of membrane targeted camk-ii in the zebrafish kupffer’s vesicle is required for left-right asymmetry. Development 137, 2753–2762. doi:10.1242/dev.049627

CrossRef Full Text | Google Scholar

Garrity, D. M., Childs, S., and Fishman, M. C. (2002). The heartstrings mutation in zebrafish causes heart/fin Tbx5 deficiency syndrome. Development 129, 4635–4645.

PubMed Abstract | Google Scholar

Gennarino, V. A., Sardiello, M., Avellino, R., Meola, N., Maselli, V., Anand, S., et al. (2009). MicroRNA target prediction by expression analysis of host genes. Genome Res. 19, 481–490. doi:10.1101/gr.084129.108

PubMed Abstract | CrossRef Full Text | Google Scholar

Goetz, S. C., Brown, D. D., and Conlon, F. L. (2006). Tbx5 is required for embryonic cardiac cell cycle progression. Development 133, 2575–2584. doi:10.1242/dev.02420

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, H., Ingolia, N. T., Weissman, J. S., and Bartel, D. P. (2010). Mammalian microRNAs predominantly act to decrease target mRNA levels. Nature 466, 835–840. doi:10.1038/nature09267

PubMed Abstract | CrossRef Full Text | Google Scholar

Hendrickson, D. G., Hogan, D. J., McCullough, H. L., Myers, J. W., Herschlag, D., Ferrell, J. E., et al. (2009). Concordant regulation of translation and mRNA abundance for hundreds of targets of a human microRNA. PLoS Biol. 7:e1000238. doi:10.1371/journal.pbio.1000238

PubMed Abstract | CrossRef Full Text | Google Scholar

Horb, M. E., and Thomsen, G. H. (1999). Tbx5 is essential for heart development. Development 126, 1739–1751.

PubMed Abstract | Google Scholar

Huang, J. C., Babak, T., Corson, T. W., Chua, G., Khan, S., Gallie, B. L., et al. (2007). Using expression profiling data to identify human microRNA targets. Nat. Methods 4, 1045–1049. doi:10.1038/nmeth1130

PubMed Abstract | CrossRef Full Text | Google Scholar

Jazbutyte, V., and Thum, T. (2010). MicroRNA-21: from cancer to cardiovascular disease. Curr. Drug Targets 11, 926–935. doi:10.2174/138945010791591403

PubMed Abstract | CrossRef Full Text | Google Scholar

Kertesz, M., Iovino, N., Unnerstall, U., Gaul, U., and Segal, E. (2007). The role of site accessibility in microRNA target recognition. Nat. Genet. 39, 1278–1284. doi:10.1038/ng2135

PubMed Abstract | CrossRef Full Text | Google Scholar

Kimmel, C. B., Ballard, W. W., Kimmal, S. R., Ullmann, B., and Schilling, T. F. (1995). Stages of embryonic-development of the zebrafish. Dev. Dyn. 203, 253–310. doi:10.1002/aja.1002030302

CrossRef Full Text | Google Scholar

Kwee, L., Baldwin, H. S., Shen, H. M., Stewart, C. L., Buck, C., Buck, C. A., et al. (1995). Defective development of the embryonic and extraembryonic circulatory systems in vascular cell adhesion molecule (VCAM-1) deficient mice. Development 121, 489–503.

PubMed Abstract | Google Scholar

Lewis, B. P., Burge, C. B., and Bartel, D. P. (2005). Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets. Cell 120, 15–20. doi:10.1016/j.cell.2004.12.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, M. (2011). Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.J. 17, 10–12. doi:10.14806/ej.17.1.200

CrossRef Full Text | Google Scholar

McCarthy, D. J., and Smyth, G. K. (2009). Testing significance relative to a fold-change threshold is a treat. Bioinformatics 25, 765–771. doi:10.1093/bioinformatics/btp053

PubMed Abstract | CrossRef Full Text | Google Scholar

McCue, H. V., Haynes, L. P., and Burgoyne, R. D. (2010). Bioinformatic analysis of CaBP/calneuron proteins reveals a family of highly conserved vertebrate Ca2+-binding proteins. BMC Res. Notes 3:118. doi:10.1186/1756-0500-3-118

PubMed Abstract | CrossRef Full Text | Google Scholar

Miano, J. M. (2003). Serum response factor: toggling between disparate programs of gene expression. J. Mol. Cell. Cardiol. 35, 577–593. doi:10.1016/S0022-2828(03)00110-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Monteiro, R., van Dinther, M., Bakkers, J., Wilkinson, R., Patient, R., ten Dijke, P., et al. (2008). Two novel type II receptors mediate bmp signalling and are required to establish left-right asymmetry in zebrafish. Dev. Biol. 315, 55–71. doi:10.1016/j.ydbio.2007.11.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Mori, A. D., Zhu, Y., Vahora, I., Nieman, B., Koshiba-Takeuchi, K., Davidson, L., et al. (2006). Tbx5-dependent rheostatic control of cardiac gene expression and morphogenesis. Dev. Biol. 297, 566–586. doi:10.1016/j.ydbio.2006.05.023

CrossRef Full Text | Google Scholar

Mosimann, C., Panáková, D., Werdich, A. A., Musso, G., Burger, A., Lawson, K. L., et al. (2015). Chamber identity programs drive early functional partitioning of the heart. Nat. Commun. 6, 8146. doi:10.1038/ncomms9146

PubMed Abstract | CrossRef Full Text | Google Scholar

Muniategui, A., Pey, J., Planes, F. J., and Rubio, A. (2013). Joint analysis of miRNA and mRNA expression data. Brief. Bioinformatics 14, 263–278. doi:10.1093/bib/bbs028

PubMed Abstract | CrossRef Full Text | Google Scholar

Nazarov, P. V., Reinsbach, S. E., Muller, A., Nicot, N., Philippidou, D., Vallar, L., et al. (2013). Interplay of microRNAs, transcription factors and target genes: linking dynamic expression changes to function. Nucleic Acids Res. 41, 2817–2831. doi:10.1093/nar/gks1471

PubMed Abstract | CrossRef Full Text | Google Scholar

Peterson, S. M., Thompson, J. A., Ufkin, M. L., Sathyanarayana, P., Liaw, L., and Congdon, C. B. (2014). Common features of microRNA target prediction tools. Front. Genet. 5:23. doi:10.3389/fgene.2014.00023

PubMed Abstract | CrossRef Full Text | Google Scholar

Pi-Roig, A., Martin-Blanco, E., and Minguillon, C. (2014). Distinct tissue-specific requirements for the zebrafish tbx5 genes during heart, retina and pectoral fin development. Open Biol. 4, 140014. doi:10.1098/rsob.140014

PubMed Abstract | CrossRef Full Text | Google Scholar

Qu, X., Jia, H., Garrity, D. M., Tompkins, K., Batts, L., Appel, B., et al. (2008). Ndrg4 is required for normal myocyte proliferation during early cardiac development in zebrafish. Dev. Biol. 317, 486–496. doi:10.1016/j.ydbio.2008.02.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Raver-Shapira, N., Marciano, E., Meiri, E., Spector, Y., Rosenfeld, N., Moskovits, N., et al. (2007). Transcriptional activation of miR-34a contributes to p53-mediated apoptosis. Mol. Cell 26, 731–743. doi:10.1016/j.molcel.2007.05.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Reichenbach, B., Delalande, J.-M., Kolmogorova, E., Prier, A., Nguyen, T., Smith, C. M., et al. (2008). Endoderm-derived sonic hedgehog and mesoderm Hand2 expression are required for enteric nervous system development in zebrafish. Dev. Biol. 318, 52–64. doi:10.1016/j.ydbio.2008.02.061

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. E., Silver, J., Oshlack, A., Holmes, M., Diyagama, D., Holloway, A., et al. (2007). A comparison of background correction methods for two-colour microarrays. Bioinformatics 23, 2700–2707. doi:10.1093/bioinformatics/btm412

PubMed Abstract | CrossRef Full Text | Google Scholar

Satoh, K., Fukumoto, Y., and Shimokawa, H. (2011). Rho-kinase: important new therapeutic target in cardiovascular diseases. Am. J. Physiol. Heart Circ. Physiol. 301, H287–H296. doi:10.1152/ajpheart.00327.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Schindler, Y. L., Garske, K. M., Wang, J., Firulli, B. A., Firulli, A. B., Poss, K. D., et al. (2014). Hand2 elevates cardiomyocyte production during zebrafish heart development and regeneration. Development 141, 3112–3122. doi:10.1242/dev.106336

PubMed Abstract | CrossRef Full Text | Google Scholar

Small, E. M., and Olson, E. N. (2011). Pervasive roles of microRNAs in cardiovascular biology. Nature 469, 336–342. doi:10.1038/nature09783

PubMed Abstract | CrossRef Full Text | Google Scholar

Smyth, G. K. (2004). Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat. Appl. Genet. Mol. Biol. 3, Article3. doi:10.2202/1544-6115.1027

PubMed Abstract | CrossRef Full Text | Google Scholar

Subtelny, A. O., Eichhorn, S. W., Chen, G. R., Sive, H., and Bartel, D. P. (2014). Poly(a)-tail profiling reveals an embryonic switch in translational control. Nature 508, 66–71. doi:10.1038/nature13007

PubMed Abstract | CrossRef Full Text | Google Scholar

Takeuchi, J. K., Lou, X., Alexander, J. M., Sugizaki, H., Delgado-Olguín, P., Holloway, A. K., et al. (2011). Chromatin remodelling complex dosage modulates transcription factor function in heart development. Nat. Commun. 2, 187. doi:10.1038/ncomms1187

PubMed Abstract | CrossRef Full Text | Google Scholar

Tamura, M., Amano, T., and Shiroishi, T. (2014). The Hand2 gene dosage effect in developmental defects and human congenital disorders. Curr. Top. Dev. Biol. 110, 129–152. doi:10.1016/B978-0-12-405943-6.00003-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Togi, K., Yoshida, Y., Matsumae, H., Nakashima, Y., Kita, T., and Tanaka, M. (2006). Essential role of Hand2 in interventricular septum formation and trabeculation during cardiac development. Biochem. Biophys. Res. Commun. 343, 144–151. doi:10.1016/j.bbrc.2006.02.122

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W.-C., Lin, F.-M., Chang, W.-C., Lin, K.-Y., Huang, H.-D., and Lin, N.-S. (2009). miRExpress: analyzing high-throughput sequencing data for profiling microrna expression. BMC Bioinformatics 10:328. doi:10.1186/1471-2105-10-328

CrossRef Full Text | Google Scholar

Weiss, C. S., Ochs, M. M., Hagenmueller, M., Streit, M. R., Malekar, P., Riffel, J. H., et al. (2013). Dyrk2 negatively regulates cardiomyocyte growth by mediating repressor function of GSK-3 on eIF2B. PLoS ONE 8:e70848. doi:10.1371/journal.pone.0070848

CrossRef Full Text | Google Scholar

Westerfield, M. (1993). The Zebrafish Book: A Guide for the Laboratory Use of Zebrafish (Brachydanio rerio). Oregon: University of Oregon Press Eugene.

Google Scholar

Yao, Y., Ma, L., Jia, Q., Deng, W., Liu, Z., Zhang, Y., et al. (2014). Systematic characterization of small RNAome during zebrafish early developmental stages. BMC Genomics 15:117. doi:10.1186/1471-2164-15-117

PubMed Abstract | CrossRef Full Text | Google Scholar

Yates, L. A., Norbury, C. J., and Gilbert, R. J. C. (2013). The long and short of microRNA. Cell 153, 516–519. doi:10.1016/j.cell.2013.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Yelon, D., Ticho, B., Halpern, M. E., Ruvinsky, I., Ho, R. K., Silver, L. M., et al. (2000). The bHLH transcription factor hand2 plays parallel roles in zebrafish heart and pectoral fin development. Development 127, 2573–2582.

PubMed Abstract | Google Scholar

Yue, D., Liu, H., and Huang, Y. (2009). Survey of computational algorithms for microRNA target prediction. Curr. Genomics 10, 478–492. doi:10.2174/138920209789208219

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X., Azhar, G., Furr, M. C., Zhong, Y., and Wei, J. Y. (2003). Model of functional cardiac aging: young adult mice with mild overexpression of serum response factor. Am. J. Physiol. Regul. Integr. Comp. Physiol. 285, R552–R560. doi:10.1152/ajpregu.00631.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X., Azhar, G., Zhong, Y., and Wei, J. Y. (2004). Identification of a novel serum response factor cofactor in cardiac gene regulation. J. Biol. Chem. 279, 55626–55632. doi:10.1074/jbc.M405945200

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, Y., Gramolini, A. O., Walsh, M. A., Zhou, Y.-Q., Slorach, C., Friedberg, M. K., et al. (2008). Tbx5-dependent pathway regulating diastolic function in congenital heart disease. Proc. Natl. Acad. Sci. U.S.A. 105, 5519–5524. doi:10.1073/pnas.0801779105

PubMed Abstract | CrossRef Full Text | Google Scholar