**4. Conclusion**

In the last years, the R&D of SPADs technology has proceeded at great sped. The main challenge towards imaging devices is geometry and process optimization to yield high compactness while not excessively degrading performance.

On the other hand to make a Time Resolved Camera a compact object available on the market, many efforts have still to be made, in terms of identifying the best readout strategies and ancillary circuitries.


Author details

Sara Lehmann-Konera1 , Marek Ruman2,3\*, Krystyna Kozioł 1 , Grzegorz Gajek4 and Żaneta Polkowska1 \*

\*Address all correspondence to<EMAIL_ADDRESS><EMAIL_ADDRESS>
1 Department of Analytical Chemistry, Faculty of Chemistry, Gdansk University of Technology, Gdansk, Poland

2 Faculty of Earth Sciences, University of Silesia, Sosnowiec, Poland

3 Centre for Polar Studies, Leading National Research Centre (KNOW), Sosnowiec, Poland

4 Department of Geology and Lithosphere Protection, Faculty of Earth Sciences and Spatial Management, Maria Curie-Skłodowska University in Lublin, Lublin, Poland


For the evaluation of our method, we applied a 10-fold crossvalidation procedure repeated 30 times to obtain more reliable results. Each fold is built in the following way. Given the tripartite graph, we selected all possible pairs of ncRNA-disease interactions. Then, we randomly partitioned them into each fold. We make sure that the tripartite network generated from each fold is not disconnected. *ncPred* makes predictions only on connected networks. We considered the following four metrics (Alaimo et al., 2013) to assess the performance of our method: precision and recall enhancement, recovery, personalization, and Surprizal. The first two establish the ability of the method to recover the interactions of the test set, therefore, obtaining biologically relevant predictions. The other two measure the ability of the method to propose unexpected interactions, which may lead to novel insights onto ncRNA functions. Special care should be given to the precision and recall enhancement metrics. They measure the reliability of the prediction algorithm by comparing the standard precision and recall with a null model. Such a model is defined as a methodology that randomly assigns ncRNA-disease pairs. This implies that values greater than one are to be considered synonymous of higher quality and, therefore, reliability.

# **3. RESULTS**

As stated earlier, to evaluate the power of our method, we applied a 10-fold cross-validation procedure repeated 30 times and averaged results to obtain more reliable estimates. In **Table 2**, we illustrate the behavior of *ncPred*, comparing it with Yang et al. (2014), in terms of precision and recall enhancement. The results demonstrate that *ncPred* clearly outperforms its competitor. In particular, we can see that while Yang et al. (2014) obtains a recall close to the null model, *ncPred* has much better results. This is crucial since the recall measures the ability of the algorithm to recover existing interactions in the network, and is therefore a sign of their reliability, namely their biological relevance.

In **Figure 3**, we report the receiver operating characteristic (ROC) curves computed on both datasets. The simulations were repeated 30 times and their results were averaged to obtain a more accurate evaluation. Both methods show a high true positive rate against low false positive rate, although *ncPred* is clearly able to achieve better results. This is also shown in **Table 2**, where we can see a significant increase in the average area under the ROC curve (AUC). Such a significance is further proved by the results shown in **Table 3**. By applying the Friedman rank sum test,we determined that the performance improvement achieved by our algorithm is

**Table 2 | Comparison of ncPred andYang et al. (2014) through the precision and recall enhancement metric, and the average area under ROC curve (AUC) calculated for each of the two datasets listed in Table 1**.

The results were obtained using the optimal values for λ1 and λ2 parameters as shown in**Table 3**.

quality of the algorithms in terms of false positives rate against true positives

**Table 4**.

statistically significant (i.e., the *p*-value is close to zero on both datasets).

Regarding the parameters λ1 and λ2, we performed a comprehensive analysis to establish the relationship between them and the prediction quality. In the supporting materials, we report the results of such analysis. The results indicate that there is no specific law, which governs their behavior. The peculiar characteristics of each dataset greatly affect the performances and, consequently, the parameters. It is, therefore, necessary to perform an *a priori* analysis in order to determine, which values give the best results. In our experiments, we used such an analysis to determine the best parameters in terms of precision and recall enhancement (see **Table 4** for details on their values). By looking at the characteristics of our data sets, the values obtained from such an analysis allowed us to suppose that the two parameters are close to zero in Helwak et al. (2013) dataset because of the greater density. This implies that to maintain high quality predictions it is necessary to reduce their number to avoid the introduction of noise. On the other hand, the Chen et al. (2013) dataset has a lower density. This allows us to produce a higher number of predictions before they start losing quality. Therefore, this explains the lambda values closer to one. It is important to point out that in order to determine the best parameters an analysis was performed considering only precision and recall enhancement, since they are closely related to the biological significance of the predictions. In this context, we report in **Table 2** only precision and recall enhancement and the AUC, ignoring the other metrics, which are available in the supporting materials.

Finally, assuming that the number of targets dominates the ncRNA one, we can state that the computational complexity of our method is *O*(*m*2 *p*). However, it is quite straightforward to implement parallelization and optimization techniques to make the computation faster.

#### **3.1. CASE STUDIES**

The analysis of the predictions for each non-coding showed that *ncPred* is able to find exactly the same predictions provided byYang et al. (2014). The main difference between the two algorithms lies in the different scores given to each prediction. As highlighted in the previous section, ncPred is clearly able to provide more substantially accurate predictions.

# **Table 3 | Friedman rank sum test applied to establish the statistical significance in the performance improvement of ncPred compared to Yang et al. (2014)**.

## **Table 4 | Optimal values of** λ**1 and** λ**2 parameters for the datasets used in our experiments**.

To further demonstrate the ability of our method, we reviewed in detail the results of five diseases (i.e., Alzheimer's Disease, Myocardial Infarction, Pancreatic Cancer, Parkinson's disease, and Gastric Cancer) as case studies. The top 10 predicted genes for each case are listed in **Table 5**. **Table 5** also shows the rank obtained by applying on our dataset, theYang et al. (2014) method. In this context, the two datasets were taken together in order to start from a wider knowledge base.

#### **3.1.1. Alzheimer's disease**

Alzheimer's disease (AD) is one of the most common forms of dementia (Hebert et al., 2003). Recent studies indicate that it affects approximately 0.40% of the world population (Brookmeyer et al., 2007). The disease is, at present, untreatable, and it is characterized by a progressive loss of mnemonic, cognitive, and intellectual capacity, which ultimately leads to the death of the patient. Among the first 10 ncRNAs, we find *PVT1* a lncRNA,

**Table 5 | List of top 10 predictions computed by ncPred and their rank obtained withYang et al. (2014) for five case studies (Alzheimer's Disease, Myocardial Infarction, Pancreatic Cancer, Parkinson's Disease, and Gastric Cancer)**.

