Circulating microRNAs are associated with Pulmonary Hypertension and Development of Chronic Lung Disease in Congenital Diaphragmatic Hernia

Pulmonary hypertension (PH) contributes to high mortality in congenital diaphragmatic hernia (CDH). A better understanding of the regulatory mechanisms underlying the pathology in CDH might allow the identification of prognostic biomarkers and potential therapeutic targets. We report the results from an expression profiling of circulating microRNAs (miRNAs) in direct post-pulmonary blood flow of 18 CDH newborns. Seven miRNAs differentially expressed in children that either died or developed chronic lung disease (CLD) up to 28 days after birth, compared to those who survived without developing CLD during this period, were identified. Target gene and pathway analyses indicate that these miRNAs functions include regulation of the cell cycle, inflammation and morphogenesis, by targeting molecules responsive to growth factors, cytokines and cellular stressors. Furthermore, we identified hub molecules by constructing a protein-protein interaction network of shared targets, and ranked the relative importance of the identified miRNAs. Our results suggest that dysregulations in miRNAs let-7b-5p, -7c-5p, miR-1307-3p, -185-3p, -8084, -331-3p and -210-3p may be detrimental for the development and function of the lungs and pulmonary vasculature, compromise cardiac function and contribute to the development of CLD in CDH. Further investigation of the biomarker and therapeutic potential of these circulating miRNAs is encouraged.

Circulating miRNAs associated with severe outcomes in CDH. Through array profiling of mature miRNAs in blood of CDH newborns, collected 24 h after birth, we identified 33 circulating miRNAs that were significantly changed (p < 0.05) in those children who died or developed CLD within 28 days (Death/CLD group) (Fig. 1A, Supplementary Table 2), compared to those who survived with no signs of CLD until up to 28 days of life (No-CLD group). Considering the reduced size of our sample, we used a combined endpoint of death and CLD as a representation of severe disease outcomes. Seven of the 33 significant miRNAs survived correction for multiple comparisons (BH-adjusted p < 0.05, Table 2, Fig. 1B) with an expression change of at least 1, and were thus considered differentially expressed (DE). We observed down-regulation of most of the significant miRNAs in the Death/CLD group. Most miRNAs also appear to have, in average, a reported medium expression in the lung, with considerable numbers of target mRNAs and proteins in this organ, according to the IMOTA database. Interestingly, miR-6511b-3p (r = 0.621) and miR-25-5p (r = −0.77) showed a relatively high correlation with death, although the number of death cases in our study is too low to consider this a truly significant finding, warranting further exploration. Additional potentially interesting relationships include the negative correlation of miR-8084 expression with moderate to severe PH (r = −0.549), using an established classification of PH severity in this population, and the positive correlation of miR-744-5p levels with PaO 2 measures at 24 h (r = 0.741) ( Fig. 2A). As expected, most significant miRNAs correlated with the outcome, represented as either the binomial No-CLD or Death/CLD classification, or as a graded No-CLD → CLD → Death classification. Heatmap of significantly changed circulating mature miRNAs in CDH-PH children who died or developed CLD during the first 28 days after birth (Death/CLD at 28 days, n = 8), compared to those children who survived without developing CLD during this period (No-CLD at 28 days, n = 10). Blood samples were drawn 24 h after birth. Statistical significance was obtained through a moderated t-test (unadjusted p < 0.05). Hierarchical clustering was applied to samples and miRNAs using Euclidean distances with complete linkage. Furthermore, some DE miRNAs mildly correlated as well with the use of ECMO support and PaO 2 , including let-7b-5p (r ECMO = 0.727, r PaO2 = −0.601), let-7c-5p (r ECMO = 0.755, r PaO2 = −0.64) and miR-8084 (r ECMO = −0.693, r PaO2 = 0.537). However, it is difficult to tell whether these relationships could be the result of ECMO-induced changes already present by the time of sample collection, or could truly result from the abnormal biological processes that ultimately determined the need for ECMO support.
Four major clusters of significant miRNAs can be appreciated from the correlation patterns amongst their expression levels, with the seven DE miRNAs distributed in three of these clusters (Fig. 2B), which might relate to differential functions between clusters and/or to different sources of the circulating miRNAs.

