Candidate Signature miRNAs from Secreted miRNAome of Human Lung Microvascular Endothelial Cells in Response to Different Oxygen Conditions: A Pilot Study

Oxygen conditions in the lung determine downstream organ functionality by setting the partial pressure of oxygen, regulating the redox homeostasis and by activating mediators in the lung that can be propagated in the blood stream. Examples for such mediators are secreted soluble or vesicle-bound molecules (proteins and nucleic acids) that can be taken up by remote target cells impacting their metabolism and signaling pathways. MicroRNAs (miRNAs) have gained significant interest as intercellular communicators, biomarkers and therapeutic targets in this context. Due to their high stability in the blood stream, they have also been attributed a role as “memory molecules” that are able to modulate gene expression upon repeated (stress) exposures. In this study, we aimed to identify and quantify released miRNAs from lung microvascular endothelial cells in response to different oxygen conditions. We combined next-generation sequencing (NGS) of secreted miRNAs and cellular mRNA sequencing with bioinformatic analyses in order to delineate molecular events on the cellular and extracellular level and their putative interdependence. We show that the identified miRNA networks have the potential to co-mediate some of the molecular events, that have been observed in the context of hypoxia, hyperoxia, intermittent hypoxia and intermittent hypoxia/hyperoxia.


Introduction
Since their discovery in 1993 in the worm Caenorhabditis elegans, the research field of miRNAs has hugely expanded.These small RNAs have been shown to regulate the majority of cellular gene expression in a nucleotide sequence-specific way [1].Due to their short length (mostly between 19 and 22 nucleotides), and the even shorter seed region, a single miRNA is capable of regulating the expression of many genes thereby orchestrating complex gene expression patterns.They are often secreted from cells and can travel in circulation over longer distances, well protected either in vesicles or bound to lipoproteins, before they can be taken up again by distal tissues to also regulate gene expression in remote organs [2,3].Transcription of miRNAs is somewhat cell-type and tissue-specific and circulating levels have been shown to have predictive potential for organ (dys)function.Samples can be obtained from a non-invasive blood collection, and miRNAs can be specifically quantified using standard analytical procedures (qPCR).Therefore, miRNAs have attracted significant interest as biomarker candidates in the context of several pathophysiological conditions, such as cancer, neurodegeneration and diseases of the heart, kidney, liver and lung [3,4].In addition, they are regarded as promising therapeutical tools in precision medicine [5].
Acute respiratory distress syndrome (ARDS) can be elicited by different triggers, such as infection, aspiration, etc., and frequently has a fatal outcome due to its heterogenous clinical manifestation and the lack of effective therapeutic interventions.Mechanical ventilation using supplemental oxygen is necessary in many cases in order to ensure proper tissue oxygenation but can, per se, have detrimental impact on the lungs, including volu-, baro-and atelectotrauma.The worst danger, however, comes from biotrauma, which is still not understood in every detail.Hyperoxia induces oxidative stress and the formation of oxygen-derived radicals, as well as impacts on normal signaling pathways, that can be regulated dependent on the partial pressure of oxygen.Impaired and mechanically ventilated lungs are further prone to forming atelectasis resulting in oxygen patterns that oscillate during the breathing cycle between hypoxic and hyperoxic levels.In animal models, it has been shown that these oxygen oscillations are transmitted via the circulation to remote organs and might have pathophysiological consequences [6].This condition of intermittent hypoxia/hyperoxia has been shown to result in the activation of biological pathways different from constant hypoxia or hyperoxia [7].Oscillations in the hypoxic range are observed in obstructive sleep apnea syndrome (OSAS), commonly designated as "intermittent hypoxia", which underlies part of the pathophysiology of the condition.Apart from reactive oxygen species (ROS) miRNAs might also be involved in regulating these effects in an auto-, para-or endocrine way.Therefore, in this study, we primarily aimed to analyze miRNAs secreted from human lung microvascular endothelial cells in a well-defined culture system in response to different oxygen conditions.We aimed to predict potential effects on nearby or remote cells by using their known target genes in enrichment analyses and tested possible autocrine effects in their mother cells.The identification of miRNAs in the lung endothelial secretome as a signature for diverse oxygen conditions could help to delineate molecular effects on target mRNAs, elucidate potential impact of oxygen conditions on lung injury in the context of pathological states such as ARDS or OSAS and could enable us to propose candidate miRNAs that could serve as biomarker (panels) originating from the lung endothelium.

Workflow
Human lung microvascular endothelial cells (HMVEC-Ls) from six supposedly healthy Donors (but unknown anamnesis) were exposed to constant and oscillating oxygen (O 2 ) conditions.Oxygen conditions were chosen to represent moderate and mild hypoxia (5% and 10% O 2 ), strong hyperoxia (95% O 2 ) and intermittent hypoxia (0-21% O 2 ) or intermittent hypoxia/hyperoxia (0-95% O 2 ) relative to normoxic culture conditions (21% O 2 = reference condition in all calculations).Supernatants were concentrated via ultrafiltration in order to enrich them for vesicle-or lipid-bound miRNAs.miRNAs in concentrates and mRNAs in lysates of mother cells were analyzed via sequencing followed by bioinformatic data analyses (Figure 1).The aim of this study was to possibly identify O 2 condition-specific changes in secreted miRNAs and their putative autocrine or para-and endocrine effects.