which regulates the transcription of *MYC* on the long distance (Carramusa et al., 2007). In Jiang et al. (2013),*MYC* has been characterized as the source of the main pathway substantially active in AD, thus having an important role in disease progression. Such a discovery confirms that *PVT1* could play a key role in the progress of AD. We have also identified the lncRNA *MEG3* that activates *TP53* and improves its binding affinity to target gene promoter (Liao et al., 2011). *TP53* was identified in Tan et al. (2012) as potential biomarker for AD. Therefore, further analysis to confirm *MEG3* role in AD are needed.

# **3.1.2. Myocardial Infarction**

Myocardial infarction (MI) is a heart condition that occurs when the proper flow of blood to a part of the heart stops, and the heart muscle is damaged due to lack of sufficient oxygen. Genome-wide association studies have identified 27 epigenetic factors that are associated with an increased risk of MI (Feero et al., 2011). For example, the genomic locus 9p21 has one of the strongest associations with the pathology (Feero et al., 2011). The majority of such factors have been identified in regions implicated in other heart diseases (Feero et al., 2011). Among our predictions, we identified the lncRNA *SRA1* that Friedrichs et al. (2009) found crucial in cardiomyopathies. This leads us to assume a possible link with MI. In the top 10 predictions we also found the lncRNA *7SL*, which, by hybridizing to the *reverse-Alu-element-containing 3* 0*UTR of MnSOD* gene, represses its expression (Lipovich et al., 2010). Overexpression of *MnSOD* has been identified as a possible protection against MI in transgenic mice (Chen et al., 1998). This could be a cue for further investigations to understand the role such a lncRNA.

# **3.1.3.  The analysis of the predictions obtained by our algorithm has provided the association with lncRNA *HOTAIR*, whose overexpression has been associated with a poor prognosis in pancreatic cancer, as well as show a pro-oncogenic activity (Kim et al., 2012). A further lncRNA is *Airn*. The deletion of its promoter in paternal allele results in aberrant activation of *IGF2R* (Nagano and Fraser, 2009), whose polymorphisms are associated with an increased risk of pancreatic cancer (Dong et al., 2012).

# **3.1.4. Parkinson's disease**

Parkinson's disease (PD) is a degenerative disorder of the central nervous system. The main cause of the disease is the death of dopamine-generating cells in the substantia nigra. The cause of this death is still unknown, nevertheless, the process of aging and metabolic stress are its common triggers (Parlato and Liss, 2014). It is interesting to note that the response to stress conditions and mechanisms for quality control are compromised in patients with PD. The reduction in the transcription of rRNA (ribosomal ribonucleic acid) is an important strategy to maintain cellular homeostasis under stress. An altered transcription is associated with neurodegenerative disorders. There are many triggers for nucleolar stress, but they seem to depend on the extitp53 protein (Parlato and Liss, 2014). Our algorithm is able to identify two probable lncRNA associated with this function: *PVT1*, also

associated with AD, whose gene locus is a target of *p53* (Barsotti et al., 2012), and *MEG3* that promotes the expression of *Tp53* and increases the binding affinity to the promoters of its target (Liao et al., 2011).

# **3.1.5. Gastric cancer**

Gastric cancer is a disease typically characterized by an overall 5 years survival rate lower than 10%, mainly due to the plurality of common symptoms that lead to treatments only in advanced disease stages (Orditura et al., 2014). Among our predictions, we find the lncRNA *Xist*. In Weakley et al. (2011), it was identified as differentially expressed in stomach preneoplastic cells, which could be a symptom of gastric cancer. Another factor could be the lncRNA *Evf2*, which is a direct putative positive regulator of transcription factor *Dlx-2* (Lipovich et al., 2010). Increased expression of *Dlx-2* was correlated with more advanced stages of the disease (Tang et al., 2013).

# **4. DISCUSSION**

In this paper, we propose *ncPred* to predict novel associations between ncRNAs and diseases. The aim is to compute ncRNAdisease association's prediction starting from a tripartite network. Such a network integrates information on ncRNAs, targeting (i.e., those genes, microRNAs, proteins whose activity is affected by non-coding RNA), and their associations with diseases in order to improve prediction quality and accuracy.

Our experimental analysis shows that our approach predicts more biologically significant associations with respect toYang et al. (2014). This assertion is confirmed by the results obtained in terms of recall, which as described above measures biological quality of results. The use of Friedman rank sum test also showed that the difference between our predictions and those of Yang et al. (2014) is not random but due to a better interpretation of available information. The results showed that our method could provide interesting suggestions in the study of the implications between ncRNA and pathologies. However, as stated in the introduction, the method can only help to make such a search more targeted and less expensive, offering a ranking of associations from more probable to less probable. Determine whether those associations are useful still remains within the competence area of bio-physicians that can provide conclusive evidence by identifying suitable patients and documenting such cases.

Despite what stated earlier, our method still has some limitations that should be taken into account. Firstly, ncRNA-target associations are still too small in number. It may be necessary to resort to additional targeting prediction techniques so as to expand original knowledge base. Secondly, the methodology does not use the biological information accompanying each association (e.g., type of ncRNA-target interaction, conditions in which the targetdisease association was detected, tissues in which associations have significance). For this reason, it may be useful to further expand the methodology by using such additional information, which could make the methodology more reliable in terms of significant predictions.

# **SUPPLEMENTARY MATERIAL**

In the Supplementary Material (Data Sheet 1.pdf) we report the ncPred parameter tuning further details concerning the comparison with Yang et al. (2014). The Supplementary Material for this article can be found online at http://www.frontiersin.org/Journal/ 10.3389/fbioe.2014.00071/abstract

# **REFERENCES**

**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.

*Received: 29 September 2014; accepted: 25 November 2014; published online: 12 December 2014.*

*Citation: Alaimo S, Giugno R and Pulvirenti A (2014) ncPred: ncRNA-disease association prediction through tripartite network-based inference. Front. Bioeng. Biotechnol. 2:71. doi: 10.3389/fbioe.2014.00071*

*This article was submitted to Bioinformatics and Computational Biology, a section of the journal Frontiers in Bioengineering and Biotechnology.*

*Copyright © 2014 Alaimo, Giugno and Pulvirenti. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited andthatthe original publication inthis journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Discovery of protein–lncRNA interactions by integrating large-scale CLIP-Seq and RNA-Seq datasets

#### **Jun-Hao Li † , Shun Liu† , Ling-Ling Zheng† , JieWu,Wen-Ju Sun, Ze-LinWang, Hui Zhou, Liang-Hu Qu\* and Jian-HuaYang\***

RNA Information Center, Key Laboratory of Gene Engineering of the Ministry of Education, State Key Laboratory for Biocontrol, Sun Yat-sen University, Guangzhou, China

## **Edited by:**

Alessandro Laganà, The Ohio State University, USA

#### **Reviewed by:**

Igor B. Rogozin, National Institutes of Health, USA Thiruvarangan Ramaraj, National Center for Genome Resources, USA

