Identification of microRNAs associated with allergic airway disease using a genetically diverse mouse population

Background Allergic airway diseases (AADs) such as asthma are characterized in part by granulocytic airway inflammation. The gene regulatory networks that govern granulocyte recruitment are poorly understood, but evidence is accruing that microRNAs (miRNAs) play an important role. To identify miRNAs that may underlie AADs, we used two complementary approaches that leveraged the genotypic and phenotypic diversity of the Collaborative Cross (CC) mouse population. In the first approach, we sought to identify miRNA expression quantitative trait loci (eQTL) that overlap QTL for AAD-related phenotypes. Specifically, CC founder strains and incipient lines of the CC were sensitized and challenged with house dust mite allergen followed by measurement of granulocyte recruitment to the lung. Total lung RNA was isolated and miRNA was measured using arrays for CC founders and qRT-PCR for incipient CC lines. Results Among CC founders, 92 miRNAs were differentially expressed. We measured the expression of 40 of the most highly expressed of these 92 miRNAs in the incipient lines of the CC and identified 18 eQTL corresponding to 14 different miRNAs. Surprisingly, half of these eQTL were distal to the corresponding miRNAs, and even on different chromosomes. One of the largest-effect local miRNA eQTL was for miR-342-3p, for which we identified putative causal variants by bioinformatic analysis of the effects of single nucleotide polymorphisms on RNA structure. None of the miRNA eQTL co-localized with QTL for eosinophil or neutrophil recruitment. In the second approach, we constructed putative miRNA/mRNA regulatory networks and identified three miRNAs (miR-497, miR-351 and miR-31) as candidate master regulators of genes associated with neutrophil recruitment. Analysis of a dataset from human keratinocytes transfected with a miR-31 inhibitor revealed two target genes in common with miR-31 targets correlated with neutrophils, namely Oxsr1 and Nsf. Conclusions miRNA expression in the allergically inflamed murine lung is regulated by genetic loci that are smaller in effect size compared to mRNA eQTL and often act in trans. Thus our results indicate that the genetic architecture of miRNA expression is different from mRNA expression. We identified three miRNAs, miR-497, miR-351 and miR-31, that are candidate master regulators of genes associated with neutrophil recruitment. Because miR-31 is expressed in airway epithelia and is predicted to target genes with known links to neutrophilic inflammation, we suggest that miR-31 is a potentially novel regulator of airway inflammation. Electronic supplementary material The online version of this article (doi:10.1186/s12864-015-1732-9) contains supplementary material, which is available to authorized users.


Background
Post-transcriptional control of gene expression is critical for the proper control of inflammation [1]. microRNA (miRNA) mediated transcript degradation and translational inhibition represent two important mechanisms of post-transcriptional control, and a large body of work demonstrates the complex role miRNAs play in finetuning the regulatory networks that govern inflammation and innate and adaptive immunity [2]. Allergic airway diseases (AADs) such as allergic rhinitis and asthma are characterized by aberrant adaptive immune responses to allergens that results in airway inflammation, mucus hyper-secretion, and airway hyper-responsiveness [3]. Evidence is accruing that miRNAs play a role in these diseases [4][5][6]. In particular, several miRNAs have been shown to regulate pathways underlying AAD pathogenesis including T-lymphocyte polarization and function [7][8][9], regulatory T-cell activity [10], and airway smooth muscle proliferation and contractility [11][12][13]. Previous studies have also identified miRNAs that are differentially expressed in airway epithelia of individuals with AAD [14][15][16] and in mouse models of AAD [8,9,17]. However, the question of whether genetic variation in miRNA expression affects AAD phenotypes has not been addressed.
We previously characterized the genetics of lung mRNA expression in a mouse model of AAD [18] using the Collaborative Cross (CC) mouse genetics reference population. The CC is a mammalian systems genetics resource comprised of recombinant inbred lines derived from eight-way crosses using five classical inbred strains (C57BL/6 J, 129S1/SvImJ, A/J, NOD/ShiLtJ, and NZO/ H1LtJ) and three wild-derived inbred strains (WSB/EiJ, PWK/PhJ, and CAST/EiJ) [19,20]. Using incipient (i.e., not yet fully inbred) CC lines, which we refer to as "pre-CC" mice, we identified more than 6,000 gene expression quantitative trait loci (eQTL), the majority of which were located near the cognate gene. We used these eQTL to identify candidate genes within QTL for inflammation phenotypes [21]. For example, we identified a locus on chromosome 7 that controls the expression of a gene that is causally related to neutrophil recruitment [22], thereby linking sequence variation, gene expression, and a downstream phenotype. Those analyses were limited to mRNAs. Thus, the role of miR-NAs in this experimental system remains unknown.
To identify miRNAs that play a role in AAD, we leveraged the genotypic and phenotypic diversity of preCC mice and utilized two complementary analytical approaches. In the first approach we assessed whether regulatory variation at miRNA loci underlies variation in allergic inflammation. We identified robustly expressed miRNAs that vary by strain, mapped the genetic loci that correlate with expression levels, and then determined whether these loci are associated with AAD phenotypes. In the second approach, we combined bioinformatic sequence analysis with statistical simulations to identify miRNAs that serve as candidate master regulators of genes altered during airway inflammation.