Characterization of Secreted miRNA during Different Oxygen Conditions
Analysis of secreted miRNAs (primarily from extracellular vesicles, but also lipidbound) was performed using clarified and concentrated conditioned medium from HMVEC-Ls.
Prior to next-generation sequencing of miRNAs, 42 miRNAs were compared via qPCR to identify a potential impact of concentration via ultrafiltration on miRNA levels.The data showed that most miRNAs were increased in concentration by a mean factor of about 12.5, indicating uniform recovery of miRNAs after ultrafiltration.Spike-in controls also indicated an absence of PCR inhibition after ultrafiltration, which was relevant for additional miRNA quantifications via qPCR (Supplemental Data S1A).
The hierarchical mapping approach of quality-and size-filtered sequencing of raw data against various transcriptome databases revealed comparable RNA composition across all samples despite differences in absolute read counts (Figure 2A). Figure 2B shows the 21 top abundant miRNAs in supernatants, and the whiskers indicate variability between the six Donors from this study.A heatmap created from data of all detected miR-NAs and principal component analysis reveal strong Donor-specific characteristics (Supplemental Data S1B).

Characterization of Secreted miRNA during Different Oxygen Conditions
Analysis of secreted miRNAs (primarily from extracellular vesicles, but also lipidbound) was performed using clarified and concentrated conditioned medium from HMVEC-Ls.
Prior to next-generation sequencing of miRNAs, 42 miRNAs were compared via qPCR to identify a potential impact of concentration via ultrafiltration on miRNA levels.The data showed that most miRNAs were increased in concentration by a mean factor of about 12.5, indicating uniform recovery of miRNAs after ultrafiltration.Spike-in controls also indicated an absence of PCR inhibition after ultrafiltration, which was relevant for additional miRNA quantifications via qPCR (Supplemental Data S1A).
The hierarchical mapping approach of quality-and size-filtered sequencing of raw data against various transcriptome databases revealed comparable RNA composition across all samples despite differences in absolute read counts (Figure 2A). Figure 2B shows the 21 top abundant miRNAs in supernatants, and the whiskers indicate variability between the six Donors from this study.A heatmap created from data of all detected miRNAs and principal component analysis reveal strong Donor-specific characteristics (Supplemental Data S1B).

Identification of Secreted miRNA Subsets in Response to Different Oxygen Conditions
Differential expression of miRNAs in different groups (with 21% O 2 as reference treatment) was investigated using edgeR.We obtained a list of miRNAs, that were changed with p-values < 0.05 relative to normoxic incubation; however, after corrections for multiple testing according to the Benjamini-Hochberg procedure, only hsa-miR-181b-5p achieved an FDR <0.1 in the 95% O 2 treatment group (Table 1).In addition to a Donor effect, which was identified using unsupervised analysis (Supplemental Data S1B-Heatmap, PCA and t-SNE [8]), treatment-induced changes in miRNA abundance usually are subtle, and a biological impact is mainly achieved through cooperative networks of several miRNAs.We therefore selected miRNAs with p < 0.05 (graphically depicted in Volcano plots, Figure 3) and extended the number of Donors for verification of some miRNAs (top 3-5 miRNAs with smallest p-value from NGS) via qRT-PCR (Table 1).across all samples despite differences in absolute read counts (Figure 2A). Figure 2B