#### **\*Correspondence:**

Liang-Hu Qu and Jian-Hua Yang, Biotechnology Research Center, Sun Yat-sen University, Guangzhou 510275, China e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
†Jun-Hao Li, Shun Liu and Ling-Ling Zheng have contributed equally to this work.

Long non-coding RNAs (lncRNAs) are emerging as important regulatory molecules in developmental, physiological, and pathological processes. However, the precise mechanism and functions of most of lncRNAs remain largely unknown. Recent advances in high-throughput sequencing of immunoprecipitated RNAs after cross-linking (CLIP-Seq) provide powerful ways to identify biologically relevant protein–lncRNA interactions. In this study, by analyzing millions of RNA-binding protein (RBP) binding sites from 117 CLIP-Seq datasets generated by 50 independent studies, we identified 22,735 RBP–lncRNA regulatory relationships. We found that one single lncRNA will generally be bound and regulated by one or multiple RBPs, the combination of which may coordinately regulate gene expression. We also revealed the expression correlation of these interaction networks by mining expression profiles of over 6000 normal and tumor samples from 14 cancer types. Our combined analysis of CLIP-Seq data and genome-wide association studies data discovered hundreds of disease-related single nucleotide polymorphisms resided in the RBP binding sites of lncRNAs. Finally, we developed interactive web implementations to provide visualization, analysis, and downloading of the aforementioned large-scale datasets. Our study represented an important step in identification and analysis of RBP–lncRNA interactions and showed that these interactions may play crucial roles in cancer and genetic diseases.

**Keywords: long non-coding RNA, RNA-binding protein, GWAS, CLIP-Seq, RNA-Seq**

# **INTRODUCTION**

Mammalian genomes encode thousands of long non-coding RNAs (lncRNAs) (Wang and Chang, 2011; Guttman and Rinn, 2012). lncRNAs play important roles in a variety of biological processes that have been implicated in regulating tumorigenesis through interaction with RNA-binding proteins (RBPs) (Konig et al., 2011; Wang and Chang, 2011; Guttman and Rinn, 2012; Ulitsky and Bartel, 2013). However, for the majority of lncRNAs, the mechanism underlying their interaction with RBPs remains unknown (Konig et al., 2011; Wang and Chang, 2011; Guttman and Rinn, 2012; Ulitsky and Bartel, 2013).

The control and function of lncRNA are governed by the specificity of RBPs (Wang and Chang, 2011; Guttman and Rinn, 2012). Increasing evidence suggests that many RBP–lncRNA interactions play important roles in correct transcriptional regulation (Konig et al., 2011; Wang and Chang, 2011; Guttman and Rinn, 2012; Ulitsky and Bartel, 2013). One emerging theme that many lncR-NAs regulate gene expression by directing chromatin modificators to specific target regions (Ulitsky and Bartel, 2013). Significant fractions (20% in human) of lincRNAs are interacted with PRC2 and other chromatin-modifying complexes (Khalil et al., 2009; Guttman et al., 2011). The functional outcomes of some binding events have been revealed. For example, HOTAIR, which is transcribed from human HOX locus, guides repressor PRC2 to specific mammalian loci to silence gene expression and to promote cancer metastasis (Rinn et al., 2007; Wang et al., 2011). Besides, many lncRNAs have been shown to interact with other types of RBPs, including DNA methyltransferases (Schmitz et al., 2010; Di Ruscio et al., 2013), transcription factors (Wang et al., 2014), and splicing factors (Tripathi et al., 2010; Gong and Maquat, 2011; Yin et al., 2012). However, deciphering the interactions between hundreds of RBPs and thousands of lncRNAs remains a daunting challenge.

Genome-wide association studies (GWAS) have identified thousands of common genetic variants related to specific traits or disease phenotypes, and many of these variants (about 88%) lie in non-coding regions, which could potentially influence processing and expression of ncRNAs (Sethupathy and Collins, 2008; Hindorff et al., 2009; Ryan et al., 2010; Cabili et al., 2011; Kumar et al., 2013; Ning et al., 2014). For example, single nucleotide polymorphism (SNP) within miR-125a gene alters the processing of pri-miRNA by DGCR8 and causes recurrent pregnancy loss in a Han-Chinese population (Duan et al., 2007; Hu et al., 2011). Another study found that a papillary thyroid carcinoma-associated SNP, rs944289 affects the expression of lncRNA PTCSC3 by changing the binding activity of C/EBPα transcription factor (Cabili et al., 2011; Jendrzejewski et al., 2012). Although the genetic variants in interaction sites of RBP–lncRNA may interfere lncRNA functions and affected the susceptibility to human diseases, the relationships between genetic variants and interaction sites were yet unexplored.

Recent advances in high-throughput sequencing of RNA isolated by cross-linking immunoprecipitation (HITS-CLIP, CLIP-Seq, PAR-CLIP, CLASH, iCLIP) have provided powerful ways to identify RBP-associated RNAs and map such interactions in the genome (Chi et al., 2009;Hafner et al., 2010;Konig et al., 2011;Helwak et al., 2013; Fu, 2014; Fu and Ares, 2014). The application of CLIP-Seq methods has reliably identified Argonaute (Ago) binding sites and miRNA-target interactome (Chi et al., 2009; Hafner et al., 2010; Helwak et al., 2013). In fact, many more studies to date have been focused on understanding the function of RBPs in RNA metabolism (Konig et al., 2011; Fu, 2014), such as pre-mRNA splicing (Fu and Ares, 2014). While an increasing number of RBPs have been explored using CLIP technologies, binding peaks mapped to non-protein-coding genes have been routinely discarded and not further analyzed. However, this data will be a rich trove well worthy of mining RBP–lncRNA relationships.

In this study, we performed a large-scale integration of public RBP binding sites generated by high-throughput CLIP-Seq technology and identified thousands of RBP–lncRNA interactions. Furthermore, by combining GWAS and RNA-Seq data, we explored clinically relevant RPB–lncRNA interactions that may facilitate the translation of genetic studies of complex diseases into therapeutics.

# **MATERIALS AND METHODS**

# **INTEGRATION OF RBP BINDING SITES FROM PUBLISHED CLIP DATA**

HITS-CLIP, PAR-CLIP, and iCLIP binding clusters/peaks data were retrieved from the gene expression omnibus and sequence read archive (SRA) (Barrett et al., 2013), the supplementary data of original references or directly from authors upon request. All binding sites coordinates were converted to hg19 and mm10 assemblies using the UCSC LiftOver Tool (Meyer et al., 2013).

# **RBP TARGET SITES SCANNING IN ANNOTATED lncRNA TRANSCRIPTS**