Methods
Our overall study design and workflow is depicted in Additional file 1: Figure S1.

Mice
We obtained 129 male preCC mice (ages 10-14 weeks) from Oak Ridge National Laboratory [21,23,24]. Each mouse was from an independent CC line that had undergone five to fourteen generations of inbreeding. We also obtained four male mice of each of the eight CC founder strains from The Jackson Laboratory. All mice were singly housed, with alpha-dri bedding, under normal 12-h light/dark cycles. All experiments conducted with mice in this study were compliant with an Institutional Animal Care and Use Committee protocol at an animal facility approved accrediated by the Association for Assessment and Accreditation of Laboratory Animal Care International.

Phenotyping protocol
We employed a house dust mite (HDM) model of AAD that produces hallmark disease phenotypes including T h 2biased airway inflammation, elevated serum IgE, mucous cell metaplasia, and airway hyper-responsiveness in a strain-dependent fashion [18]. Mice were sensitized with 10 μg of the immunodominant allergen from the Dermatophagoides pteronyssinus species of HDM, Der p 1, by intra-peritoneal injection on days 0 and 7, followed by challenge with 50 μg of Der p 1, administered by oropharyngeal aspiration, on day 14. On day 17, mice were euthanized, followed by collection of whole lung lavage fluid with two successive volumes of 0.5 and 1.0 ml PBS. No perfusion was performed, so vascular contents were still present in the lungs. Following lavage, lung tissue was snap frozen. Cells in lavage fluid were isolated by centrifugation; eosinophil and neutrophil counts were then manually determined using cytospins and morphologic criteria.

miRNA expression analysis
For the eight CC founder strains (n = 4/strain), we isolated total lung RNA by Trizol extraction. RNA quality was assessed using an Agilent Bioanalyzer. With one exception, all samples had RNA integrity numbers greater than 7. RNA samples were then processed and hybridized to Affymetrix miRNA 2.0 arrays (GSE63954). We limited our analysis to the 723 probe sets on the array that are specific to mouse miRNAs. Manual inspection of results from principle component analysis of miRNA expression revealed that two samples were outliers (one A/J and one NZO/H1LtJ sample); these two samples were removed prior to subsequent analyses. We used an ANOVA model to identify miRNAs that were differentially expressed by strain at a false discovery rate (FDR) q-value < 0.05. We further filtered this list of miRNAs to those that were highly expressed, which we defined as within the top 20 th percentile of miRNA expression in any one CC founder strain, and those that varied by at least 1.5-fold between the highest and lowest-expressing strains. This limited the list to 38 miRNAs. We also selected one miRNA, miR-17* (now known as miR-17-3p), as a reference miRNA for normalization in experiments with preCC mice. miR-17* was selected because its expression was near the mean of all miRNAs that met our selection criteria (mean of all miRNAs = 7.3; mean of miR-17* = 7.1) and it was not differentially expressed by strain (p = 0.89). From these ANOVA models described above, we estimated broad-sense heritability (H 2 ) by calculating the interclass correlation (r 1 ) and the coefficient of genetic determination (g 2 ) [25]: where MSB and MSW are the mean squares between and within, respectively, from the eight-way ANOVA model described above and n is the number of mice per strain. g 2 is a slightly modified estimate of heritability that accounts for the doubling of the additive genetic variance with inbreeding and is calculated as: For the preCC mice, the same RNA isolation protocol was used. We then used Exiqon locked nucleic acid (LNA)based qRT-PCR assays for each of the 38 miRNAs, except miR-805 (because an assay was not available from Exiqon). We added miR-148b and miR-182 due to published evidence indicating roles in asthma or allergic inflammation in the lung [10,26], and miR-17* and U6 for normalization. Thus in total we measured the expression of 40 miRNAs and one snRNA (U6). Since the time of ordering the Exiqon LNA qRT-PCR assays, three of the measured miRNAs have had name changes: miR-17* is now miR-17-3p; miR-193* is now miR-193-5p; and miR-322* is now miR-322-3p. We maintained the old nomenclature throughout the manuscript since that is consistent with Exiqon product information at the time of ordering. Exiqon miRCURY LNA Universal RT was used to generate cDNA, after which miRNA-specific primers were used to amplify each miRNA. For each miRNA for each mouse, we calculated the delta Ct relative to miR-17*. Across all samples and miRNAs, the delta Ct relative to miR-17* and U6 snRNA (another commonly used control for normalization) were highly correlated (Pearson r = 0.81, p < 5 × 10 −16 ). qPCR data is provided as Additional file 2.