Identification of Secreted miRNA Subsets in Response to Different Oxygen Conditions
Differential expression of miRNAs in different groups (with 21% O2 as reference treatment) was investigated using edgeR.We obtained a list of miRNAs, that were changed with p-values < 0.05 relative to normoxic incubation; however, after corrections for multiple testing according to the Benjamini-Hochberg procedure, only hsa-miR-181b-5p achieved an FDR <0.1 in the 95% O2 treatment group (Table 1).In addition to a Donor effect, which was identified using unsupervised analysis (Supplemental Data S1B-Heatmap, PCA and t-SNE [8]), treatment-induced changes in miRNA abundance usually are subtle, and a biological impact is mainly achieved through cooperative networks of several miRNAs.We therefore selected miRNAs with p < 0.05 (graphically depicted in Volcano plots, Figure 3) and extended the number of Donors for verification of some miR-NAs (top 3-5 miRNAs with smallest p-value from NGS) via qRT-PCR (Table 1).In order to identify potentially more robust miRNA hits, we used the R-package LIMMA as an additional data analysis approach.Originally developed for the detection of differentially expressed genes in microarrays, LIMMA has developed as a widely used standard in the analysis of transcriptomics and proteomics data in experiments with complex experimental designs.Briefly, LIMMA fits a regression-like linear model to the data and experiment design, thus allowing the estimation of the logFC, the p-value and the adjusted p-value.LIMMA analysis was followed by weighted gene coexpression network analysis (WGCNA) according to Langfelder and Horvath [9].Briefly, in WGCNA, coexpression networks are calculated using Pearson s correlation coefficient.These In order to identify potentially more robust miRNA hits, we used the R-package LIMMA as an additional data analysis approach.Originally developed for the detection of differentially expressed genes in microarrays, LIMMA has developed as a widely used standard in the analysis of transcriptomics and proteomics data in experiments with complex experimental designs.Briefly, LIMMA fits a regression-like linear model to the data and experiment design, thus allowing the estimation of the logFC, the p-value and the adjusted p-value.LIMMA analysis was followed by weighted gene coexpression network analysis (WGCNA) according to Langfelder and Horvath [9].Briefly, in WGCNA, coexpression networks are calculated using Pearson's correlation coefficient.These coexpression networks are segmented into groups of genes (modules) via hierarchical clustering of the coexpression matrix and cutting the resulting dendrogram using dynamic tree cutting.Thus, groups of genes with similar expression patterns are obtained, which can be put into biologic context using term enrichment analysis.WGCNA reduces the problem of post hoc correcting of data for multiple testing.Highly interconnected genes are clustered in "modules" to summarize genes with high absolute correlations.Characteristics of such networks are the following: (1) the "connectivity" (or "degree"), which is a measure for the connection strength with other network genes; (2) the "intramodular connectivity k within ", which is a measure for connectivity of a gene to the genes of a particular module; (3) the module "eigengene", which is the first principal component of a module and therefore the main representative of the gene expression profile in this module; (4) the "eigengene significance", which is the correlation coefficient of the module eigengene with the outcome (sample trait); (5) the "module membership" (correlation of the expression profile of a gene with the module eigengene of a given module); ( 6) "hub genes" representing genes inside modules with high connectivity; ( 7) "gene significance", which provides a measure of how significant a certain gene is for the biological question and ( 8) "module significance", which is the average gene significance of all genes in a given module.
In summary, the workflow comprises the construction of a gene coexpression network via correlation analysis, the identification of modules through hierarchical clustering, putting modules in relation to biological information, reduction of biological data to eigengenes and the identification of key drivers in interesting modules (intramodular connectivity) in order to identify biomarker candidates.
The WGCNA of our data from secreted miRNAs in this study resulted in 15 modules (M0-M14) with p-values not below 0.05 (Table 3; Figure 4A), which certainly limits any valid interpretation.Underlying causes are the small sample size (due to limited availability of different Donor cells with appropriate growth characteristics) and the rather small treatment effect on levels of individual secreted miRNAs.In order to deduce a possible trend that can be a starting point for further studies, we make an interpretation attempt: modules exhibiting correlations with the lowest p-value were M8 (k = 0.393; p = 0.149 for 5% O 2 ; k = −0.293;p = 0.295 for 0-21% O 2 oscillations), M13 (k = 0.358; p = 0.193 for 0-95% O 2 oscillations), M14 (k= 0.324; p = 0.243 for 95% O 2 ), M0 (k = −0.271;p = 0.334 for 0-21% O 2 ) and M2 (k = −0.251;p = 0.372 for 0-21% O 2 ).Interestingly, we found the miRNAs from the alternative arm of the NGS-identified precursor miRNAs let-7i-3p, miR-1304-5p, miR-125a-3p, miR-224-3p, miR-629-3p, miR-134-3p and miR-129-1-3p in modules M0 and M2 for 0-21% O 2 .Correlation coefficients and p-values of all modules are shown in Table 3.In order to predict putative biological processes mediated via miRNA-targeted genes, we performed enrichment analysis of miRNA targets from individual modules.Detailed results from the enrichment analysis are listed in Supplemental Data S2.Items with low p-values represent possible candidates as a basis for further mechanistic studies.
Table 4 shows connectivity values of miRNAs central in the respective modules, indicating especially important miRNAs for the processes represented by the module.The gray module combines all remaining mRNAs with no coexpression behavior.
2.4.RNA Sequencing from the Mother Cells of the Secretome RNA sequencing was primarily performed for the secreted miRNAs, but also for mRNAs of the mother cells (HMVEC-Ls) from three Donors under the same O 2 exposure condition.This was realized in order to get an estimate of gene expression regulation on the level of the mother cells.Cellular miRNA levels will be different from secreted miRNAs, but secreted miRNAs can also influence their mother cells (=autocrine effect).The autocrine effect of secreted miRNAs, however, is difficult to dissect from intrinsic miRNA-and other regulative mechanisms of mRNA expression and stability, but we attempted to deduce possible parallel tendencies.
According to data analysis of secreted miRNAs as described above, LIMMA-WGCNA was applied to the RNAseq data from cell lysates resulting in 95 modules with similar coexpression (Figure 4B).Correlations of all modules are shown in Table 5.The modules of cellular mRNAs exhibiting statistically significant positive or negative correlations with O 2 conditions are highlighted in Table 5.
Detailed data from the enrichment analysis of all modules are extensive and are therefore shown in Supplemental Data S3.Modules with strongly significant correlations (p < 0.0001) and their key enriched biological processes (GOBP terms) are summarized in Table 6.

