Circulating MicroRNA: Incident Asthma Prediction and Vitamin D Effect Modification

Of children with recurrent wheezing in early childhood, approximately half go on to develop asthma. MicroRNAs have been described as excellent non-invasive biomarkers due to their prognostic utility. We hypothesized that circulating microRNAs can predict incident asthma and that that prediction might be modified by vitamin D. We selected 75 participants with recurrent wheezing at 3 years old from the Vitamin D Antenatal Asthma Reduction Trial (VDAART). Plasma samples were collected at age 3 and sequenced for small RNA-Seq. The read counts were normalized and filtered by depth and coverage. Logistic regression was employed to associate miRNAs at age 3 with asthma status at age 5. While the overall effect of miRNA on asthma occurrence was weak, we identified 38 miRNAs with a significant interaction effect with vitamin D and 32 miRNAs with a significant main effect in the high vitamin D treatment group in VDAART. We validated the VDAART results in Project Viva for both the main effect and interaction effect. Meta-analysis was performed on both cohorts to obtain the combined effect and a logistic regression model was used to predict incident asthma at age 7 in Project Viva. Of the 23 overlapped miRNAs in the stratified and interaction analysis above, 9 miRNAs were replicated in Project Viva with strong effect size and remained in the meta-analysis of the two populations. The target genes of the 9 miRNAs were enriched for asthma-related Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathways. Using logistic regression, microRNA hsa-miR-574-5p had a good prognostic ability for incident asthma prognosis with an area under the receiver operating characteristic (AUROC) of 0.83. In conclusion, miRNAs appear to be good biomarkers of incident asthma, but only when vitamin D level is considered.