Human gene annotations were acquired from GENCODE Version 17 (Harrow et al., 2012). Mouse gene annotations were extracted from Ensembl Gene Release 72 (Hubbard et al., 2009) and LiftOver to mm10 assembly. lncRNAs were further filtered to remove the transcripts overlapping with protein-coding genes. The aforementioned RBP CLIP clusters were used to intersect with the coordinates of all annotated transcripts to find their RBP binding sites, which were fed to Circos (Krzywinski et al., 2009) for visualization.

# **TCGA TUMOR EXPRESSION DATA AND EXPRESSION CORRELATION OF RBPs AND lncRNAs**

The Cancer Genome Atlas (TCGA) RNA-Seq expression datasets (level 3, IlluminaHiSeq\_RNASeqV2) for 14 cancer types and gene annotation file (TCGA.hg19.June2011.gaf) were downloaded from TCGA Data Portal (Cancer Genome Atlas Research Network, 2008). Expression of 397 known lncRNAs can be measured in TCGA level 3 RNA-Seq data. Expression correlation (Pearson correlation coefficient) between lncRNAs and RBPs was estimated using co-expression program (the program is available from the authors upon request), which was written in C language and ALGLIB library, and *p*-value was adjusted with the false discovery rate (FDR) correction (Benjamini and Hochberg, 1995).

# **IDENTIFICATION OF DISEASE-RELATED SNPs IN RBP BINDING SITES ASSOCIATED WITH lncRNAs**

Disease/phenotype associated SNPs were curated from published GWAS data provided by the NHGRI GWAS Catalog (Welter et al., 2014), Johnson and O'Donnell (2009), dbGAP (Mailman et al., 2007), and GAD (Becker et al., 2004). Additional SNPs in linkage disequilibrium (LD) with reported disease-related loci were selected with the criteria requiring an *r* 2 value over 0.5 in at least one of the four populations (CEU, CHB, JPT, and YRI) genotype data of the HapMap project (release 28) (International HapMap 3 Consortium et al., 2010). For each SNP, rs ID were lifted to dbSNP build 141 based on the "RsMergeArch.bcp" and "SNPHistory.bcp" table from dbSNP, and genomic coordinates were lifted to the hg19 assembly using the UCSC LiftOver tool. All these disease-related SNPs or LD SNPs were mapped to exons and splicing sites (2 nt in the intron that is close to an exon) of the annotated lncRNA transcripts and further examined whether they were located in any RBP binding clusters.

# **DATA VISUALIZATIONS**

RNA-binding protein–lncRNA interactions were deposited in our starBase V2.0 (Li et al., 2014) under the "Protein–RNA" section1 . For each interaction, we provided links to our enhanced deepView genome browser2 , which was written using a GD graphics library for PHP, to visualize RBP binding sites, lncRNAs, and other annotation tracks in an integrated display style similar to that of UCSC genome browser.

# **RESULTS**

# **THE GENOME-WIDE BINDING MAP OF RNA-BINDING PROTEINS AND THE ANNOTATION OF RBP–lncRNA INTERACTIONS**

We curated 117 published CLIP-Seq datasets to profile the genome-wide binding maps of 65 RBP. Unique binding sites of distinct RBPs varied from thousands to millions, and the genomic context distributions of binding sites for different RBPs distinguished from each other (**Figure 1**; Table S1 in Supplementary Material). For example, PUM2, a translational repressor during embryonic development and cell differentiation (Huang et al., 2011), predominately bound to 30UTR regions of protein genes, while another translation inhibitor FMRP (Napoli et al., 2008) tended to interact with CDS. The discrepancy in binding context preferences for RBPs could root from different amounts of available datasets, usages of various variants of CLIP-Seq, varying sequencing depth, and/or genuine distinctions in the underlying recognition mechanism of RBPs.

Despite that the majority of RBP binding sites were mapped to protein-coding genes, on average 1.1% of RBP binding sites lay within exons of human lncRNAs. In total, 21,073 and 1,662 RBP–lncRNA interactions were identified in human and mouse, respectively (**Table 1**). It is noteworthy that most well-studied lncRNAs interacted with chromatin modificators, acting as tethers or scaffolds (Khalil et al., 2009;Kung et al., 2013). Thus, we considered the binding features of Ezh2, a subunit of PRC2 complexes, by analyzing the CLIP-Seq data from mouse embryonic stem cells

1http://starbase.sysu.edu.cn/rbpLncRNA.php

2http://starbase.sysu.edu.cn/browser.php

# **Table 1 | The summary of CLIP-Seq datasets used in this study and the resulting RBP–lncRNA interactions**.

(Kaneko et al., 2013). Our results demonstrated that Ezh2 interacted with 35 lncRNAs including many imprinted RNAs, such as Tsix, Meg3, Rian, and Pvt1 (**Figure 2**), which was consistent with the epigenetic features of PRC2 (Zhao et al., 2010).

# **EXPLORING COMBINATORIAL EFFECTS AMONG RBPs**

For the 12,255 human lncRNAs, 56.8% were found bound to at least 1 RBP. Surprisingly, 16 lncRNAs, including GAS5 and NEAT1, harbored binding sites of over 30 RBPs (**Figure 3**; Table S2 in Supplementary Material), indicating their diverse roles in biological processes when accompanied with different RBPs. Since one lncRNA could interact with multiple RBPs, it could be expected that some RBP binding sites were overlapped with each other. Therefore, we explored combinatorial effects among RBPs by employing integrated CLIP-Seq datasets. For example, we utilized PAR-CLIP data generated in HEK293 and intersect binding sites of three RNA destabilizer HuR, Ago2, and MOV10. The results showed that tens of lncRNAs, including cancer-related lncRNAs TUG1, DLEU2, and GAS5, were bound by at least two of the

three RBPs at identical binding sites (**Figure 4**). This phenomenon suggested that the stabilities of these lncRNAs were likely under joint control of these three RBPs, which could be explained by their confirmed interplays in HEK293 (Chendrimada et al., 2007) and Hela cells (Kim et al., 2009).

# **EXPRESSION ASSOCIATION OF RBP–lncRNA INTERACTIONS**

To realize the roles of RBP–lncRNA interactions in cancer, we preformed co-expression analysis between RBPs and lncRNAs by virtue of 90 human CLIP-Seq datasets and expression data from more than 6,000 tumor samples in 14 types of cancer. Up to 583 pairs concerning 47 RBPs and 49 lncR-NAs showed strong correlation at expression levels in at least 1 cancer type (**Figure 5A**). Marvelously, PUM2 and TUG1 involved with cell cycle regulation (Khalil et al., 2009; Huang et al., 2011) showed significant positive expression correlation (*p* < 0.05) in all 14 cancer types (**Figure 5B**). Two potential PUM2 binding sites on TUG1 have the consensus recognition

motif UGURUAUA, which was highly conserved in mammals (**Figure 5C**).

# **PREDICTING GWAS-ASSOCIATED RBP BINDING SITES IN lncRNAs**