Comparison of Cellular Signaling in Response to O 2 Oscillations with Different Amplitudes
It is especially interesting to use these data to compare molecular mechanisms between oxygen oscillations with the same frequency but different amplitudes (0-21% O 2 and 0-95% O 2 ).Outcome parameters from previous studies from our research group have indicated that molecular signaling pathways most likely are different in several aspects [7,10].The structures of modules with significant p-values (as indicated by correlation coefficients) are overlapping in some aspects but also reversed in other details (see Table 5).Generally, stronger correlations with higher statistical significance are observed in the samples exposed to 0-21% O 2 compared to 0-95% O 2 indicating a stronger expression of trait.Modules with significant positive correlation in both oscillation treatments (M47, M77 and M80) are associated with SMAD and BMP signaling, negative regulation of TGFβR pathway, cytoprotection by HMOX1, fucosylation, HIF1 targets, organization of cell-cell junction, positive regulation of exocytosis and excretion (including extracellular vesicles), glycolysis and gluconeogenesis, as well as glucose metabolism.Modules with significant negative correlation in both oscillation treatments (M14, M32, M42 and M66) are associated with GPI anchor biosynthesis, ceramide pathway, regulation of cell-cell adhesion, PTEN pathway, metabolism of porphyrins, nucleotide excision repair, double strand break repair via NHEJ, positive regulation of receptor mediated endocytosis, programmed cell death, protein targeting to peroxisomes.Positive correlation in 0-21% O 2 treatment and negative correlation in 0-95% O 2 treatment (M8, M54, M84 and M91) are associated with regulation of ROS biosynthesis, cytochrome complex assembly, degradation of extracellular matrix, TNFR2 non-canonical NF-κB pathway, mitochondrial respiratory chain complex assembly, cytokine production, calcium signaling using intracellular calcium, cell-cell communication, regulation of the nitric oxide metabolic process, positive regulation of stress fiber assembly, reactive nitrogen species (RNS) metabolic process, interferon-gamma production, IL-1, IL-6, IL-12 and TNFα (generally cytokine) production, positive regulation of TLR and pattern-recognition receptor signaling pathway and collagen degradation.Negative correlation in 0-21% O 2 treatment and positive correlation in 0-95% O 2 treatment (M7, M22, M23, M56 and M72) are associated with morphogenesis (tube formation), cell division and mitotic cell cycle, positive regulation of TORC1 and TOR signaling, positive regulation of response to DNA damage stimulus, negative regulation of extrinsic apoptotic signaling pathways, ribosome and translation initiation, response to starvation, the PDGF pathway, telomere maintenance via telomere lengthening, cellular senescence, lipid modification, fatty acid metabolism and protein localization to the ER.

Putative Autocrine Effects of Secreted miRNAs on Cellular mRNA Targets
The in vitro system of this study uses only one cell type (HMVEC-L); therefore, target cell effects of secreted miRNAs can only be directly referred to the mother cells (=autocrine effect).In order to relate changes in miRNA levels in the secretome to alternated mRNA levels of putative miRNA targets in the cells, we applied iTALK, an open source R package that can be used to characterize and illustrate intercellular communication: https://github.com/Coolgenome/iTALK (Accessed on 7 March 2024) [11,12].The software was designed to visualize ligand-receptor-mediated intercellular cross-talk but can also be adapted to miRNA-mRNA interactions from sequencing data.
Output data were filtered and ranked according to q-value and interaction pairs showing an inverse expression trend between miRNA and putative target mRNA were visualized in Circos plots.Clear trends were only observed with treatment groups of 0-21% O 2 and 95% O 2 (Figure 5).Identified target genes were analyzed with the Reactome database (www.reactome.org,Accessed on 7 March 2024) in order to deduce potential autocrine effects (Table 7).

95% O2
Figure 5. Visualization of iTALK outputs as Circos plots.miRNAs with altered abundance in the secretome after 48 h treatment with 0-21% O2 and 95% O2 are plotted and put into relation with putative target genes from cellular RNAseq data with purple arrows indicating downregulated miRNAs and upregulated target genes and beige arrows linking upregulated miRNAs and the respective downregulated mRNA targets.

Putative Para-and Endocrine Effects
In order to visualize differences of secreted miRNA/target mRNA interactions under O 2 oscillations with different amplitudes, we constructed secreted miRNA-putative target mRNA networks using www.mirnet.casoftware miRNet 2.0 for 0-21% O 2 and 0-95% O 2 oscillations (Figure 6).Enrichment analysis of target genes in such networks for constant and intermittent O2 conditions was performed in order to predict possible biological effects that might be induced in other target cells via cooperative action of their associated miRNAs (Table 8).Enrichment analysis of target genes in such networks for constant and intermittent O 2 conditions was performed in order to predict possible biological effects that might be induced in other target cells via cooperative action of their associated miRNAs (Table 8).

