Tracheal aspirate miRNA signatures in preterm infants with severe bronchopulmonary dysplasia

Bronchopulmonary dysplasia (BPD) is a form of chronic lung disease that develops in neonates as a consequence of preterm birth and arrested fetal lung development. The incidence of BPD remains on the rise, as a result of increasing survival of extremely preterm infants. Severe BPD contributes to significant health care costs and is associated with prolonged hospitalizations, respiratory infections, and neurodevelopmental deficits. In this study, we aimed to detect novel biomarkers of severe BPD. We collected tracheal aspirates (TA) from preterm babies with mild/moderate (n = 8) and severe (n = 17) BPD, and we profiled the expression of 1048 miRNAs using a PCR array. Associations with biological pathways were determined with the Ingenuity Pathway Analysis (IPA) software. We found 31 miRNAs differentially expressed between the two disease groups (2-fold change, FDR < 0.05). Of these, 4 miRNAs displayed significantly higher expression levels, and 27 miRNAs had significantly lower expression levels in the severe BPD vs. the mild/moderate BPD group. IPA identified cell signaling and inflammation pathways associated with miRNA signatures. We conclude that TAs of extreme premature infants contain miRNA signatures associated with severe BPD. These signatures may serve as biomarkers of disease severity in infants with BPD.


Introduction
Bronchopulmonary dysplasia (BPD) is a form of chronic lung disease that develops in neonates as a consequence of preterm birth and arrested fetal lung development [1]. BPD is associated with significant human and public health burden, with reported cases varying significantly among centers world-wide [2]. The incidence of BPD remains on the rise with recently stated global incidence ranging from 10% to up to 89% in one report [2,3]. The disease was first described by Northway and colleagues in 1967 in preterm infants receiving invasive mechanical ventilation with characteristic clinical, chest radiological, and pathologic findings [4].
The etiology and pathogenesis of BPD is complex. The contributing factors are multi-factorial and include genetic predisposition, epigenetic factors, arrest of lung development, chronic inflammation, mechanical ventilation, and oxygen toxicity [1,5]. Severe BPD is particularly prevalent in extremely low birth weight infants, a fast-growing population. Fetal growth restriction and male gender are additional risk factors for severe BPD [5,6]. Neonates with severe BPD are at risk for adverse short-and long-term outcomes that are potentially life-long [7]. Preterm infants who are

Tracheal aspirate collection
Following informed consent by the parents, we collected TAs from enrolled patients following routine suctioning. We also obtained pertinent clinical information from electronic medical records, such as sex, race/ethnicity, method of delivery, antenatal steroids, gestational age (GA) at birth, birth weight, GA at the time of sample collection, day of life, and fraction of inspired oxygen (FiO2) at the time of sample collection. Immediately after collection, TAs were transported to the laboratory and stored at -80°C until analysis.

MicroRNA purification
Starting from 500l of TAs, we purified miRNAs using the Norgen miRNA Purification Kit, after addition of a spike-in control (cel-mir-39, QIAGEN). The miRNA fraction was then quantified using a Nanodrop, and its purity was assessed with a Bioanalyzer at the Penn State College of Medicine Genome Sciences Core Facility. For miRNA profiling, a total of 250 ng of RNA were first retrotranscribed using the miScriptII RT kit (QIAGEN), following the manufacturer's protocol.

MicroRNA arrays
The expression of 1,066 human miRNAs was quantified using the miScript miRNA PCR array kit (QIAGEN, MIHS-3216Z) on a QuantStudio 12K Flex system (Life Technologies, Carlsbad, CA). Following the PCR reaction, cycle threshold (Ct) values were extracted, and normalized by global means. Differential expression was calculated using the 2 -ΔΔCT method and a mild/moderate sample as a calibrator [23]. Fold changes in miRNA expression were obtained using the Bioconductor limma package on R [24]. MiRNAs with undetectable expression (Ct values > 40) in more than 90% of samples of each experimental group were excluded from the analysis. Datasets (metadata, nonnormalized, and normalized values) were uploaded to the Gene Expression Omnibus under project GSE165828 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE165828).