Although GWAS over the years have revealed a significant number of genetic variants related to diseases or phenotypes, a considerable portion of these identified loci are not within protein-coding genes and therefore not functionally explained to date (Hindorff et al., 2009). Here, we tried to fill this gap by connecting RBP binding sites in lncRNAs and potential disease-related SNPs.

Altogether, 87,677 unique disease-related SNPs were collected from four public GWAS data source (Table S3 in Supplementary Material, detailed in Section "Materials and Methods"). Considering that additional SNPs in LD with reported disease-related loci may also map to RBP binding sites in lncRNAs, we perform LD analysis to extracted SNPs that had high LD relationship with disease-related SNPs using a threshold of *r* 2 > 0.5 in at least one population from the HapMap CEU, CHB, JPT, and YRI genotype data, which yielded a total of 895,968 disease-related or LD SNPs.

We found that 2431 of these SNPs were mapped to the exons of 2089 transcripts of 1489 lncRNA genes, among which 162 SNPs were also located in at least 1 binding sites of 29 RBPs (Table S4 in Supplementary Material). For example, three disease-related SNPs, namely, rs16902485, rs10283090, and rs2720659, resided in the exons of lncRNA PVT1. According to the GWAS annotations of Johnson and O'Donnell (2009), the latter two of the three SNPs were associated with "type II diabetes mellitus," which was in good accordance with the recent reports showing that PVT1 may contribute to diabetic nephropathy (Hanson et al., 2007; Alvarez and DiStefano, 2011; Alwohhaib et al., 2014). These SNPs were also overlapped with binding sites of U2AF65, HuR, and eIF4AIII, respectively (**Figure 6**), suggesting variants in these sites might result in impaired binding of these RBPs to PVT1, which thereby might lead to the development of corresponding diseases.

Next, we checked whether disease-related SNPs might be located in the splicing sites of lncRNAs and affect the alternative splicing of lncRNAs. We defined a splicing site as the 2 nt within an intron close to the exon–intron junction. As a result, we found that only 24 SNPs lay within lncRNA splicing sites (Table S4 in Supplementary Material), among which only 1 SNP, rs17207481, was overlapped with binding sites of FUS and HuR. These results suggested that SNPs exerted limited effects on disease occurrence through the mechanism of disturbing alternative splicing of lncRNAs.

# **COMPARATIVE ANALYSIS OF RBP TARGETS USING THE deepView GENOME BROWSER**

To facilitate comparative analysis of the CLIP-Seq datasets and exploration of RBP–lncRNA interactions, we developed the improved deepView Genome Browser2 in starBase V2.0 (Li et al., 2014). In the query page of the browser, users can input one interested genomic region in the "search term" and select corresponding genome assembly to gain an integrated view of various genomic features. Information on binding sites of RBPs, predicted miRNA-target sites overlapped with CLIP-Seq data, as well as gene annotations from RefSeq and Ensembl were provided in toggleable tracks. The image of the browser will be updated immediately by clicking the "refresh tracks" button when users change track options. **Figure 7** illustrated the visualization of FUS–MEG3 interactions with deepView. Users can click the"zoom in"or"zoom out"

button at the top to shrink or extend on the center of the annotation tracks window by 1.5-, 3-, or 10-folds. Clicking the "View region at UCSC" button will redirect users to the UCSC page and exhibit the current region on the UCSC genome browser. To explore RBP binding sites on a particular gene, users can type its gene symbol in the position textbox and then click the "GO" button to update the display image to determine, which RBPs might participate in regulating the gene. Our visualization method allows a direct comparison of binding patterns of different RBPs, binding preferences of one particular RBP in different cell lines and tissues, and genomic contexts of RBP binding sites.

# **DISCUSSION**

Although a few dozen lncRNAs have been characterized to some extent and reported to function in important cellular processes, the functions of most annotated lncRNAs are

unknown (Guttman and Rinn, 2012; Ulitsky and Bartel, 2013). Several bioinformatics resources and tools have made efforts to functionally annotate lncRNAs (Da Sacco et al., 2012), such as fRNAdb (Kin et al., 2007) and ncFANs (Liao et al., 2011). These tools mainly inferred lncRNA function by their differential expression in distinct biological states or their co-expression patterns with protein-coding genes, but little attention was paid to the relationship of lncRNAs and their bounded proteins. In this study, by analyzing a large set of RBP binding sites derived from all available CLIP-Seq experimental techniques (PAR-CLIP, HITS-CLIP, iCLIP, CLASH), we have shown extensive and complex RBP–lncRNA interaction networks (**Figure 1**).

Recent studies have revealed that many lncRNAs function through specific interactions with RBPs, but whether these interactions are direct and specific remains controversial. RBP–lncRNA interactions identified by low stringent immunoprecipitation of non-cross-linked RNA–protein complexes, such as RIP-Chip and RIP-Seq, may contain indirect binding relationships (Konig et al., 2011). In comparison to previously reported significant fractions (10% in mouse) of PRC2-associated lncRNAs (Zhao et al., 2008), we found that a relatively small fraction (~1%) of lncR-NAs were bound by Ezh2 in mouse (**Figure 2**). Therefore, we provide enhanced resolution to determine lncRNA functional networks based on RBP–lncRNA interactions supported by highthroughput CLIP-Seq data. More than 80,000 binding clusters identified from 65 different RBPs represent a valuable resource

for resolving some obstacles that have arisen in efforts to understand lncRNA action. Nevertheless, although CLIP-Seq is designed to detect direct binding events of proteins and RNAs, the resulting data might still contain false positives and false negatives, which may root from every cumbersome step of this technique. To minimize the impact of such false discoveries, we filtered the origin results by the reported FDR and provided evidences such as number of CLIP reads and number of supporting experiments, which may help users to gain RBP–lncRNA interactions of high-confidence.

By cross analysis of binding maps for multiple RBPs, this study offers a new resource to understanding joint control of target lncRNA expression. While only 65 RBPs were analyzed, we found that many of the RBPs bound to the same lncRNA (**Figure 3**). This is consistent with the compelling idea that lncRNAs can serve as scaffolds that assemble many relevant RBPs to regulate gene expression (Wang and Chang, 2011; Ulitsky and Bartel, 2013). At the same time, we also identified hundreds of identical binding sites that bound by multiple different RBPs in lncRNAs (**Figure 4**), probably reflecting competition among RBPs that binding on a given lncRNA.

Our combined analysis of CLIP-Seq data and GWAS data revealed hundreds of disease-related SNPs resided in the RBP binding sites of lncRNAs (Table S4 in Supplementary Material). Unlike the sporadic attempts on simply finding genetic variants associated with disease susceptibility within lncRNA genes

(Bochenek et al., 2013; Mirza et al., 2014), our approaches focused on SNPs that might impact on the binding events between RBPs and lncRNAs. Since most lncRNAs fulfill their roles through by forming complex with their protein partners, our results provide insights on the functions of lncRNAs from the perspective of RBP binding malfunction in diseases, which in turn may contribute to disease etiology.