Discussion
The original purpose of our study was to aim for the identification of marker miRNAs that are secreted from the human pulmonary microvasculature in response to different oxygen conditions or patterns and play a role in the clinical setting, such as different hypoxic conditions as in tumors or due to O 2 undersupply, hyperoxia elicited by supplemental O 2 therapy, intermittent hypoxia as in OSAS or also in tumors and intermittent hypoxia/hyperoxia due to atelectasis in ventilated patients with higher FiO 2 .Influencing factors in this context are highly complex: Donor characteristics such as age, gender, ethnicity and confounding factors (co-morbidities, smoker and environmental factors), cell type, duration of gas exposure and timepoint of sampling and oxygen concentration, just to name a few.The task becomes even more difficult if we want to conclude downstream effects of such specific signature miRNAs, as they will have multiple potential target mR-NAs engaging in a plethora of sometimes even contrasting signaling pathways.miRNAs secreted into the circulation may be taken up by cells of nearby or remote tissues, and it is still unclear if these uptake mechanisms are highly specific or partly promiscuous.Bioinformatic analysis methods help to sort extensive data, portray an overall picture of the setting and can contribute to narrowing the hypothesis.
In our pilot study, we tried to simplify the system as much as possible.We chose as a model system pure preparations of supposedly healthy human pulmonary microvascular endothelial cells in culture from Donors of different gender and ethnicity, in order to identify treatment-specific effects independent from the individual.These Donor cells were exposed in parallel to the different O 2 conditions for 48 h.Exposure duration was chosen to achieve the best trade-off between the amount of secreted miRNAs (protein-bound or in extracellular vesicles-EVs) and stability of EVs.
Pure preparations of HMVEC-Ls from healthy human lungs of sufficient cell counts to conduct these experiments represent very precious material and availability is very limited.It is generally difficult to obtain effects with statistical significance after post hoc corrections (Benjamini-Hochberg) when testing for hundreds of individual items (=miRNAs in NGS), especially because several miRNAs usually show small to moderate changes upon treatment but act in cooperation to alter downstream mRNA expression resulting in biological effects.
As part of a pilot study, our results were extracted from the NGS analysis of six Donors and a fraction of miRNA alterations could be verified in a further six Donors via qRT-PCR.Data analysis and prediction of potential functional implications were performed with the help of bioinformatic methods.Most identified extracellular miRNAs have been annotated to occur in extracellular vesicles and some miRNAs have been already described in the literature in the context of endothelial function (for example: miR-92a, miR-126, miR-98, and miR-181).Many miRNAs have been mainly investigated in the context of cancers; therefore, the role of these miRNAs in non-cancerous cells needs to be interpreted with caution.
The exact experimental condition is crucial for the interpretation of the outcome.The literature contains a large number of different studies that have investigated, for instance, hypoxia-regulated miRNAs (HRM).Their significance, however, is strictly dependent on the cell type, the actual O 2 concentration and the duration of exposure.Central regulators of the hypoxic cell response are a family of hypoxia-inducible transcription factors (HIFs).Part of HIF-dependent signaling is regulated by various miRNAs, such as HIF biogenesis, activation, degradation and isotype-switch and translation of HIF-dependent proteins [13].Moreover, 5% O 2 is relatively moderate but not severe hypoxia and might be above the threshold of induction for hypoxia-induced transcription factors of the HIF family.HIF1A is even a target gene of the upregulated miR-19b, as well as the DNA-methyltransferase DNMT3A, which regulates DNA methylation under hypoxia.miR-19b together with miR-23a has been identified as a hub miRNA linking hypoxia and DNA methylation in thrombotic conditions [14].We identified miR-19b-3p and miR-23a-3p as miRNAs with high intramodular connectivity (k within ) in M1 in WGCNA, which also reveals their function as hub miRNAs in this module.miR-19b also targets PTEN, a negative regulator of the Akt-pathway and thereby enhances proliferation and inhibits apoptosis, which is a common observation under mild hypoxia.Similarly, miR-21 (in M3) targets PTEN and impacts DNA methylation [15].Downregulation of the tumor suppressor miR-26a and let-7 family member let-7d-3p will further support cell proliferation [16].miR-181c has been associated with COPD [17].miR-1307 is an oncogenic miRNA associated with proliferation and metastasis and has been shown to play an important role in cellular infection by SARS CoV-2 virus and disease progression [18].miR-92a belongs to the miR-17-92 cluster and is used as exosomal biomarker in several cancers.Levels of miR-92a in exosomes can indicate the level of activation of the endothelium and endothelial dysfunction.In cardiovascular medicine it is used as a marker for the distinction of acute myocardial infarction and stable coronary heart disease [19].Pathway analysis of these extracellular miRNAs revealed their active function in a typical and previously described hypoxic response.Enriched pathways of miRNAs secreted under 5% O 2 (= moderate hypoxia) included the following: TGFß and VEGF signaling as a typical hypoxic response [20], ATF4 induces VEGF-A and accumulates in ER under hypoxia [21], increased cell proliferation supported by decreased cellular senescence and DNA damage/telomer stress, a tendency to pro-thrombotic events via fibrin clot formation [22], decreased collagen degradation in the ECM [23], smooth muscle contraction (possibly as part of the pulmonary hypoxic vasoconstriction), hypoxic p53 signaling [24], calcium-activated K+ channels [25], glucose transporter translocation [26] and regulation of lipid metabolism [27].Other miRNAs that are frequently named in the context of hypoxic regulation, such as miR-210, miR20a/b, miR-199a, miR-200b/c, miR-424 and miR-429, did not stand out in this context, either because they strongly varied between Donors (i.e., hsa-miR-424 and hsa-miR-199a) or reads were extremely low or zero in the secretome (i.e., hsa-miR-20a/b, hsa-miR-210, hsa-miR-200b/c and hsa-miR-429).
Furthermore, 10% O 2 is an even milder hypoxic condition and revealed only a downward trend of miR-4492, possibly reducing cellular senescence.
Additionally, 95% O 2 is a very strong hyperoxic condition that induces cell stress, DNA damage by ROS, inflammation via NF-KB activation and senescence.ROS mainly occurs in the mitochondria and through NADPH oxidase (NOX).The miRNA-181a and 181b couple target several proteins involved in mitochondria biogenesis, turnover and clearance as well as detoxification of ROS.Downregulation of miR-181a and 181b has proven to ameliorate disease conditions associated with mitochondrial dysfunction, such as retinal diseases, where oxidative stress is an issue [28].miR-155 has a critical role in oncogenesis, inflammation and immunity [29].
The 0-21% O 2 condition, or intermittent hypoxia, is a model for oxygen conditions occurring in tumors or in obstructive sleep apnea (OSAS) and has been shown in the literature to activate very specific signaling pathways.miR-30b-5p has been shown to play an important role in the suppression of lysosmal biogenesis and autophagy [30].It is a tumor suppressor in the context of lung cancer and enhances apoptosis [31].Experimentally validated targets are CCNE1, SOCS1, SMAD1, CAT and BCL6 indicating a role in the regulation of cytokine signaling, apoptosis and oxidative stress.Validated targets of miR-129-5p are CAMTA1 (calmodulin-binding transcription factor), indicating a role of calcium under this condition, SOX4 (regulation of Wnt pathway), GALNT1 (modifying glycan structures), EIF2C3 and HMGB1, which has been shown to protect against ischemia/reperfusion injury [32].Apart from the enriched pathways that have miR-125a-5p in common with other miRNAs, an interesting finding comes from the literature.Hirota et al. showed that miR-125a-5p and its family member miR-125b-5p target ATRAP (AT1R-associated protein), thereby activating angiotensin II-AT1R signaling and facilitating hypertension, which is also a phenotype in OSAS [33].miR-654-3p was found to reduce inflammation by targeting ADAM10 and RAB22A (both promoters of atherosclerosis) [34] and TNFRSF9 (TNF receptor family member 9) [35].Activated NF-KB can upregulate miR-224-5p [36], which was shown to be a circulating indicator of patients' responsiveness to glucocorticoid therapy as it downregulates glucocorticoid receptor degradation by inhibiting GSK-3ß [37].miR-134-5p has been shown to play a role in the regulation of vascular smooth muscle cell phenotypic switch, vascular development, platelet activation, regulation of the excretion of matrix metalloproteinases of the ADAMTS family, inflammation and calcium signaling.Overexpression seems to suppress pathological vascular remodeling [38].Downregulated miRs under 0-21% O 2 were miR-410-3p, let-7i-5p, miR-92a-3p, miR-4497 and miR-629-5p.Downregulation of miR-410-3p has been shown to enhance adipogenesis mediated via its target IRS-1 (insulin receptor substrate-1) [39].Intermittent hypoxia is a condition of enhanced oxidative stress.Li et al. have shown that exosomes released from endothelial cells under oxidative stress reveal downregulated miR-92a-3p, which might stimulate angiogenesis via its target gene Tissue factor [40].In a previous study, our group was able to show that Tissue factor expression was enhanced under intermittent hypoxia and was increasingly released within extracellular vesicles [41].Direct targets of miR-4497 are NF-KB2 and SP1 transcription factor.Therefore, down regulated miR-4497 might enhance inflammation, apoptosis and the DNA damage response [42].miR-629-5p has been shown to be upregulated in the plasma of patients with pulmonary arterial hypertension and also in hypoxic pulmonary artery smooth muscle cells, promoting vascular remodeling via FOXO3 and PERP [43].Intermittent hypoxia seems to induce downregulation of miR-629-5p.
Hypoxic/hyperoxic oscillations (0-95% O 2 ) have been shown in previous studies to use different signaling pathways than hypoxic oscillations or constant hypoxic or hyperoxic conditions [7,44].We found, as exosomal miRNAs upregulated under this condition, miR-574-3p, which has been shown to be linked to HIF/VEGF signaling at least under hypoxia [45], and miR-23a-3p (upregulated as under 5% O 2 ).miR-155-5p was similarly downregulated as under constant strong hyperoxia (95% O 2 ).Let-7d-3p was downregulated analogous to hypoxic (5% O 2 ) treatment.An interesting target of identified miR-4492 is CAPON (nitric oxide synthase 1 adapter protein), which interacts with NOS1 and plays a role in psychiatric as well as heart diseases.Downregulated miR-4492 might increase CAPON expression, which is linked to induction of cardiac disease (sudden cardiac death) [46].CAPON also plays a role in diabetes.It has a functional role in NOS1 anchoring and activation.In a previous study from our group, we found increased levels of nitric oxide-derived peroxynitrite in pulmonary endothelial cells exposed to hypoxic/hyperoxic oscillations [7].Whether this is due to altered NO production and/or oxidative stress has not been delineated yet.miR-301b-3p plays a role in different types of cancers.Functional enrichment of oxygen-dependent proline hydroxylation points to a role of HIF regulation, similar as with miR-130a-3p.Downregulation of miR-130a-3p seems to have a pro-inflammatory effect on endothelial cells [47].miR-7-5p targets are experimentally well validated and include PAK1, EGFR, IRS, RAF, CKAP and others.Functional studies of altered expression have mostly been conducted in the context of tumors.miR-98-5p levels have been found to be decreased in the serum of stroke patients and mice with ischemic conditions.Yu et al. found that miR-98-5p protects against oxidative stress (possibly also in the context of I/R injury) and inhibits apoptosis [48].
Analysis of mRNA expression in secreting mother cells was limited to a still smaller sample size (three Donors) but revealed molecular characteristics that are known to be associated with the respective O 2 condition.An attempt was made to correlate changes in mRNA levels as targets of altered secreted miRNAs (= a possible autocrine effect).In samples from 0-21% O 2 and 95% O 2 , we could identify target mRNAs of altered miRNAs in the secretome that were inversely regulated.A true autocrine effect cannot be discerned from independent effects of cellular miRNAs regulating expression of cellular mRNAs.However, a recent study shows that endothelial cells are capable of reprograming vascular cells via their EV-cargo (miRNAs and proteins) which can have different apical and basolateral loading mechanisms [49].
Again, a major limitation of our study is the small number of different Donor cells due to a very limited availability of healthy human lung tissue.This fact restricts the statistical power of this study, and, additionally, relatively small treatment effects impede the identification of miRNAs that could serve as O 2 -condition-specific markers from the lung microvascular endothelium.However, some identified miRNAs can be regarded as candidates and should be validated in further studies in the future.