Ingenuity pathway analysis
The IPA software (QIAGEN, www.qiagen.com/ingenuity) was used to identify miRNA interaction networks, predicted target genes, and molecular functions based on prediction scores [25]. The IPA core analysis functionality was performed for both direct and indirect relationships using the Ingenuity Knowledge Base. The miRNAs of interest were analyzed using the microRNA target filter feature. Only experimentally observed or predicted with high confidence targets in humans were considered.

Statistical analysis and code availability
For miRNA arrays, statistical analyses were performed using the Bioconductor limma package on R using the parametric empirical Bayes method [24]. Differential expression was defined as Benjamini-Hochberg False Discovery Rate (FDR) below 0.05. Heatmaps were generated with the nonnegative matrix factorization (NMF) package on R [27], and volcano plots were generated using the GraphPad Prism software. For miRNA validation assays, differences in miRNA expression between the two groups were determined by t-test using the GraphPad Prism software and considered statistically significant at p<0.05.

Code availability
The code used for data analysis in the current study can be found at the Silveyra lab repository, http://psilveyra.github.io/silveyralab/

Patient Demographics
We obtained TAs from 25 mechanically ventilated preterm infants with mild/moderate BPD (n=8) and severe BPD (n=17). A summary of the patient's clinical information is shown in Table 1, and individual data for main variables are visualized in Figure 1. As expected, FiO2 at the time of sample collection was significantly higher in the severe BPD vs. the mild/moderate BPD group (p<0.05). No statistically significant differences were found for the rest of the measured variables, including GA at birth and at sample collection time, birth weight, day of life at sample collection time, male sex, antenatal steroid exposure. Regarding the racial and ethnic composition of the patient cohort, there was a significant overrepresentation of Hispanic and mixed-race subjects in the severe BPD group than the mild/moderate BPD group (p<0.05).  (3) 0.070 More than one race 0 (0) 6 (1) 0.008 1 Any antenatal steroids exposure (1 or more doses received by mother before birth)

MiRNA expression in TAs
A heatmap of the expression of the 1048 miRNAs expressed in TA samples from mild/moderate and severe BPD infants is shown in Figure 2. A volcano plot indicating FDR vs fold change (FC) of expression between severe BPD and mild/moderate BPD samples is also shown in Figure 3.  After accounting for FDR < 0.05, the PCR array analysis revealed differential expression of 31 miRNAs with |FC| > 2 between the mild/moderate and severe BPD samples (Table 2). Of these, 4 miRNAs had significantly higher expression in the severe BPD group (log FC > 0.3, FDR < 0.05) and 27 had significantly higher expression in the mild/moderate BPD group (log FC < -0.3, FDR < 0.05) ( Table 2).

Pathway analysis
Ingenuity pathway analysis (IPA) of differentially expressed miRNAs in the two BPD groups indicated significant associations with molecular and cellular functions including cell signaling, DNA replication, and cell cycle, as well inflammatory response and disease, summarized in Table 3. We also found that the top regulatory networks associated with the differentially expressed miRNAs included cell cycle, control of gene expression, cell signaling, and cellular function, assembly and maintenance (Table 4). Of the 31 miRNAs identified in the array, 5 were associated with the most significant cellular and molecular networks found by IPA. These include 3 miRNAs with significantly lower expression in severe vs. mild/moderate BPD: miR-15* (mir-15a-3p), miR-615-3p, and miR-1255-5p, and 2 miRNAs with significantly higher expression in severe BPD: miR-185-5p (miR-185) and miR-9718 (which has the same seed sequence as miR-628-5p: UGCUGAC) ( Table 4). A diagram of the identified networks, miRNA-gene targets interactions, and connections among different networks based on experimental evidence is shown in Figure 5. The miRNA targets identified by IPA include cyclins, transcription factors, hormone receptors, and cell signaling enzymes (Table 4).  Next, to expand our analysis beyond experimentally validated interactions, we used the IPA miRNA target filter to identify predicted miRNA-mRNA target relationships, based on miRNA seed sequence foreseen interactions. Among the top networks associated with predicted target genes, we identified organismal injury and respiratory disease pathways, as well as cellular trafficking, immunological disease, and cell and tissue development gene networks (Table 5). Graphical representations for the top associated networks are also shown in supplementary figures S1 and S2.