Functional implications of miRNA dysregulation in the development of CLD in CDH children.
We desired to use the miRSystem online tool to investigate target genes of our DE miRNAs, because of the advantages provided by its integrative analysis. However, only four of the seven DE miRNAs (let-7b-5p, let-7c-5p, miR-210-3p and miR-331-3p) were present in this system, for which the target genes for the remaining DE miR-NAs (miR-185-3p, miR-1307-3p and miR-8084) were obtained from independent queries in miRTargetLink. Our integrated list, consisting of 857 genes with validated status or a minimum of three hits from miRSystem, and genes from miRTargetLink with some degree of evidence, resulted in the identification of 411 genes predicted to be targeted by more than one of our seven DE miRNAs (Supplementary Table 3). No common target for all DE miRNAs was identified; the largest hits were observed for SEMA4G, which was predicted to be targeted by four DE miRNAs, while the highest observed/expected (O/E) ratios, based on the miRSystem integrative analysis, were observed for NRTN and ESPL1, targeted by the let-7 family of DE miRNAs (Table 3).
We further explored the functional interactions between common targets of DE miRNAs, as well as their functional implications, by constructing a PPI network and performing pathway analyses on it. This allowed us to identify not only functional types of interaction (i.e. activation/catalysis, inhibition, complex formation/ binding, or predicted interaction), but also functionally related target clusters (network modules), as well as proteins likely linking different biological processes/pathways (network hubs, i.e. proteins showing the larger numbers of interactions within the network). The resulting full network consisted of 227 nodes and was enriched for pathways highly related to responses to growth factors, adhesion molecules and immunological stimulation, as well as morpho-/organogenesis, including the PI3K-Akt, MAPK, FoxO, Ras, IL4-mediated, TGF-β, integrin, Jak-STAT and p53 signaling pathways, focal adhesion and ECM-receptor interactions, among the most significantly enriched terms, and cell cycle regulation by BTG proteins, signaling by activin, Fas, SMAD2/3, EPHA, growth hormone and TNF receptors, regulation of hypoxia and oxygen homeostasis by HIF-1-α, and oxidative stress response, among the pathways with the higher proportions of proteins from the network in the pathway (Supplementary Table 4, Module: All).
The main network (Fig. 3, Supplementary Table 5) was conformed by 208 nodes, 480 edges (interactions), and 13 modules (clusters). As network hubs, we identified MAPK1 (38 interactions), STAT3 (31 interactions), ITGB3 and SMAD2 (22 interactions). In a lesser extent, NRAS and TP53 (18 interactions), as well as CHUK and PDGFB (16 interactions), were also highlighted. With the exception of MAPK1, which was a predicted target of miR-1307-3p and miR-185-3p, these network hub molecules were predicted targets for let-7b/c-5p miRNAs ( Table 3). Pathway enrichment analysis of the network modules revealed that the largest clusters of proteins are involved in the responses to immune stimulation and growth signals, cell cycle regulation and tissue/organ morphogenesis, while smaller clusters play roles in hypoxia, oxytocin signaling, platelet homeostasis and cardiac conduction ( Table 4, Supplementary Table 4). All these modules are closely linked through the aforementioned hub proteins. To illustrate this functional connectivity further, we mapped those miRNA targets observed within the Pathways in Cancer KEGG diagram as, due to the integrative nature of this pathway and its relation to cellular proliferation, growth and differentiation, we observed the highest proportion of proteins from the network in this term ( Supplementary Fig. 1). In summary, the results from our analysis support the involvement of the TGF-β family of growth and differentiation factors, inflammation and semaphorin signaling, which may impair organ fetal organogenesis and postnatal lung and cardiac functions, in the severe pathophysiology of CDH with poor outcomes. Circulating DE miRNAs as potential prognostic biomarkers and/or therapeutic targets for CLD in CHD children. Finally, because of the high overlapping and redundant nature of miRNA target genes and pathways, we wished to investigate the relative importance of dysregulations at 24 h of the circulating DE miR-NAs for the outcomes of Death/CLD at 28 days after birth in CDH children, in order to rank our DE miRNAs for future studies on their potential as prognostic biomarkers for CLD in CDH, and/or therapeutic targets. For this, we used C&RTs and found that most Death/CLD cases (7/8) could be differentiated from those No-CLD cases through the expression values of miR-1307-3p (importance rank: 79, split constant: 5.31), while miR-185-3p was important to separate the remaining set of Death/CLD cases initially classified as No-CLD within the first split (importance rank: 100, split constant: 3.305) to achieve a complete (100%) classification accuracy (Fig. 4A), using the lower cost (cross-validation cost: 0.167, resubstitution cost: 0), most suitable tree (Fig. 4B). From this algorithm, we obtained then a relative higher ranking of miR-185-3p and miR-1307-3p, closely followed by miR-210-3p (importance rank: 77), and the let-7b/c-5p miRNAs (importance rank: 76) (Fig. 4C), even though these latter were not included in the split criteria.