Cell Culture
Primary human lung microvascular endothelial cells were purchased from LONZA (HMVEC-L cat.nr.CC-2527) and were propagated in EGMTM-2MV microvascular endothelial cell growth medium-2 including growth factor supplements (cat.nr.CC-3202) (LONZA, Basel Switzerland).For details on Donor characteristics (age, gender and ethnicity) see Supplemental Data S5.For the experiments, cells were plated into gelatine-coated 6-well cell culture plates with a gas-permeable fluorocarbon membrane as growth support ("Imaging Plates", Zellkontakt, Nörten-Hardenberg, Germany) using full growth medium.Prior to gas exposure, the medium was replaced with 2 mL/well serum-free EGMTM-2MV medium.

Processing of Cell Culture Supernatants
After 48 h, gas exposure cell culture supernatants were collected and were centrifuged for 10 min at 200 g to remove cells and larger debris and for 10 min at 2000 g to remove smaller parts such as apoptotic bodies.Enrichment of extracellular vesicles was achieved via ultrafiltration using Ultra-15 filters (30 kDa, MERCK cat.nr.UFC903024, miRNet 2.0) at 3200 g.Supernatants were concentrated from 12 mL to a final volume of 500 µL.Then, 200 µL concentrate was mixed with 1 mL Qiazol lysis reagent (Qiagen, Hilden, Germany) and left at room temperature for 10 min and then frozen at −70 • C. Additionally, 300 µL were directly frozen at −70 • C in low-binding tubes without prior lysis (Eppendorf, Hamburg, Germany) until further processing.