Introduction
Asthma is a chronic inflammatory disease characterized by airway hyper-responsiveness, mucus hypersecretion, and airway remodeling [1]. According to the Centers for Disease Control and Prevention (CDC) report in 2018, about 25 million Americans have asthma, including 7.5% of children (age < 18 years) and 7.4% of adults (age 18+ y) [2]. Asthma is more common in children than adults and the number of cases with childhood asthma increased from 2014 to 2018 by 744,397 in the US [2]. The cost of asthma in direct and indirect medical care is estimated to exceed 18 billion dollars each year [3].
Recurrent wheeze in early childhood is a significant risk factor for asthma, but not all children who wheeze go on to develop asthma. According to epidemiologic studies, 60% of young wheezers are symptom-free at the age of 6 y, and a majority of them remain asymptomatic at the age of 11 and 16 [4][5][6]. Although wheezing in early childhood is common, finding a way to predict who will continue to wheeze and develop asthma is important for asthma diagnosis and management.
Vitamin D is a fat-soluble secosteroid that can be produced in response to the skin being exposed to sunlight but is present in the diet. Vitamin D is essential for bone growth and remodeling and is also involved in immune system function [7,8]. Currently, epidemiologic and experimental evidence has suggested an association between higher vitamin D levels in pregnant women and reduced asthma in their offspring, but the results of these observational studies differ in the results [9]. To address this problem, the Vitamin D Antenatal Asthma Reduction Trial (VDAART) was developed, as a randomized, doubleblind, placebo-controlled trial with the primary outcome to determine whether prenatal vitamin D supplementation in pregnant women can decrease asthma incidence in the women's offspring [10].
MiRNAs are small non-coding RNAs with a length of 21-24 bp and they can be loaded into RISC (RNA-induced Silencing Complex) to bind to the 5 -UTR of mRNAs and thus down-regulate gene expression [11]. MiRNAs also exist in a variety of biofluids (extracellular RNAs), such as plasma, serum, sputum, and urine [12]. With the development of next-generation sequencing, an increasing number of studies have profiled small RNA-Seq for miRNA expression [13,14] and miRNAs have been utilized as the non-invasive biomarkers of both diseases and treatment responses [13,[15][16][17][18]. Some miRNAs have been identified to associate with childhood asthma in several studies [19][20][21][22].
In this study, we sought to identify circulating miRNAs that can predict incident childhood asthma and assess that prediction in relation to vitamin D treatment.

Participant Selection
We selected 75 participants from the VDAART (ClinicalTrials.gov Identifier: NCT00920621) based on the presence of recurrent wheezing occurring at 3 years of age. These participants provided plasma samples at age 3 and were followed up at 5 years of age when their asthma status was diagnosed by physicians. Two participants were not followed up at age 5 years and were removed in the following analysis.

Small RNA Sequencing
Total RNA was isolated from samples by the Qiagen miRNAeasy Serum/Plasma extraction kit and QIAcube automation. Small RNA sequencing libraries were prepared using the Norgen Biotek Small RNA Library Prep Kit and then sequenced on the Illumina NextSeq 500 platform at 51 bp single-end reads.
We used the exRNA-processing toolkit (exceRpt) to assess the read quality and profile the miRNAs [23]. Raw read counts were log-transformed and quantile normalized. Read counts less than 5 and miRNAs with coverage less than 80% of all samples were removed.

Statistical Analysis
Statistical analysis was performed using R software version 3.6.3 (Bell Laboratories, Murray Hill, NJ, USA). Logistic regression was used to assess the association between the miRNAs normalized count and asthma status (asthma: Y = 1 and no asthma: Y = 0) at age 5 years for participants in VDAART. We first examined the association in all candidate participants and then in the interaction and stratified analyses by the vitamin D treatment status (high vitamin D treatment group (4400 international units (IU)/day): TG = 1 and low vitamin D treatment group (400 IU/day): TG = 0).

Validation
We selected 20 participants in Project Viva (ClinicalTrials.gov Identifier: NCT02820402) a pre-birth cohort established to examine the association between prenatal diet and other factors with maternal and child health [24]. All of these 20 participants had a history of recurrent wheeze at 3 years old when the blood draw was taken. These plasma samples were prepared and sequenced in the same way as the VDAART samples, as well as having the miRNA profiling performed.
The vitamin D concentration was measured by both an automated chemiluminescence immunoassay and a manual radioimmunoassay at 16-26 weeks' gestation from each mother. We used the average of the two values to estimate the maternal 25(OH)D concentration [25] and divided mother-child pairs into adequate (maternal levels ≥ 20 ng/mL) and inadequate (maternal levels < 20 ng/mL) groups according to the report from the Institute of Medicine (IOM) [26]. For the 20 mother-child pairs in Project Viva, 12 children were from the adequate group and the other 8 children were from the inadequate group of maternal 25(OH)D concentration. The primary outcome was the asthma status at 7 years of age. We validated our results in the maternal 25(OH)D adequate group as well as considering the interaction by maternal vitamin D status.

Meta-Analysis
To combine the results from both VDAART and Project Viva, we performed a metaanalysis of the two data sets. The R package "metafor" was employed to combine the odds ratio (OR) values via the restricted maximum-likelihood estimator [27].

MiRNA-Target Gene Network and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Enrichment Analysis
We built the miRNA-target genes interactions network and performed Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis through miRNet 2.0 (McGill University, Montreal, Quebec, Canada) [28]. MiRNA-target interactions (MTIs) came from miRTarBase release 8.0 [29]. Nominal p values of pathways were adjusted for false discovery rate (FDR).

Prediction
We built a logistic regression model to predict asthma status at age 7 years in the Project Viva adequate vitamin D group. The area under the receiver operating characteristic (AUROC) was calculated to measure the performance of our logistic regression model in the prediction of asthma status.

Baseline Characteristics
A total of 75 participants in VDAART were selected, of whom 31 were children with asthma at age 5 years, 42 were healthy controls and 2 were lost follow-up. The baseline characteristics of the 73 subjects with available outcome for asthma status at age 5 years were shown in Table 1. Of the 73 participants, 33 were in the treatment arm of the VDAART trial where the children were born to mothers who had received 4400 IU of vitamin D3 per day during pregnancy and 40 were in the placebo arm (control group) where the children were born to mothers who had received regular multivitamins containing 400 IU of vitamin D3 per day. Sex, race, and other characteristics were evaluated and found not to be different between treatment arms.

Main Analysis
We investigated the association between miRNAs at age 3 years and asthma status at age 5 years in all 73 VDAART participants and then examined the results in Project Viva. The results are shown in Supplemental Table S1. Six miRNAs were nominally significant (p < 0.05) in VDAART and three were present in Project Viva. Only miRNA hsa-miR-548k could be validated as being in the same direction of effect in VDAART and Project Viva. Given the marginal main effects, we performed interaction and stratified analyses on our dataset by intervention arm. Data presented as n (%) or mean ± standard deviation (SD); IU: international unit; * p value from Fisher's exact test.

Significant miRNAs in Treatment Interactions Analysis
We identified 38 miRNAs that were associated with asthma at age 5 years in the interaction analysis (nominal p < 0.05) in VDAART. The results were ranked by effect size (OR) and shown in Table 2. The effect modification by vitamin D status was seen in both the positive and negative directions, with strong effect estimates noted. Sixteen miRNAs had OR values greater than 3 and five miRNAs had OR less than 0.33. Hsa-miR-3942-5p had the largest OR value of 178.73 [95% confidence interval (CI), 4.6 to 6947.95], a high risk of incident asthma and hsa-miR151a-5p had the smallest OR value of 0.05 [95% CI, 0.01 to 0.42], a strong protective effect of incident asthma.
Compared with the 32 significant miRNAs in the high vitamin D treatment group, there were only 6 significant miRNAs in the low vitamin D treatment group (control group), of which hsa-miR-505-3p and hsa-miR-340-3p had a higher risk and the others having a protective effect. Four miRNAs, including hsa-miR-7-1-3p, hsa-miR-505-3p, hsa-miR-193b-5p and hsa-miR-5010-5p, overlapped with the significant miRNAs in the interaction analysis. The results are shown in Supplemental Table S2.

Validation in Project Viva
We performed a replication analysis in Project Viva on the 23 significant miRNAs identified in both the high vitamin D treatment group and the interaction analysis in VDAART. Because of the small sample size in Project Viva, we focused our validation efforts on both the directionality and effect size. We kept miRNAs with OR ≥ 2 or OR ≤ 0.5 in VDAART and kept miRNAs with OR ≥ 1.5 or OR ≤ 0.67 in Project Viva. Following this methodology, 9 miRNAs were successfully validated in Project Viva and the results are shown in Supplemental Table S3.
Of the 9 validated miRNAs, 6 miRNAs increased the risk of incident asthma and 3 miRNAs had a protective effect on decreasing the risk of incident asthma. For example, hsa-miR-574-5p (Supplemental Figure S1A) Figure S1B) was associated with a protective effect for incident asthma in the high vitamin D treatment group in VDAART with an OR value of 0. 15