eQTL approach Genotyping and eQTL mapping
We genotyped each mouse at the University of North Carolina at Chapel Hill, using one of two Affymetrix SNP arrays (A or B) that were produced during the development of the Mouse Diversity array (MDA) [27]. After removing uninformative and poorly performing SNPs, these arrays contained 181,752 (A-array) and 180,976 (B-array) SNP assays, and the set of SNPs on each array did not overlap. Most mice (83 %) were genotyped on the B-array and the remaining were genotyped on the A-array. These training arrays were annotated to NCBI Build 36 of the mouse genome, but we mapped QTL boundaries to Build 37 positions to integrate with other resources. We report NCBI Build 37 positions in our results. We estimated the most probable ancestor for each SNP in each mouse using the GAIN algorithm [28], and reconstructed founder haplotypes based on these results. We then merged the non-overlapping SNP datasets from arrays A and B by imputing unobserved genotypes based on inferred founder haplotype. For QTL mapping, we used HAPPY [29] to infer ancestry matrices for an additive genetic model. For computational efficiency, we then averaged the matrices across SNPs between which GAIN inferred no recombination in the population, and this reduced the mapping dataset to 27,059 intervals. Genotype data are provided as Additional file 3. We used BAGPIPE [30] to fit a regression model and calculate LOD scores. Significance thresholds were determined by permutation (n = 250 permutations per miRNA), and we used the 1.5 LOD drop method to approximate confidence intervals for QTL [31]. We defined local eQTL as eQTL located within 10 Mb of the cognate gene based on prior work on mRNA eQTL [21]. Narrow sense heritability (h 2 ) for eQTL was estimated by regression of miRNA expression on CC founder haplotype probabilities. For example, the expression of a given miRNA Y is where each strain term represents the estimated number of that strain's haplotypes present at the QTL for a given mouse (this estimate being the posterior expectation derived from the probabilistic haplotype reconstruction using HAPPY [29]), and where A/J is designated as the reference strain. We used the model R 2 to describe the percent of phenotypic variation accounted for by the locus.

Identification of regions of shared ancestry and phylogenetic analysis
We used a comparative genomic approach to narrow QTL regions and identify candidate genes, as described previously [24]. Using the allele effect plots as a guide, we grouped founder alleles into two groups by effect (e.g., high vs. low allele effect). We divided the two groups based on the greatest difference between ordered allele effects estimated at the QTL peak. We then used the Collaborative Cross Viewer at the University of North Carolina [20,32] and genome sequences from the Wellcome Trust Sanger mouse genomes project (MGP) [33,34] to identify regions in the confidence interval in which grouped strains have shared SNP genotypes. Positions identified as heterozygous and low confidence genotype calls were omitted from the analysis.