RNA Isolation for Next-Generation Sequencing
Total RNA extraction was performed from 200 µL ultrafiltrate using the miRNeasy mini kit (Qiagen, Hilden, Germany).Synthetic oligonucleotides obtained from the miR-CURY Spike-In kit (Qiagen, Hilden, Germany) were added to the Qiazol lysis buffer before homogenization.Glycogen (5 mg/mL) was added to the chloroform extract at a 1:100 dilution to enhance precipitation.All other steps were performed according to the recommendations of the manufacturer.Total RNA was eluted in 30 µL nuclease-free water and stored at −80 • C in low-binding tubes until further analysis.

Small RNA-Sequencing Analysis (Next-Generation Sequencing)
Equal volumes of total RNA were used for small RNA library preparation using the CleanTag library preparation kit (TriLink Biotechnologies, San Diego, CA, USA).Adapterligated libraries are amplified using barcoded Illumina reverse primers in combination with the Illumina forward primer.Pools of 16-24 smallRNA sequencing libraries were prepared at equimolar rates on the basis of DNA-1000 high-sensitivity bioanalyzer results (Agilent, Santa Clara, CA, USA).Sequencing was performed on an Illumina HiSeq 2500 System with 50 bp single-end runs.Sequencing Reads were adapter trimmed and filtered for low-quality reads (Q < 30).MicroRNA annotation was performed using the miND pipeline [51,52].

RNA Sequencing (RNAseq) of Cellular mRNAs
Total RNA was isolated from attached cells after removing supernatants using the RNeasy mini plus kit (Qiagen, Hilden, Germany).Sequencing libraries for total RNA were prepared at the Core Facility Genomics, Medical University of Vienna, using the QuantSeq 3 ′ FWD protocol version 2 with unique dual indices (Lexogen, Vienna BioCenter), following the low input branch of the protocol.A total of 20 PCR cycles were used for library prep, as determined via qPCR according to the library prep manual.Libraries were quality-checked on a Bioanalyzer 2100 (Agilent) using a High-Sensitivity DNA Kit for correct insert size and quantified using Qubit dsDNA HS Assay (Invitrogen).Pooled libraries were sequenced on a NextSeq500 instrument (Illumina) in 1x75bp single-end sequencing mode.An average of 6 million reads per sample were generated.Reads in fastq format were generated using the Illumina bcl2fastq command line tool (v2.19.1.403).Reads were trimmed and filtered using cutadapt version 2.8 to trim polyA tails, remove reads with N's and trim bases with a quality of less than 30 from the 3 ′ ends of the reads.On average, 4.8 million reads were left after this procedure.Trimmed reads in fastq format were aligned to the human reference genome version GRCh38 with Gncode 29 annotations using STAR aligner [53] version 2.6.1a in 2-pass mode.Raw reads per gene were counted with STAR.4.9.Bioinformatical Processing of Data 4.9.1.Extracellular miRNAs NGS sample data were processed according to the miRNA next-generation-sequencing discovery assay (miND) [51].Differential expression of miRNAs in the endothelial secretome was analyzed by the established toolkit "edgeR" from NGS data [54].Identified miRNAs found to be significantly changed in their abundance in the cell culture supernatants in response to different oxygen conditions were uploaded to the web-based platform miRNet 2.0 (https://www.mirnet.ca) in order to construct miRNA-target mRNA networks and to perform pathway enrichment analysis using Reactome database.