Discussion
In this study, we identified seven miRNAs differentially expressed at 24 h of life in the blood of CDH newborns who developed CLD or died within 28 days after birth, and showed that these not only correlate with different measures of pulmonary dysfunction, but also interconnect to participate in processes crucial to establish proper morphology and function of the heart and lungs. Our results suggest that a subset of these miRNAs might hold the potential to serve as prognostic circulating biomarkers for the development of CLD in CDH newborns, which should be properly evaluated in an extended cohort, and implicate TGF-β and semaphorin signaling, as well as inflammatory responses, in the development of CLD in CDH. The rate of mortality and disability in the long-term among CDH patients depends, in great extent, on the presence and severity of PH 19 . Persistent problems in surviving children, such as long-term oxygen dependency, need for mechanical ventilation, persistent wheezing, and increased risk for pulmonary infections [20][21][22] , suggest that a combination of mechanical ventilation-induced damage and primary defects of the lung and pulmonary vascular structures are responsible for the cardiopulmonary abnormalities observed in CDH children 2 . The contributions of miRNAs to lung morphogenesis during fetal development, and to inflammatory processes during pulmonary disease, have been attractive topics for research in past years. This led to the identification of a number of these molecules, including miR-210 and let-7 family members, as important regulators of normal and abnormal lung development and function 23 , and of pulmonary arterial remodeling in response to vessel injury and cellular stressors, such as hypoxia 24 . Consistently, members of the let-7 family of miRNAs and miR-1307 have been implicated in the severity of PH in systemic scleroderma 25 , miR-185 has been shown to contribute to the oxidative stress-mediated, hyperoxia-induced death of lung epithelial cells during acute lung injury and acute respiratory distress syndrome 26 and, together with miR-210, contribute to the rapid progression of pulmonary fibrosis 27 . Further, miR-210, together with miR-331, has been associated with cardiovascular disease and markers of systemic inflammation in human immunodeficiency virus (HIV)-1 infection 28 .
The hypoxia-induced miR-210 has been previously implicated in PH, other lung diseases and ischemic heart disease [29][30][31][32] . We suggest that this miRNA also contributes to the development of CLD and the severity of PH as outcomes in CDH newborns. Interestingly, though, miR-210 presented a lower importance than miR-185 and miR-1307 to differentiate CDH cases within the Death/CLD group from those in the No-CLD group. We can speculate that these latter miRNAs regulate distinctive pathways that contribute to progression of CLD in a  greater extent than the former, which might be implicated in pathways contributing to PH pathogenesis. Similarly, the previously reported PH-associated miR-21, miR-204, and clusters miR-143/145, miR-17-92 and miR-130/301, as well as the CDH-associated miR-200b and miR-10a might not have shown significance in our study because of a reduced relevance for severity and progression of CLD in CDH newborns, compared to the pathogenesis of each individual disease. Moreover, most observations on miRNA dysregulations in CDH and PH come from studies in animal or cellular models, which may not fully translate to the human pathobiology. It is well-known that the TGF-β superfamily plays crucial roles in pre-and post-natal lung development, importantly shaping alveolarization and controlling the extracellular matrix composition and tissue homeostasis, among other functions. Not surprisingly, its involvement in pulmonary and cardiovascular diseases has been largely described [33][34][35][36] . Because the TGF-β canonical and non-canonical signaling pathways overlap with a number of pathway terms highly enriched in our PPI network 34,37 , our results provide evidence for the involvement of the TGF-β superfamily in the development of CLD in CDH newborns, providing support to the hypothesis that a primary disruption in lung vasculature and airway development accompany the diaphragmatic defect in CDH, and that the extent of the resulting cardiopulmonary dysfunction contributes to the development of CLD in CDH newborns. Further support to this hypothesis is found through our observation that major genes Figure 3. Protein-protein interaction network generated for shared miRNA target genes. After integration of miRSystem and miRTargetLink target gene lists, those genes predicted to be targeted by more than one of the differentially expressed miRNAs served as input for the network creation using the Reactome FI app for Cytoscape. Disconnected nodes and isolated clusters of less than 4 proteins were excluded from the network view and module enrichments. Information on paired protein interactions is available in the Supplementary Table 3. Top results from the pathway enrichment analyses for the network and network modules can be found in Table 3, and Supplementary Table 4 participating in the semaphorin signaling, involved in development and regulation of immune responses through guidance cues 38,39 , are predicted targets of the miRNAs we found dysregulated in blood of CDH children with severe outcomes.
The O/E ratio of fetal lung volume measured by ultrasound or magnetic resonance imaging (MRI) has been proposed to serve as prognostic biomarker for mortality, the need for ECMO and development of CLD 40,41 . Our study is, to our knowledge, the first to assess the differential miRNA blood profiles related to death and the development of CLD in CDH, and provide an extensive analysis of functional implications and hierarchies of miRNA dysregulations which not only suggests potential for measurements of circulating miRNAs as prognostic biomarkers, but also uncovers crucial regulators with potential for therapeutic purposes. It has already been demonstrated that less severe lung hypoplasia in CDH is associated with compensatory upregulation of miR-200b, and that treatment with this miRNA at an early developmental stage improves lung growth in the nitrofen-induced CDH rodent model 42 .
We acknowledge the small sample size of our study is an important limitation, providing small power for statistical testing. Our miRNA candidates would require validation in a larger, well-characterized cohort. Another limitation in our study is the lack of repeated measures in a time series, which would have allowed us to identify changes induced by ECMO and other external factors that influence the development of CLD in CDH newborns. Similarly, the unavailability of tissue-specific miRNA expression limits our ability to interpret the biological meaning of our pathway and network analyses. However, all of these pathways have been described in the context of lung development, PH, and inflammation, all being a part of the pathophysiology of CDH. Considering array measures of small non-coding RNAs may be biased only to those miRNAs and other small RNAs that are part of the array, it is possible that using RNA-seq would have led to discovery of additional, novel circulating miRNAs associated with severe outcomes in CDH. This is an important avenue of future research. Finally, we also acknowledge the potential overfitting of our C&RTs and, although we obtained similar results using two different algorithms (data not shown) for this same reason, the aim of this analysis was not to provide proof of the biomarker performance of our identified miRNAs, but to help us in ranking them for future follow-up studies.
In conclusion, we identified several miRNAs differentially expressed at 24 h of life in the blood of CDH newborns who developed CLD or died within 28 days after birth. Whether these miRNAs might serve as prognostic circulating biomarkers for the development of CLD in CDH newborns, and might be useful therapeutic targets, needs to be investigated in a larger cohort. Despite the aforementioned limitations, we believe our study provides encouraging evidence to ensure further validation of our present findings.