Overall, our studies and the accompanying datasets demonstrated that one single lncRNA will generally be bound and regulated by one or multiple RBPs, the combination of which may coordinately determine the final regulatory outcome. We have also shown that an exhaustive and high-resolution RBP–lncRNA interaction map will help to discover genetic variations that contribute to complex genetic diseases by affecting post-transcriptional gene regulation.

# **AUTHOR CONTRIBUTIONS**

Jian-Hua Yang, Liang-Hu Qu, and Jun-Hao Li conceived the project. Jun-Hao Li, Shun Liu, Ling-Ling Zheng, and Jian-Hua Yang performed the computational and statistical analysis. Jun-Hao Li, Shun Liu, Ling-Ling Zheng, Jian-Hua Yang, Liang-Hu Qu, Jie Wu, Wen-Ju Sun, Ze-Lin Wang, and Hui Zhou wrote the manuscript. Liang-Hu Qu and Jian-Hua Yang supervised the project. All authors read and approved the final manuscript.

on MEG3 are visible in the RBP binding sites track. In this track, the

multiple RBPs.

# **ACKNOWLEDGMENTS**

This research is supported by the Ministry of Science and Technology of China, National Basic Research Program (No. 2011CB811300); the National Natural Science Foundation of China (No. 91440110, 31230042, 31370791, 31471223, 31401975); the funds from Guangdong Province (No. S2012010010510, S2013010012457); the project of Science and Technology New

Star in ZhuJiang Guangzhou city (No. 2012J2200025); Fundamental Research Funds for the Central Universities (No.<PHONE_NUMBER>161070, 14lgjc18); China Postdoctoral Science Foundation (No. 200902348). This research is supported in part by the Guangdong Province Key Laboratory of Computational Science and the Guangdong Province Computational Science Innovative Research Team.

# **SUPPLEMENTARY MATERIAL**

The Supplementary Material for this article can be found online at http://www.frontiersin.org/Journal/10.3389/fbioe.2014.00088/ abstract

**Table S1 | The summary of binding sites distribution across genomic features for 47 human RBPs**.

**Table S2 | Counts of binding RBPs in the 12,255 human lncRNAs**.

**Table S3 | Disease/trait related SNPs collected from four public GWAS databases**.

**Table S4 | Disease/trait related SNPs overlapped with RBP binding sites in lncRNAs**.

# **REFERENCES**

**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.

*Received: 19 October 2014; accepted: 22 December 2014; published online: 14 January 2015.*

*Citation: Li J-H, Liu S, Zheng L-L, Wu J, Sun W-J, Wang Z-L, Zhou H, Qu L-H and Yang J-H (2015) Discovery of protein–lncRNA interactions by integrating large-scale CLIP-Seq and RNA-Seq datasets. Front. Bioeng. Biotechnol. 2:88. doi: 10.3389/fbioe.2014.00088*

*This article was submitted to Bioinformatics and Computational Biology, a section of the journal Frontiers in Bioengineering and Biotechnology.*

*Copyright © 2015 Li, Liu, Zheng , Wu, Sun, Wang , Zhou, Qu and Yang . This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

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

*Romina D'Aurizio1 , Francesco Russo1,2, Elena Chiavacci3 , Mario Baumgart4 , Marco Groth4 , Mara D'Onofrio5 , Ivan Arisi5 , Giuseppe Rainaldi1 , Letizia Pitto3 \* and Marco Pellegrini1 \**

#### *Edited by:*

*Christian M. Zmasek, J. Craig Venter Institute, USA*

# *Reviewed by:*

*Eirini Trompouki, Max Planck Institute for Immunobiology and Epigenetics, Germany Claus Jürgen Scholz, University of Würzburg, Germany*

#### *\*Correspondence:*

*Letizia Pitto<EMAIL_ADDRESS>Marco Pellegrini<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Bioinformatics and Computational Biology, a section of the journal Frontiers in Bioengineering and Biotechnology*

*Received: 18 January 2016 Accepted: 24 June 2016 Published: 14 July 2016*

#### *Citation:*

*D'Aurizio R, Russo F, Chiavacci E, Baumgart M, Groth M, D'Onofrio M, Arisi I, Rainaldi G, Pitto L and Pellegrini M (2016) Discovering miRNA Regulatory Networks in Holt–Oram Syndrome Using a Zebrafish Model. Front. Bioeng. Biotechnol. 4:60. doi: 10.3389/fbioe.2016.00060*

*1 Laboratory of Integrative Systems Medicine (LISM), Institute of Informatics and Telematics (IIT), Institute of Clinical Physiology (IFC), National Research Council (CNR), Pisa, Italy, 2Department of Computer Science, University of Pisa, Pisa, Italy, 3 Institute of Clinical Physiology (IFC), National Research Council (CNR), Pisa, Italy, 4 Leibniz Institute on Ageing, Fritz Lipmann Institute (FLI), Jena, Germany, 5Genomics Facility, Fondazione EBRI Rita Levi-Montalcini, Roma, Italy*

MicroRNAs (miRNAs) are small non-coding RNAs that play an important role in the post-transcriptional regulation of gene expression. miRNAs are involved in the regulation of many biological processes such as differentiation, apoptosis, and cell proliferation. miRNAs are expressed in embryonic, postnatal, and adult hearts, and they have a key role in the regulation of gene expression during cardiovascular development and disease. Aberrant expression of miRNAs is associated with abnormal cardiac cell differentiation and dysfunction. Tbx5 is a member of the T-box gene family, which acts as transcription factor involved in the vertebrate heart development. Alteration of Tbx5 level affects the expression of hundreds of genes. Haploinsufficiency and gene duplication of Tbx5 are at the basis of the cardiac abnormalities associated with Holt–Oram syndrome (HOS). Recent data indicate that miRNAs might be an important part of the regulatory circuit through which Tbx5 controls heart development. Using high-throughput technologies, we characterized genome-widely the miRNA and mRNA expression profiles in WT- and Tbx5-depleted zebrafish embryos at two crucial developmental time points, 24 and 48 h post fertilization (hpf). We found that several miRNAs, which are potential effectors of Tbx5, are differentially expressed; some of them are already known to be involved in cardiac development and functions, such as miR-30, miR-34, miR-190, and miR-21. We performed an integrated analysis of miRNA expression data with gene expression profiles to refine computational target prediction approaches by means of the inversely correlation of miRNA–mRNA expressions, and we highlighted targets, which have roles in cardiac contractility, cardiomyocyte proliferation/apoptosis, and morphogenesis, crucial functions regulated by Tbx5. This approach allowed to discover complex regulatory circuits involving novel miRNAs and protein coding genes not considered before in the HOS such as miR-34a and miR-30 and their targets.

Keywords: zebrafish, heart, microRNA, NGS, microarray, Holt–Oram, data integration

# 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 oneto-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.11 , 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 (<IP_ADDRESS>), 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