Cellular mRNAs
Raw data of cellular mRNA gene expression as obtained from RNAseq analysis were analyzed using the R-packge "LIMMA" followed by gene set variation analysis (GSVA) [55].A coexpression network was generated using weighted gene coexpression network analysis (WGCNA) [9] (see below), and modules of correlated genes with similar coexpression patterns were identified.Pathway enrichment analysis was performed on genes belonging to individual modules as well as correlation analysis in order to identify modules which were, in a statistical sense, significantly outstanding with a positive or negative correlation coefficient compared to all other treatment groups.
WGCNA was performed using variance-stabilized transformed counts matrices and the WGCNA package of R using signed networks, Pearson correlation and a beta of 8.The biologic context of the resulting modules was assessed using the ClusterProfiler package [56] and the molecular signature database, collection Geneontology-biological process (GoBP) [57].
miRNA-RNA interactions were determined using the iTALK package of R and the molecular signature database (collection C2) to determine mRNA targets of miRNAs [11,57].

Conclusions
This pilot study let us identify some candidate miRNAs that have the potential to serve as indicators of the (lung) microvascular endothelium response to altered oxygen conditions and might mediate some of the downstream effects of these conditions on circulation and other organs.Owing to the rare availability of pure pulmonary microvascular endothelial cell preparations from healthy human lungs, the rather small sample size of this study does not allow for the unambiguous identification of miRNAs that function as biomarkers for different oxygen conditions.However, in the light of other data from previous studies, the predicted biological processes affected by these miRNAs fit into the picture.These results can and should be validated in future preclinical or clinical studies in order to better evaluate their potential as part of biomarker panels.

Figure 1 .
Figure 1.Workflow showing HMVEC-L exposure to different O2 conditions, biochemical processing, RNA sequencing and bioinformatic processing of data.

Figure 1 .
Figure 1.Workflow showing HMVEC-L exposure to different O 2 conditions, biochemical processing, RNA sequencing and bioinformatic processing of data.
show the 21 top abundant miRNAs in supernatants, and the whiskers indicate variability b tween the six Donors from this study.A heatmap created from data of all detected mi NAs and principal component analysis reveal strong Donor-specific characteristics (Su plemental Data S1B).

2 Figure 2 .
Figure 2. Visualization of data from next-generation sequencing of secreted extracellular (mi)RNAs.(A) Quantitative distribution of different types of RNAs detected in individual samples as obtained from raw reads without any further normalization.(B) Top 21 miRNAs with highest abundance (shown as mean and range from 6 Donors in "Reads per Million"-RPM).

O 2 Figure 3 .
Figure 3. Volcano plot showing log (fold change) of miRNA abundance relative to uncorrected pvalue.Shown are all miRNAs with 0.05 < p-value < 1.0 as black dots.The dotted line represents a threshold of p = 0.05.miRNAs with p < 0.05 are labeled and assigned blue (downregulated) or red (upregulated) symbols.

Figure 3 .
Figure 3. Volcano plot showing log (fold change) of miRNA abundance relative to uncorrected p-value.Shown are all miRNAs with 0.05 < p-value < 1.0 as black dots.The dotted line represents a threshold of p = 0.05.miRNAs with p < 0.05 are labeled and assigned blue (downregulated) or red (upregulated) symbols.

Figure 4 .
Figure 4. (A) WGCNA of miRNA identified and quantified via NGS in the HMVEC-L secretome.Clustering tree (dendrogram)generated through hierarchical clustering of adjacency-based dissimilarity detected modules with similar coexpression.The dendrogram and module distribution shows that the

Figure 5 .
Figure 5. Visualization of iTALK outputs as Circos plots.miRNAs with altered abundance in the secretome after 48 h treatment with 0-21% O 2 and 95% O 2 are plotted and put into relation with putative target genes from cellular RNAseq data with purple arrows indicating downregulated miRNAs and upregulated target genes and beige arrows linking upregulated miRNAs and the respective downregulated mRNA targets.

4. 8 .
Data Analysis of qRT-PCR Results and Statistical Evaluation qPCR data were analyzed by normalizing Ct values towards spike-in controls, and as quantitative measure for relative expression, we used dCt = Ct (different O 2 conditions)-Ct (21% O 2 ), which is equivalent to −log(2) fold change (FC).

Table 1 .
Fold change of secreted miRNAs with p < 0.05.as identified and quantified with NGS.Some of the data were verified with qRT-PCR.

Table 2 .
Pathways associated with individual identified miRNAs and their targets as listed in the Reactome database.

Table 3 .
Correlation coefficients of individual modules from quantified miRNAs after WGCNA.

Table 4 .
Module assignment of NGS-identified miRNAs and network connectivity.

Table 6 .
Enriched functional terms of modules with highly correlated gene expression.

Table 7 .
Pathways related to target genes conversely altered in HMVEC-Ls and secreted miRNAs (iTALK).

Table 8 .
Putative signaling pathways affected by cooperative action of miRNAs (Reactome database).