227
Signaling pathways activated in response to growth factors and adhesion molecules; tissue and organ morphogenesis.

Statistical analysis.
Statistical comparisons of demographic/clinical data and circulating miRNA expression levels between children in the Death/CLD and No-CLD groups was performed by means of two tailed Mann-Whitney U-tests, and a moderated t-test adjusted for multiple comparisons through the Benjamini-Hochberg (BH) method, respectively. Statistical significance was set to raw p < 0.05 values for demographic/clinical data and "significantly changed" miRNAs, whereas "differentially expressed" (DE) miRNAs were considered those with BH-adjusted p < 0.05 and expression change ≥1.0. Correlation coefficients between clinical parameters and expression levels of DE miRNAs were obtained through the Pearson's method.
Target genes. DE miRNAs were submitted for "miRNAs to Target Genes" analysis on miRSystem 16 , using default settings and the expression changes between Death/CLD and No-CLD groups as weights. DE miRNAs not found in the miRSystem database were queried individually on miRTargetLink 17 to obtain the evidence-based target genes. Overlaps between individual target lists were then identified. Furthermore, it was investigated whether the significantly changed miRNAs are expressed in lung tissue and/or there is evidence that they have targets in lung, by searching each miRNA in the interactive multi-omics-tissue atlas (IMOTA) 18 .

Protein-protein interactions and biological pathways.
To investigate the functional implications of miRNA dysregulations in severe CDH outcomes, we created a protein-protein interaction (PPI) network of molecules targeted by at least two of our DE miRNAs with pathway annotations in miRSystem using the Reactome FI app for Cytoscape 3.5.1. We clustered the network and fetched the annotations to show the types of functional relationships. All unconnected nodes, as well as clusters conformed by less than four proteins, were removed from the network view. The full network and network modules (clusters) not hidden from view were analyzed for enriched biological pathways. Enrichment results were filtered to keep only those terms with false discovery rate (FDR) <0.01 and a number of module genes in the term > 2.