Meta-Analysis of Validated miRNAs in Vitamin D Antenatal Asthma Reduction Trial (VDAART) and Project Viva
To combine the effect in the high vitamin D treatment group of VDAART and high maternal vitamin D group (adequate group) of Project Viva, we performed a meta-analysis on the 9 validated miRNAs. The results are illustrated in Figure 1. All the observed effects of 9 miRNAs were significant in the analysis (nominal p < 0.05). Hsa-miR-574-5p had the largest risk for incident asthma with OR = 8.45 [95% CI, 1.52 to 46.97] and hsa-miR-151a-5p had the largest protective effect with OR = 0.19 [95% CI, 0.04 to 0.82]. We also performed a meta-analysis on the effect modification and the results were displayed in Supplemental Figure S2. Both hsa-miR-574-5p and hsa-miR-151a-5p had the largest effects in opposite directions.

MiRNA-Target Gene Network and KEGG Pathway Enrichment Analysis
The miRNA-target gene network was built and displayed in miRNet 2.0. There were 42 KEGG pathways enriched to the 9 validated miRNAs with FDR less than 0.05 and these pathways are shown in Supplemental Table S4. Eight signaling pathways were curated and emphasized in Table 4 with known studies about asthma. The top three signaling pathways were Cell Cycle (FDR = 4.8 × 10 −7 ), TGF-beta (Transforming Growth Factor beta) Signaling Pathway (FDR = 4.12 × 10 −5 ) and ErbB (also known as Epidermal Growth Factor Receptor) Signaling Pathway (FDR = 1.88 × 10 −4 ). The target genes involved in these signaling pathways were visualized in Figure 2.

Prediction of Asthma Based on miRNAs
We predicted asthma status based on normalized miRNA counts in the high maternal vitamin D group in Project Viva through a logistic regression model. When using a single miRNA, hsa-miR-574-5p had the best performance with AUROC = 0.83 ( Figure 3A). When considering the combination of two miRNAs, hsa-miR-215-5p (AUROC = 0.81) and hsa-miR-370-3p (AUROC = 0.75) performed best with AUROC = 0.86 ( Figure 3B).