# 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.3 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.

# 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.

2http://www.bioconductor.com 3http://genome.ucsc.edu/

1http://www.bioinformatics.babraham.ac.uk/projects/fastqc/

# 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.

# 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

*a Data for miR-19a-3p comes from our previous published data in (Chiavacci et al., 2015).* 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).

# 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.

# 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.4

# 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.

# 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 http://journal.frontiersin.org/article/10.3389/fbioe.2016.00060

4http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE64466

# REFERENCES

**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.

*Copyright © 2016 D'Aurizio, Russo, Chiavacci, Baumgart, Groth, D'Onofrio, Arisi, Rainaldi, Pitto and Pellegrini. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*


**1. Introduction**

Nanofluids attract current predilection because of its heat conduction attributes. Changing the flow geometry, boundary conditions, or thermal conductivity of liquids can improve convective heat transfer. Over the years, researchers have tried to increase the thermal conductivity of liquids. For this purpose, with the idea of Maxwell [1] solid metal particles are introduced into the base liquid. The large micro-sized particles are used to make suspensions because the conductivity of solids is greater than that of liquids, but these particles tend to produce greater resistance to the flow of base fluid. Modern nanotechnology tends to take a new direction in this field. In 1995, Choi [2] proposed a liquid with nano-sized particles suspended in a base liquid to eliminate the disadvantages of micro-sized particles. These liquids have efficient convective heat transfer compared to pure liquids. Recently, the idea of nanofluid in peristalsis has been studied by some researchers [3–10].

Peristalsis is characterized as the extension and the arrival of a substance into a liquid that improves the formative waves that broaden the length of the conduit, blending and shipping the liquid toward the wave spread. It is a mechanism that is available in numerous organs of the human body. In some specific instruments—for example heart–lung machines, implantation gadgets, and other pumping apparatus—such types of processes are utilized. It is of specific significance in many species and especially in human body that the transportation of many tissues of the body under various conditions, for example, the sucking of blood by leaches, the heap from the kidneys to the bladder through filtration, transport of the spermatozoa to the male genital tract, the development of the bosom in the Fallopian tubes, vasomotion of little veins, just as the blending and transport of gastrointestinal entry material.

The use of heat is of particular importance in the field due to its wide scope in engineering and biomechanics. In addition, the common relationship of heat stress and peristalsis can be observed during the oxygenation process with the patient. The assessment of magnetic resonance in biological tissues has aroused great interest among researchers regarding physical problems such as blood.

The assessment of heat transfer is related to the conditions of convection used in processes such as thermal conductivity, mechanical properties, chemical reactions, and so on. Aziz [11] presented a similarity solution to incorporate the convective walls conditions for thermal boundary layer on a smooth plate. In another article, Makinde and Aziz [12] developed the MHD mixed model on a flat surface in a concise way in terms of compatibility. Makinde [13] also discussed the flow of the MHD component with the temperature and the mechanical evaluation of a plate on a flat surface with extended conditions. Merkin and Pop [14] considered the analysis of heat transfer by dynamically simulating the flow of a uniform current on a flat surface with a horizontal displacement. According to them, the heat flux near the main edge is dominated by the surface heat flux.

After knowing the significance of the above discussed phenomena, authors are keen to develop a series solution of peristaltic flow of nanofluid with Williamson fluid model as a base liquid with convective boundary conditions travelling through asymmetric channel. At least we know that this study has not been yet explored in the literature. This study will be a good base for the engineers to utilize the results in procedures like thermal energy storage, gas turbines, nuclear workshops, etc. The problem is modeled under the induction of lubrication approach. The series solutions of stream function, temperature distribution, and nanoparticle concentration are achieved by using a well-known converging method the homotopy perturbation method. The important features are analyzed more specifically by sketching graphs to estimate the impact of pertinent constant physical factors.

#### **2. Mathematical Modeling**

The incompressible Williamson model is chosen as a base fluid for nanofluid in between an asymmetric channel experiencing heat convection at the peristaltic type surfaces. The width of the channel is taken as (*d*11 + *d*12). Flow is initiated due to the propagation of curved waves travelling with uniform speed *c* towards the flow. The exchange of heat is recognized by imposing temperatures *T*0 and *T*1 at the lower and upper areas, correspondingly. To discuss nano particle phenomenon, we have taken the nanoparticle concentration *C*0 and on the lower side and upper one, accordingly (see Figure 1).

**Figure 1.** Geometry of the channel.

The magnetic field *B*0 is exerted orthogonally. The wall surfaces are taken as