Analysis of the effects of SNPs on pre-miRNA structure
To examine the consequences of SNPs on miRNA structure, we used two computational approaches. First, we used RNAfold [35] to compare the predicted structures of pre-miR-342 based on the reference allele (C57BL/6J) to a musculus-derived allele (shared by the NOD/ShiLtJ and PWK/PhJ strains) containing two SNPs (rs264778660 and rs242689107) and one indel (rs261236356) that fall within the precursor. Second, we used SNPfold [36] to formally test whether either SNP affects the ensemble of predicted RNA structures.
Bioinformatic/Statistical approach to identify putative miRNA-mRNA regulatory networks To construct putative miRNA-mRNA regulatory networks, we utilized a multistep procedure outlined in green in Additional file 1: Figure S1. We first identified mRNA transcripts linearly (positively or negatively) correlated at an FDR < 0.1 with log-transformed eosinophil or neutrophil counts that were previously quantified [21,22]. We refer to these as called quantitative trait transcripts [37]. Second, we identified quantitative trait miRNAs, which we defined as miRNAs that were linearly correlated with inflammation phenotypes at an FDR < 0.1. For each phenotype (eosinophil or neutrophil counts), we paired positive quantitative trait transcripts with negative quantitative trait miRNAs and vice-versa. We then used a Monte Carlo simulation strategy called miRhub [38,39] to identify the miRNAs that may act as regulatory hubs for a given set of oppositely correlated quantitative trait transcripts. miRhub identifies miRNAs that are predicted to regulate a target gene list/network more than expected by chance. To accomplish this the algorithm utilizes the following inputs: (1) the strength and clustering of miRNA target sites in the gene set (as predicted by the TargetScan algorithm [40]), these predicted interactions can be filtered to require a minimum level of conservation across species -we specified conservation between mouse and one other species (rat, human, or chicken); (2) the centrality of target genes in a network, which is based on high confidence protein-protein interactions listed in the STRING 9.0 database (http://string-db.org/) -we did not employ this feature; and (3) the expected range of targeting scores in random gene lists/networks (n = 1000). Empirical p-values are generated based on the comparison of the observed distribution of target scores to the distribution found in the random gene networks, and are adjusted for multiple comparisons using a false discovery rate.

Results
As outlined in Additional file 1: Figure S1, we used two complementary approaches to identify miRNAs that may be associated with AAD, one based explicitly on genetically determined variation in miRNA expression ("eQTL approach" shown in blue), and one that does not rely on genetic diversity per se but rather combined bioinformatics and statistical analyses to construct putative miRNA-mRNA regulatory networks ("bioinformatics/statistical approach" shown in green). We begin with results from our survey of miRNA expression in CC founder strains, followed results from the eQTL and bioinformatics/statistical approaches, respectively.

Detection of differentially expressed miRNAs in the lungs of CC founder strains
We measured the expression of 723 miRNAs in whole lung RNA samples from the eight CC founder strains (n = 3-4/strain) using microarrays (Methods). We identified 92 that were differentially expressed by strain at an FDR < 0.05 (Additional file 4: Table S1), 38 of which were identified as highly expressed in at least one strain (Methods, Additional file 5: Table S2). Hierarchical clustering of these 38 miRNAs clearly demonstrated the effect of strain ( Fig. 1) and the sample clustering was consistent with the known phylogenetic relationships among these strains [41]. We estimated the broad sense heritability (H 2 ) for miRNA expression (Additional file 5: Table S2); H 2 values ranged from 0.28 for miR-200a to 0.94 for miR-342-3p.
To ensure that these results were not false positives due to altered hybridization between array probes and genetic variants, we mapped probesets of differentially expressed miRNAs to genetic variants contained in the Sanger Mouse Genomes Project data [33,34]. Three probesets aligned to regions that contain structural variants among CC founder strains (miR-148b, miR-192, and miR-194), but the observed patterns of expression were not correlated with the strain distribution of structural variants. Thus we conclude that the variation in miRNA expression is not biased by the array platform we used.

eQTL Approach
We then measured the expression of 37 out of 38 of these miRNAs (Methods), as well as miR-148b and miR-182 due to published evidence indicating roles for these miRNAs in asthma or allergic inflammation in the lung [10,26], in 129 preCC lines using Exiqon LNA-based qRT-PCR assays. We observed a broad range of expression for many miRNAs (Additional file 1: Figure S2). We note in particular the bimodal distributions for two miRNAs, miR-133 and miR-489, suggesting a single, large-effect eQTL for each.
Next we performed genome-wide scans for each miRNA to identify miRNA eQTL. In total, we identified 18 highconfidence eQTL (at p adjusted < 0.05) corresponding to 14 different miRNAs, and an additional eight potential eQTL (at p adjusted < 0.10) for seven more miRNAs (Table 1). This equates to high-confidence eQTL for 38 % of the miRNAs that were differentially and highly expressed among the founder strains (and 35 % of the total number of miRNAs studied), though surprisingly miR-133 was not one of them. Of the high-confidence eQTL, nine were located near the miRNA gene itself (within 10 Mb), which we refer to as local eQTL [21], while nine were located distal to the miRNA gene, including on other chromosomes. These results differ from our previous studies of protein coding gene eQTL, which revealed predominately (~75 %) largeeffect, local eQTL (Fig. 2) [21].
We categorized the miRNA eQTL in terms of their effect size and location in the genome. Two eQTL, for miR-489 and miR-342-3p, stood out in terms of effect size. Not surprisingly, expression of these two miRNAs also showed the highest heritability in CC founder lines (Additional file 5: Table S2). We confirmed that for both -3.24 0.00 3.24 Fig. 1 Hierarchical clustering miRNA expression among CC founder lines. Of the 92 differentially expressed miRNAs detected using a microarray platform, 38 were selected based on the expression values and these are depicted here. Note that with one exception (WSB/EiJ mouse at far right), the overall sample clustering is consistent with the phylogenetic relationships among these strains miRNAs the expression levels among preCC mice with a specific CC founder haplotype at the eQTL were concordant with expression levels in the CC founder strain (Figs. 3 and 4). We calculated the narrow sense heritability (h 2 ) to be 0.60 and 0.84 for the miR-489 and miR-342-3p eQTL, respectively. Taken together, these results indicate that cis-regulatory elements at the respective eQTL are the primary determinants of miR-489 and miR-342-3p expression in the lung. However, we note that miR-489 expression in CC founders suggests three distinct groups (high expression group composed of NOD/ShiLtJ and WSB/EiJ; medium expression group composed of A/J, C57BL/6J, 129S1SvImJ, and NZO/H1LtJ; and low expression group composed of CAST/EiJ and PWK/PhJ), but the allele effects for the eQTL on chromosome (Chr) 6 indicates that this eQTL only explains the differences between CAST/EiJ and PWK/PhJ vs. the other six CC founders. Thus it appears that one or more additional loci likely contribute to miR-489 expression; but we did not  The allele effects for miR-342-3p indicated that the NOD/ ShiLtJ and PWK/PhJ alleles were similar in terms of their effect on miRNA expression (Fig. 4b). Using publically available array-based genotype and haplotype data from the CC founder strains [20,32], we found that, like the musculus strain PWK/PhJ, the NOD/ShiLtJ founder strain has a musculus-derived haplotype in the QTL confidence interval (Chr 12: 108,668,025-110,004,769 bp). We confirmed that  Fig. 4 The miR-342-3p eQTL. a. Expression of miR-342-3p in CC founders (a) and preCC mice (b) as a function of founder haplotype at the eQTL on chromosome 12. CC founder miRNA expression was measured by microarray while preCC miRNA expression was measured by qRT-PCR and data for the latter are presented as −1*DeltaCq (normalized to miR-17*). c. Phylogeny of CC founder strains based on SNP data for the region on Chr 12 containing the miR-342-3p eQTL. Bootstrap values were greater than or equal to 96 (out of 100) for each branch of tree.The red line denotes the branch that contains the putative causal variant these two strains share a haplotype in this region using more complete single nucleotide polymorphism (SNP) data from the Sanger Mouse Genomes Project [33] (Fig. 4c). This finding facilitated the identification of putative causal variants descended from the musculus sub-species that are shared by NOD/ShiLtJ and PWK/PhJ and are different from all other founder strains. We identified 938 SNPs in the region that met these criteria, the majority of which were located in or near 11 genes (defined as +/− 5 kb from the gene), including 20 SNPs in or near the miR-342 locus. Two SNPs (rs264778660 and rs242689107) and one indel (rs261236356) were in the miR-342 precursor (pre-miR-342) and thus were considered strong candidate causal variants for the eQTL because the stemloop structure of the precursor is critical for miRNA maturation and expression. The RNAfold algorithm [35] predicted an allele-dependent effect for both SNPs (but not the indel) on the terminal loop of the pre-miR-342 stem-loop structure (Fig. 5). Subsequent analysis using SNPfold [36] provided further in silico evidence that rs264778660 significantly alters the conformation of pre-miR-342 (p = 0.01, Additional file 1: Figure S3). These results suggest that at least rs264778660 is likely to cause a change in pre-miR-342 structure, which in turn may alter processing by Dicer and Drosha, and consequently affect expression of both miR-342-3p and miR-342-5p. Indeed, among CC founder strains, we found a strong correlation between miR-342-3p and miR-342-5p expression, with NOD/ShiLtJ and PWK/PhJ strains having low expression for both arms of miR-342 (Additional file 1: Figure S4). We noted that indel (rs261236356) is predicted to cause a change in the seed sequence of miR-342-5p (Fig. 5), potentially altering the relationship between this miRNA and target genes. We did not measure the expression of miR-342-5p in the preCC because its expression levels were very low. Seven other miRNAs had local eQTL: miR-146b, miR-181d, miR-187, miR-203, miR-221, miR-322, and miR-503. The allele effects for these eQTL are shown in Additional file 1: Figure S5. Some miRNAs had both local and distant eQTL. For example, miR-322 and miR-503 each had one local eQTL (on Chr X) and one trans-eQTL on Chr 11. We built regression models for these two miRNAs and found that these two loci accounted for 44 % and 37 % of variation in miR-322 and miR-503 expression, respectively.
The trans-eQTL locus shared by miR-322 and miR-503 was also weakly associated with the expression of miR-351 and miR-497 (p adjusted < 0.1). We found that each pairwise comparison of miRNA expression for all four miRNAs was highly significant (Additional file 6: Table S3) and that the allele effects for all four of these trans-eQTL were consistent (Fig. 6). This suggests that a single variant mediates allele-specific trans-regulation of the expression of all four miRNAs. miRNA expression in mice with the 129S1/SvImJ allele was significantly different from all other strains with the exception of WSB/EiJ (p < 0.05 by Tukey's honestly significant difference test), indicating the 129S1/SvImJ allele is functionally distinct from most other CC founder strains. This region contains 100 genes and we looked for SNPs for which the 129S1/SvImJ strain is distinct from all other CC founders excluding WSB/EiJ and identified 212 SNPs in 20 genes that met this criterion: A830035A12Rik, Abca6, Amz2, C330019F10Rik, Cacng1, Cacng5, Ccdc46, Gm11650, Gm11655, Gm11674, Gm11677, Gm11678, Gm11680, Gm11681, Gm11714, Gna13, Helz, Prkca, Smurf2, and Tex2 (Additional file 7: Table S4). These genes represent causal candidates for the miRNA trans-eQTL. We identified a second set of trans-eQTL on Chr 9 at~107-108 Mb for miR-486 and miR-451, but the allele effects for these eQTL were not completely consistent (Additional file 1: Figure S6), thus we cannot conclude that these two trans-eQTL are one and the same.

The relationship between miRNAs and inflammation phenotypes
Genetic variation in a miRNA has been linked to disease phenotypes [42][43][44][45]. We asked whether the miRNA Bioinformatic/Statistical approach to identify putative miRNA-mRNA regulatory networks As shown in Additional file 1: Figure S1, we used a complementary approach to ask whether any miRNAs were candidate master regulators of genes altered during airway inflammation. We first identified miRNAs whose expression levels were linearly correlated with eosinophils and neutrophils. At an FDR < 0.1, thirteen and nineteen miR-NAs were correlated, respectively. The vast majority was negatively correlated (Additional file 8: Table S5). We then used the miRHub algorithm (Methods), which determines whether the predicted regulatory effect of any given miRNA on a specific set of genes is significantly greater than expected by chance (i.e., acts as a "regulatory hub"), to predict whether any of these miRNAs are candidate drivers of the gene expression profiles associated with the phenotypes of interest. For the set of genes that were positively correlated with neutrophils (n = 674 at FDR < 0.1), we identified miR-497, miR-351 and miR-31 as candidate regulatory hubs (Fig. 7). No miRNA hubs were identified for the sets of genes positively (n = 1802) or negatively (n = 1605) associated with eosinophils (Additional file 1: Figure S7). These results indicate that miR-497, miR-351 and miR-31 may serve as mediators of neutrophilic inflammation by targeting genes that regulate neutrophil recruitment to the airways; predicted targets of each miRNA are listed in Additional file 9: Table S6.
miR-31 has been previously reported as a marker and/ or regulator of other inflammatory conditions, such as inflammatory bowel disease [46] and psoriasis [47]. In a previously published dataset of gene expression from human keratinocytes transfected with a miR-31 inhibitor, 96 genes were up-regulated, and of these, 24 were predicted targets of miR-31 [47]. We found that the mouse orthologs of two of these genes, Oxsr1 and Nsf, were also present in our list of predicted miR-31 targets correlated with neutrophils, resulting in statistically significant enrichment (p = 4 × 10 −3 by hyper-geometric test). While an overlap of two genes is not large, it is appreciable given the differences in cell types, species, and disease processes examined. We also found that the expression of Oxsr1 and Nsf was negatively correlated with miR-31 expression (r = −0.26, p = 2.0 × 10 −3 and r = −0.23, p = 8.7 × 10 −3 , respectively). Thus we conclude that miR-31 is likely to be an important modulator of neutrophilic airway inflammation in part by targeting Oxsr1 and Nsf.

Discussion
We examined whether regulatory variation at miRNA loci underlies variation in response to allergen sensitization and challenge using a genetically diverse population of mice. The phenotypic variation across these mice enabled the identification of eQTL for 14 out of the 40 miRNAs (35 %) that we studied, as well as suggestive eQTL for an additional seven miRNAs. We characterized these eQTL in terms of effect size and location relative to the miRNA gene and found that we could explain a large fraction of the total variation in expression for several miRNAs using eQTL. Some miRNAs exhibited appreciable heritability estimates based on gene expression in the CC founder strains, but did not have detectable eQTL in the preCC population (e.g., miR-133). One possible explanation is that variation in expression is due to epistasis among two or more loci that lack sizeable individual effects. Due to the limited power afforded by the sample size of our study, we did not test for epistatic eQTL. Additionally, due to the fact that the qRT-PCR system we used was developed by a company and the primer sequences are proprietary, we cannot rule out the possibility that false negative and/ or false positive eQTL we detected in the preCC mice are due to SNPs (or other types of genetic variants) that could affect the measurement of miRNA expression.
In the case of miR-342-3p, we were able to identify highly plausible candidate causal variants using a combination of sequence data and structural modeling. We found that a SNP (rs264778660) located in the miRNA precursor is likely to alter the terminal loop of the pre-miR-342 structure, which in turn may alter processing and maturation of miR-342-3p. Liu et al. have shown that the terminal loop of a pre-miRNA affects processing by Dicer [48]. We suggest that rs264778660 affects the Out of all the other miRNAs we studied, very few had variants located in the precursor sequence and none had sequence variants located in seed regions (except for miR-342-5p which was measured in the CC founder lines but not in the preCC population). This is consistent with other population genetic studies, which have concluded that evolutionary pressures have preserved regulatory networks involving miRNAs by selecting against mutations that cause changes in miRNA targeting [49,50].
One of the interesting results from our eQTL analysis was that distal or trans-eQTL were quite common for miR-NAs. This contrasts with typical results for mRNA eQTL, including our own findings in the same set of samples [21], where local eQTL constitute the majority. Our findings are similar to those recently reported for miRNA eQTL in human peripheral blood, where miRNA eQTL were often located far upstream from the miRNA itself [51]. In one case, we identified four miRNAs whose expression was regulated by a single locus on Chr 11. Of the~20 candidate genes in the region, Prkca, which encodes Protein Kinase C-alpha, is a plausible candidate because it is a well-known regulator of immune responses in the lung [52] and therefore may regulate pathways that in turn affect miRNA expression.
Trans-eQTL were also more common than local eQTL in a previous miRNA eQTL mapping study involving an F2-intercross with diabetes resistant and susceptible strains (C57BL/6J and BTBR, respectively, each containing two mutant alleles of the leptin gene) [53]. In that study, the authors measured the expression of 220 miRNAs in adipose tissue in 290 F2 mice and identified eQTL for 21 of them (i.e.,~10 % of miRNAs), three of which were local or cis-eQTL. Thus, in comparison, we identified more eQTL (in terms of number of eQTL relative to the number measured) and a larger proportion of local eQTL. It is likely that this can be attributed to the greater genetic diversity present in the CC founder strains versus the two founder strains used in that F2 mapping study. Five miR-NAs or miRNA family members were measured in that study and were also a part of this study, namely miR-31, miR-34c, miR-148a, miR-181, and miR-451. We detected a similar local eQTL for miR-181 but there was no overlap in eQTL for the other miRNAs. These differences could be explained at least in part by different regulatory elements responsible for controlling miRNA expression in different tissues (adipose vs. lung) and in response to different perturbations (ob transgene/diet vs. allergen). In aggregate, our results and those of Zhao et al. [53] suggest that miRNA expression levels in mice are more strongly dependent on genetic variation at distal loci than mRNAs.
Additionally, the LOD scores for local miRNA eQTL were, on average, lower in magnitude than those for local mRNA eQTL we previously identified. Hence local eQTL account for a smaller proportion of miRNA expression variation than local mRNA eQTL. This too is consistent with results from a recent human miRNA eQTL study in which local eQTL accounted for~1 % of variation in miRNA expression compared to 30-50 % for mRNA eQTL [51]. Our results indicate that the genetic architecture of miRNA expression is clearly different from mRNA expression, likely involving more loci with additive or perhaps even multiplicative effects. As such, it is perhaps not surprising that we did not identify a miRNA eQTL that underlies variation in allergic inflammation in our model of asthma. For a single miRNA eQTL to underlie a trait QTL, the miRNA eQTL must have a large effect on miRNA expression and be causally linked to the trait. While we identified miRNAs that are correlated with inflammation phenotypes, none of these miRNAs had a single major effect eQTL that explains the resulting variation in inflammation. However, additive or multiplicative effects across miRNA eQTL may explain a greater proportion of the variation in inflammation. More complex statistical models will be required to characterize the effects of multiple loci, and this will inherently require greater statistical power.
The bioinformatic analysis we conducted to identify miRNAs that may act as key regulators of genes involved in granulocyte (eosinophil and neutrophil) recruitment pointed to three miRNAs of interest for neutrophils, namely miR-497, miR-351 and miR-31. Of these, miR-31 is potentially the most promising based on prior work demonstrating that miR-31 is highly expressed in airway epithelia [54], and in skin epithelia regulates the expression of pro-inflammatory chemokines by targeting genes that regulate NF-kB activity [47]. Additionally, miR-31 is not highly expressed in neutrophils at baseline or after bacterial challenge [55], suggesting that miR-31 is not simply a surrogate metric of neutrophil counts. Based on our finding of significant enrichment of predicted miR-31 targets in our data compared to dataset from human keratinocytes in which miR-31 was inhibited [47], we suggest that two genes, Oxsr1 and Nsf, may be novel regulators of neutrophilic inflammation in the airway. NSF (N-ethylmaleimide-sensitive factor) is a cytoplasmic ATPase that plays an important role in the exocytosis of inflammatory mediators by disassembling the soluble-N-ethylmaleimide-sensitive-factor accessory-protein receptor (SNARE) complex after exocytosis [56]. Using a mouse model of peritonitis, Morrell et al. showed that inhibition of NSF resulted in diminished neutrophil recruitment to the peritoneum [57]. Given these data, one simple hypothesis is that NSF in airway epithelium regulates the release of a chemokine involved in neutrophil recruitment, and this pathway is targeted by miR-31.
Our results also indicate that other genes with established relationships to neutrophil recruitment are predicted target genes of miR-31. We note three genes in particular: Ccr2, a chemokine receptor linked to monocyte and neutrophil recruitment to the lung and other organs [58][59][60], Tnip1 (Abin1), which is involved in toll-like receptor mediated innate immune responses [61], and Unc93b1, a gene involved in antigen presentation by dendritic cells [62]. In combination with prior work in human airway and skin epithelia, our result suggests that miR-31 may be targeting genes both in the airway epithelia and in white blood cells (neutrophils, monocytes, and dendritic cells), and provide rationale for additional studies on the role of miR-31 in neutrophilic airway inflammation. We note that many previous studies have documented an important role for miR-31 in multiple cancer processes [63]; thus our findings suggest that miR-31 is microRNA with pleiotropic effects.
Conclusion miRNA expression in the allergically inflamed murine lung is regulated by genetic loci that are smaller in effect size compared to mRNA eQTL and often act in trans. Thus our results indicate that the genetic architecture of miRNA expression is different from mRNA expression. We identified three miRNAs, miR-497, miR-351 and miR-31, that are candidate master regulators of genes associated with neutrophil recruitment. Because miR-31 is expressed in airway epithelia and is predicted to target genes with known links to neutrophilic inflammation, we suggest that miR-31 is a potentially novel regulator of airway inflammation.

Availability of supporting data
Data for array-based expression of miRNAs in Collaborative Cross founder strains is provided in NCBI's Gene Expression Omnibus: http://www.ncbi.nlm.nih.gov/geo/ query/acc.cgi?acc=GSE63954. Genotype data for preCC mice and qPCR-based measurement of miRNAs in preCC mice are provided as supporting data (Additional Files 2 and 3, respectively).