Discussion
In our study, we identified several miRNAs that were associated with incident asthma at age 5 years given vitamin D effect modification in VDAART; these miRNAs were replicated in another independent birth cohort, Project Viva. This finding suggests that miRNA may be a pivotal vitamin D related mediator of asthma risk in early childhood. The meta-analysis of the two cohorts strengthened our results and the top miRNA had excellent prognostic power to predict asthma at age 5 years.
We noted that hsa-miR-574-5p displayed the strongest risk effect in the high vitamin D treatment group of all the validated miRNAs (Figure 1). Sinha and colleagues found hsa-miR-574-5p over expressed by 3.3-fold (p = 0.04) in asthmatic patients compared with healthy subjects in the exhaled exosome [43]. Garbacki observed miR-574-5p was down-regulated (fold change = 0.37) with short-term exposure to an allergen but upregulated (fold change = 13.18) with long-term exposure to allergen in the mouse model of asthma, implying a cell cycle regulatory function [44]. Gomez built a miRNA and mRNA network based on the sputum of patients with asthma, and found the hsa-miR-574-5p module positively correlated with eosinophil counts (p = 0.008) and negatively correlated with bronchodilator response (p = 0.04) [45]. Moreover, several studies had shown that 25-hydroxyvitamin D (25OHD) or vitamin D3 could alter miRNA expression [46,47]. Hsa-miR-574-5p was down-regulated (fold change = 1.79) in the high vitamin D group in the plasma of pregnant mothers [48]. Together, these studies strongly suggest that hsa-miR-574-5p may be a vitamin D-related mediator and biomarker in asthma.
Our study also showed that hsa-miR-151a-5p had the strongest protective effect on incident asthma in the high vitamin D treatment group (Figure 1). Francisco-Garcia's group examined miRNAs in nanovesicles from bronchoalveolar lavage of severe asthmatic patients and they reported that hsa-miR-151a-5p positively correlated with FEV1% (prebronchodilator forced expiratory volume in one second as a percent predicted) (r = 0.48, p = 0.03), which confirmed with our analysis [49]. Jorde et al. also observed that hsa-miR-151a-5p was up-regulated in plasma after 1 year of 40,000 IU vitamin D3 per week [46]. In summary, hsa-miR-151a-5p may also be a potential mediator and biomarker.
For other validated miRNAs, in our study, many of them were associated with allergic inflammation. MiR-342-3p was associated with allergic airway disease in a mouse population [53] and was also observed to suppress inflammation response in human macrophages THP-1 cells [54]. MiR-122-5p was increased in extracellular vesicles from subjects with asthma [55] and decreased in asthmatic bronchial epithelial cells [56]. MiR-193b-5p played an important role in virus induced lung injury [57] and miR-125b-2-3p might affect the G2/M phase of the cell cycle [58].
Using our 9 validated miRNAs, we built a miRNA-target gene network using the data in miRTarBase [29] and performed a KEGG pathway enrichment analysis through miRNet [28]. Eight signaling pathways were identified with an FDR less than 0.05. Of these, cell cycle was the most significant signaling pathway (FDR = 4.8 × 10 −7 ) and some studies demonstrated that the airway smooth muscle (ASM) cell proliferation was increased in asthmatic patients [30]. The TGF-beta pathway (FDR = 4.12 × 10 −5 ) has been widely investigated and associated with the airway remodeling process in asthma [31,32]. EGFR was observed to be highly expressed in bronchial epithelial cells in asthma [33] and the ErbB signaling pathway (FDR = 1.88 × 10 −4 ) plays a key role in mediating airway hyperresponsiveness and remodeling in a chronic allergic mouse model [34]. The Wnt signaling pathway (FDR = 1.23 × 10 −3 ) was associated with impaired lung function in childhood asthma [35] and down-regulated by vitamin D, leading to alleviation of airway remodeling [36]. The JAK-STAT (Janus kinase/signal transducers and activators of transcription) signaling pathway (FDR = 2.03 × 10 −4 ) was involved in cytokine related regulation and could induce polarization of T helper cells, of which Th2 cells were believed to play an important part in initiating the airway inflammatory response [37,38]. As the central role in the p53 signaling pathway (FDR = 5.77 × 10 −3 ), p53 was reported for increased expression in bronchial smooth muscle (BSM) from asthmatic patients and associated with BSM proliferation and mitochondrial biogenesis [39]. T cell receptor signaling pathway (FDR = 3.18 × 10 −2 ) contained complex signaling cascades and was crucial for T cell development [40]. The Toll-like receptor signaling pathway (FDR = 4.55 × 10 −2 ) has been associated with innate immune system and could induce a pro-inflammatory response, resulting in airway inflammation [41,42].
Our study has several limitations most notably the small sample size of both the discovery and validation cohorts. However, through the use of meta-analysis, we demonstrated that all the target miRNAs are significant in their effects. Additionally, most of these miRNAs are associated with allergic inflammation, which is also supported by the KEGG pathway enrichment analysis. Of all the target miRNAs, hsa-miR-574-5p and hsa-miR-151a-5p are the strongest potential mediators and biomarkers of asthma and hsa-miR-574-5p has an excellent prognostic power individually.

Conclusions
In summary, our findings show that circulating microRNAs are associated with, and predictive of, incident asthma under the effect of vitamin D treatment. Circulating microRNAs may be good biomarkers and mediators of asthma.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/jpm11040307/s1, Table S1: Significant miRNAs in Main Analysis, Table S2: Stratified analysis in Placebo Group, Table S3: Validated significant miRNAs in Project Viva, Table S4: KEGG enrichment analysis of target genes of validated miRNAs. Figure   Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The datasets used and analyzed are available from the corresponding author on reasonable request.