$$Y = H\_1 = d\_{11} + a\_{11} \cos[2\pi\lambda X'] \tag{1}$$

$$Y = H\_2 = -d\_{12} - b\_{11} \cos[2\pi\lambda X' + \overline{\varphi}], \text{ where } X' = X - ct \tag{2}$$

In upper defined equations, *a*11 and *b*11 represent the wave amplitudes, λ gives the wavelength, *t* suggests the time, *X* depicts the wave's direction, and *Y* is placed normally to *X*. The range of phase variance ϕ alters as 0 ≤ ϕ ≤ π. If ϕ = 0, we meant that a symmetric dimensional channel is having waves located out of the phase and ϕ = π, suggest the waves within the phase. Moreover *a*11, *b*11, *d*11, *d*12 and ϕ overcome the following relation

$$a\_{11}^2 + b\_{11}^2 + 2a\_{11}b\_{11}\cos\overline{\varphi} \le (d\_{11} + d\_{12})^2\tag{3}$$

The mathematical models of the considered problem given as

$$
\nabla \cdot \widehat{V} = 0 \tag{4}
$$

$$
\rho \left( \frac{\partial \widehat{V}}{\partial t} + \widehat{V} \cdot \nabla \widehat{V} \right) = -\nabla \overline{P} + \nabla \cdot \mathbf{S} + \rho\_f \mathbf{g} \, a\_f \left( \overline{T} - T\_0 \right) + \rho\_f \mathbf{g} \, a\_f \left( \overline{\mathbf{C}} - \mathbf{C}\_0 \right) + \mathbf{J} \times \mathbf{B} \tag{5}
$$

$$\mathbf{E}\left(\rho c\right)\_f \left(\frac{\partial \overline{T}}{\partial t} + \widehat{\boldsymbol{V}} \cdot \nabla \overline{T}\right) = \nabla \cdot \boldsymbol{K} \nabla \overline{T} + \mathbf{S} \cdot \boldsymbol{\nabla} \overline{\boldsymbol{V}} + (\rho c)\_p \left(D\_\mathcal{B} \Big(\nabla \overline{\mathbf{C}} \cdot \nabla \overline{T}\Big) + \frac{D\_\mathcal{T}}{T\_m} (\nabla \overline{T} \cdot \nabla \overline{T})\right) \tag{6}$$

$$
\frac{\partial \widetilde{\mathbf{C}}}{\partial t} + \hat{\mathbf{V}} \cdot \nabla \widetilde{\mathbf{C}} = \nabla \cdot \left( D\_B \nabla \widetilde{\mathbf{C}} + D\_T \frac{\nabla \widetilde{T}}{T\_o} \right) \tag{7}
$$

where **g** is the gravitational body force and α*f* represents the volumetric volume distension nanofluid's coefficient. In above relations, (ρ*c*)*f* denotes the fluid's heat capacity, (ρ*c*)*p* accounts for effective nanoparticles heat capacity, **J** = σ( *V* × **B**) reveals the current density, *B* = (0, *B*0) notifies the external magnetic field, and **S** is placed for the Cauchy stress tensor for Williamson fluid and is determined as

$$\mathbf{r} = \left(\mu\_{\diamond} + (\mu\_0 + \mu\_{\diamond}) \left(1 - \Gamma \dot{\mathbf{y}}\right)^{-1}\right) \dot{\overline{\mathbf{y}}} \tag{8}$$

where · γ comprises the subsequent value

$$\dot{\overline{\mathbf{y}}} = \sqrt{\frac{1}{2} \sum\_{i} \sum\_{j} \dot{\mathbf{y}}\_{ij} \dot{\mathbf{y}}\_{ji}} = \sqrt{\frac{1}{2} \overline{\Pi}} \tag{9}$$

Here Π is the strain tensor. The velocity profile for the given problem is considered as *V* = (*U*- , *V*- ). Introducing a wavy frame we introduce the following transformations

$$\mathbf{x} = \mathbf{X}', \ y = \mathbf{Y}, \ u = \mathsf{U}' - \mathsf{c}, \ v = \mathsf{V}', \ p(\mathbf{x}) = \overline{\mathsf{P}}(\mathbf{X}, \ t) \tag{10}$$

We suggest the following dimensionless parameters to be used in the above expressions

*x* = 2π*x* λ , *y* = *y d*11 , *u* = *u c* , *v* = *v c*δ , δ = *d*11 λ , *d* = *d*12 *d*11 , *p* = *d*2 11*p* μ*c*λ , *h*11 = *H*11 *d*11 , *h*12 = *H*12 *d*12 , *a*12 = *a*11 *d*11 , *Br* = *Ec*Pr, *b* = *b*11 *d*11 , Re = ρ*cd*11 μ , ψ = ψ *cd*11 , θ = 1*T*−*T*0 *T*1−*T*0 , *Ec* = *c*2 *cp*(*T*1−*T*0), Pr = ρν*c K* , *S* = *Sd*11 μ*c* , *We* = Γ*c d*11 , ϕ = *C* 1−*C*0 *C*1−*C*0 , *Gr* = ρ*f g*α*f d*2(*T*1−*T*0) *c*μ , *Gc* = ρ*f g*α*f d*2(*C*1−*C*0) *c*μ , *Nb* = τ*DB*(*C*1−*C*0) ν , *Nt* = τ*DT*(*T*1−*T*0) τ*Tm* , *M* = σ μ *B*0*d*11, *We* = Γ*c d*11 (11)

where *M*, *We*, *Br*, Pr, *Nb*, *Nt*, *Gr*, and *Gc* represent the Hartman number, Weissenberg number, Brinkman number, Prandtl number, Brownian motion parameter, thermophoresis parameter, local temperature Grashof number, and local nanoparticle Grashof number, accordingly. After incorporating the above structured parameters and applying the conditions of large wavelength along with small Reynolds number in a wavy frame coordinates we have the final form of Equations (4)–(7)

$$\begin{aligned} \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} &= 0\\ \frac{dp}{dx} = \frac{\partial}{\partial y} \left[ \frac{\partial^2 \psi}{\partial y^2} - M^2 \psi + \mathcal{W} \mathcal{E} \left( \frac{\partial^2 \psi}{\partial y^2} \right)^2 \right] + \mathcal{G}\_\mathbf{r} \boldsymbol{\theta} + \mathcal{G}\_\mathbf{r} \boldsymbol{\phi} \end{aligned} \tag{12}$$

or

$$-\frac{\partial^2}{\partial y^2} \left[ \frac{\partial^2 \psi}{\partial y^2} + \mathcal{W} \epsilon \left( \frac{\partial^2 \psi}{\partial y^2} \right)^2 - M^2 \psi \right] + G\_r \frac{\partial \mathcal{O}}{\partial y} + G\_c \frac{\partial \psi}{\partial y} = 0 \tag{13}$$

$$\Pr\left[N\_b \frac{\partial \boldsymbol{\psi}}{\partial \boldsymbol{y}} \frac{\partial \boldsymbol{\partial}}{\partial \boldsymbol{y}} + N\_t \left(\frac{\partial \boldsymbol{\partial}}{\partial \boldsymbol{y}}\right)^2\right] + \frac{\partial^2 \boldsymbol{\partial}}{\partial \boldsymbol{y}^2} + Br\left[\left(\frac{\partial^2 \boldsymbol{\psi}}{\partial \boldsymbol{y}^2}\right)^2 + \mathcal{W}\boldsymbol{\epsilon}\left(\frac{\partial^2 \boldsymbol{\psi}}{\partial \boldsymbol{y}^2}\right)^3\right] = 0\tag{14}$$

$$\frac{\partial^2 \varrho}{\partial y^2} + \frac{N\_\text{f}}{N\_b} \frac{\partial^2 \theta}{\partial y^2} = 0 \tag{15}$$

where ψ is stream function satisfying the relations *u* = ∂δ/∂*y* and *v* = −δ∂ψ/∂*x*, The no-slip boundary conditions for velocity *u* and nanoparticles fraction ϕ and convective boundaries are taken into consideration for temperature θ which have the following dimensionless form in the wave frame [15]

$$\begin{aligned} \psi = \frac{\mathbb{F}}{2}, \quad \frac{\partial \psi}{\partial y} = -1, \text{ at } y = h\_{11}, \; \psi = -\frac{\mathbb{F}}{2}, \; \frac{\partial \psi}{\partial y} = -1, \; \text{at } y = h\_{12}, \\\ \frac{\partial \theta}{\partial y} - B\_i \theta = -B\_i \text{ at } y = h\_{11} \text{ and } \theta = 0 \text{ at } y = h\_{12} \end{aligned} \tag{16}$$

$$\begin{aligned} \varphi = 1 \qquad \text{at } y = h\_{11} \quad \text{and} \quad \varphi = 0 \text{ at } y = h\_{12}. \end{aligned} \tag{16}$$

where *h*11 = 1 + *a*12 cos *x* and *h*12 = −*d* − *b* cos(*x* + ϕ). Also *Bi* = *hf d*11/*K* is the Biot number, *hf* stands for the coefficient of convective thermal transport. The mean flow rate in dimensionless format is elaborated as

$$Q = F + 1 + d \tag{17}$$