Discussion
Many studies have focused on prevention of severe BPD, but there is paucity of literature on the biological basis of this disease's pathogenesis. There is also a lack of validated biomarkers for BPD severity and targeted pharmacological therapies to improve the disease burden [28]. The current recommendations include optimizing ventilator settings while focusing on nutrition and addressing comorbidities [29]. This mostly stems from the fact that exact underlying mechanisms of severe BPD are uncertain. In the current study, we sought to identify the molecular signatures associated with BPD severity, to uncover regulatory pathways and new biomarkers for severe BPD. Through analysis of TAs, a non-invasive method of sample collection in mechanically ventilated infants, we discovered specific miRNA profiles in age-matched preterm newborns with mild/moderate and severe BPD. Our analysis revealed 31 differentially expressed miRNAs between subject groups. Of these, 4 miRNAs displayed at least 2-fold higher expression and FDR < 0.05 in severe BPD infants, and 27 miRNAs had significantly lower expression in than in TAs from mild/moderate BPD infants. We also found associations of these miRNAs with functional gene networks related to lung inflammation, respiratory disease, oxidative stress, and cellular development, indicating that specific regulatory pathways may contribute to fundamental mechanisms of BPD severity.
Of the four miRNAs with higher expression in severe BPD infants, two (hsa-miR-545* and miR-185) have been previously implicated in mechanisms of cell cycle arrest, a known effect of prolonged mechanical ventilation and hyperoxia [30,31]. Specifically, miR-545 has been shown to suppress cell proliferation through targeting of cyclins and cyclin-related kinases [32], and miR-185 has been associated with induction of G1 cell cycle arrest and promotion of necroptosis and apoptosis via receptor-interacting kinases and caspase activity in alveolar epithelial type II cells [33,34]. Moreover, the role of mir-185 in hyperoxia-induced DNA damage and lung epithelial cell death triggered by oxidative stress has been reported using both animal models and human cells [34,35]. Importantly, extracellular vesicles containing miR-185 are significantly elevated following hyperoxia-induced cell death in alveolar type II epithelial cells [34]. While we did not characterize the composition of TAs regarding extracellular vesicles, it is possible that the increased levels of miR-185 found in TAs from severe BPD infants (who are subjected to higher and more prolonged hyperoxia than mild/moderate BPD infants) could represent a signal of hyperoxia-induced epithelial cell damage via increased extracellular vesicle miR-185 cargo. Therefore, our results could provide clinical evidence of a mechanism previously described in lung disease models.
Two additional targets of miR-185 are the Rho GTPases, CDC42 and RHOA (Table 4), implicated in cell cycle, cellular function and maintenance, and cell signaling. By targeting these mRNAs, miR-185 has been shown to inhibit cell proliferation and migration in hepatic, colorectal, and lung cancer cell models [36][37][38][39]. Importantly, miR-185 can also interact with and target the long-noncoding RNA FOXD2-AS1, known to inhibit cell proliferation, migration, and epithelial-mesenchymal transition (EMT) in glial cells and renal fibrosis models [40,41]. Regarding lung cells, Lei et al. demonstrated that miR-185 can reduce collagen V overexpression and EMT in pulmonary fibrosis [42]. While promotion of EMT in alveolar epithelial type II cells is known to alter normal alveolar development processes and contribute to BPD phenotypes, the effects of miR-185 in EMT processing during lung development and BPD progression have not been yet explored [43,44]. Finally, miR-185 has also been found to target the SOX9 gene and regulate Wnt signaling [45], a key regulator of lung organogenesis and development [46,47].
The remaining two miRNAs found overrepresented in TAs of severe BPD infants were miR-378 and miR-628-5p. With the exception of lung cancer [48], no studies to date have assessed the role of miR-378 in the lung. However, this small RNA has been implicated in several functions associated with BPD, including mitochondrial metabolism, autophagy, and oxidative stress responses in a variety of tissues [49,50]. On the other hand, miR-628-5p has been implicated in developmental functions such as embryonal implantation, as well as mechanisms of immune modulation, viral infection, and cell proliferation [51][52][53]. One study identified the fibroblast growth factor receptor 2 (FGFR2), the receptor for FGF10, as a target of miR-628-5p in ovarian cancer cells [54]. This gene has been identified as a key factor in early lung development [55], as well as alveolar epithelial type II cells homeostasis [56]. Interestingly, miR-628-5p has also been postulated as potential biomarker for a variety of cancers in adults [57]. While no studies have yet assessed a direct role of this miRNA in lung development and disease, it is possible it could serve as a biomarker for severe BPD in neonates.
Our analysis identified 27 miRNAs with low expression in severe vs. mild/moderate BPD. Of these, 3 were implicated in the top significant pathways identified by IPA (hsa-miR-15a*, hsa-miR-1255b-5p, and hsa-miR-615-3p). Specifically, these miRNAs were implicated in cell cycle and cell signaling functions, as well as disease processes in extrapulmonary tissues (Table 4). Regarding miR-15a, an avian model of long-term hypoxia stress revealed that its levels are influenced by hypoxiainduced factor 1 (HIF-1) during lung development [58]. Hypoxic episodes are common stressors associated with BPD secondary to immature respiratory center drive. Both acute and chronic hypoxic events associated with reduced peripheral chemosensitivity and pro-inflammatory responses occur during the course of neonatal BPD, and are associated with poorer neurodevelopmental outcomes [59] . In addition, miR-15a has been shown to inhibit lung fibrosis by targeting the yes-associated protein 1 (YAP1), a key downstream effector of the Hippo pathway [60]. Therefore, it is possible that lower expression of this miRNA in severe BPD infants is associated with fibrotic phenotypes in severe BPD. It is worthwhile to mention that the Hippo pathway is a powerful regulator of lung development, cell differentiation, and tissue regeneration and homeostasis, and has been found dysregulated in lymphocytes of patients with BPD [61,62]. Therefore, the role of miR-15a and the Hippo pathway in mechanisms leading to BPD severity warrants further investigation.
We detected and validated downregulation of miR-1255b-5p in TAs of severe BPD infants. This miRNA was previously noted to regulate vascular endothelial growth factor A (VEGFA) expression in liver cirrhosis and hepatocellular carcinoma [63]. Abnormal VEGF and disrupted angiogenesis are described in lung tissues of infants that died of severe BPD [64]. In addition, miR-1255b-5p is downregulated in hypoxia-induced lung A549 cells [65] as well as in models of high altitude retinopathy [66], and colorectal cancer, where it also regulates telomerase activity and suppresses EMT [67]. Notably, hypoxic stress in the presence of hyperoxia results in altering alveologenesis in mouse models of neonatal BPD [68].
Through IPA analysis, we found that the downregulated miRNA miR-615-3p directly targets the Ligand Dependent Nuclear Receptor Corepresor (LCOR) and Androgen Receptor (AR) transcripts, and indirectly targets the nuclear receptor PPARG (peroxisome proliferator activated receptor gamma) ( Table 4 and Figure 5). The AR is expressed in both male and female lung epithelial cells and can bind to DNA sequences to regulate the expression of cell differentiation genes [69,70]. Studies have also shown that the AR can mediate androgen-induced delays in fetal lung maturation [71][72][73]. For example, maternal treatment with dihydrotestosterone inhibits surfactant phospholipid production in the fetal lung, while treatment with the AR selective antagonist flutamide enhances surfactant phospholipid production [74]. Thus, downregulation of miR-615-3p (a negative regulator of AR expression) in severe BPD could result in increased levels of AR, contributing to delays in lung maturation and surfactant expression. Similarly, PPARG contributes to homeostasis maintenance of epithelial-mesenchymal interactions, which are key players of lung organogenesis. Disruption of this process results in trans-differentiation of lung alveolar lipofibroblasts to myofibroblasts, contributing to the development of pulmonary fibrosis [75,76]. Importantly, PPARG agonists have been suggested as therapeutic targets for BPD, due to their implications in the Wnt/-catenin and TGF- pathways [77]. Both pathways are induced by hypoxia, leading to decreased levels of PPARG and the subsequent differentiation of fibroblasts, leading to pulmonary fibrosis. Finally, because miR-615-3p is known to repress LCOR in macrophages [78], its downregulation in severe BPD may result in high LCOR levels, which also contribute to decreased PPARG. Overall, the miRNA signatures identified in severe BPD TAs can be associated with known mechanisms of lung inflammation, oxidative stress, developmental delay, and pulmonary fibrosis. More research is needed to explore their contributions to these mechanisms, and their therapeutic potential.
The strengths of this study include the prospective nature of the design, using a cohort of preterm infants that is matched in an age, sex, birth weight, method of delivery, and antenatal steroid exposure, and the identification of novel miRNA expression signatures associated with BPD phenotypes in TAs. An advantage of these samples is that they are readily available in intubated patients and may contain cellular makers that are very similar to those noted in the lungs. This study also excluded infants with pulmonary hypertension, which is a severe complication of BPD, hence providing more specificity to the signatures identified in the severe BPD cohort.
We have also identified the following study limitations. First, this was a very small but representative sample of preterm neonates admitted to a single center which affects the generalizability of the findings. It is also possible that with such a small sample size, the differential miRNA expression could be specific to the patients studied. However, the study design was exploratory in nature, utilizing convenient sampling in a small cohort of preterm infants receiving invasive mechanical ventilation. Second, given that the samples were collected at a single time point of the disease process, we are unable to determine whether the observed miRNA signatures are representative of the disease status or are only present at a specific time in the course of the pathology. This limits our ability to determine if there is a predictive potential to these markers, or whether they represent a reflection of ongoing molecular pathways inherent to the disease severity, or a trigger for molecular events related to disease phenotypes. Additionally, some key clinical variables were not fully explored in this study due to the small sample size. For example, the influence of sex differences, birth weight, fetal growth restriction, and exposure to antenatal corticosteroids on miRNA profile changes were not directly evaluated. Third, given that we studied preterm neonates of different developmental trajectories, it is possible that the differences we noted in the miRNA expression may reflect changes in cell populations present in the TA [79,80]. Lastly, it is becoming evident in recent studies that the microbiota in TAs can affect miRNA expression [81,82]. Because pulmonary inflammation secondary to perinatal and postnatal infection has been implicated as a key contributor in the pathogenesis of BPD, changes in the TA microbiome may be reflected in the observed miRNA signatures. Despite these limitations, this study demonstrates miRNA profiling in TAs of preterm neonates can be used successfully as a potential biomarker for BPD severity.
In summary, we describe miRNA expression signatures in TA samples obtained from extreme preterm infants with severe BPD, compared to mild to moderate BPD. We identified associations of these molecules with pathways involved in the disease phenotype, including lung inflammation, arrested cell division and development, and oxidative stress. Future studies are needed to evaluate the significance of miRNA regulation in severe BPD, as well as their potential use as predictive and diagnostic tools for disease severity and progression.

Supplementary Materials:
The following are available online at www.mdpi.com/xxx/s1, Figure S1: Respiratory disease, cancer, organismal injury and abnormalities network, Figure S2: Cellular movement, immune cell trafficking, immunological disease network.