Differential Expression of Serum MicroRNAs Supports CD4+ T Cell Differentiation into Th2/Th17 Cells in Severe Equine Asthma

MicroRNAs (miRNAs) regulate post-transcriptional gene expression and may be exported from cells via exosomes or in partnership with RNA-binding proteins. MiRNAs in body fluids can act in a hormone-like manner and play important roles in disease initiation and progression. Hence, miRNAs are promising candidates as biomarkers. To identify serum miRNA biomarkers in the equine model of asthma we investigated small RNA derived from the serum of 34 control and 37 asthmatic horses. These samples were used for next generation sequencing, novel miRNA identification and differential miRNA expression analysis. We identified 11 significantly differentially expressed miRNAs between case and control horses: eca-miR-128, eca-miR-744, eca-miR-197, eca-miR-103, eca-miR-107a, eca-miR-30d, eca-miR-140-3p, eca-miR-7, eca-miR-361-3p, eca-miR-148b-3p and eca-miR-215. Pathway enrichment using experimentally validated target genes of the human homologous miRNAs showed a significant enrichment in the regulation of epithelial-to-mesenchymal transition (key player in airway remodeling in asthma) and the phosphatidylinositol (3,4,5)-triphosphate (PIP3) signaling pathway (modulator of CD4+ T cell maturation and function). Downregulated miR-128 and miR-744 supports a Th2/Th17 type immune response in severe equine asthma.


Introduction
Despite intensive efforts, the prevalence of allergies and asthma is still increasing worldwide [1]. Although the therapy for asthma has improved over the years, asthmatic condition can still lead to a sudden death [2,3]. Asthma is a complex disease and hundreds of candidate genes have been proposed for this condition [4,5]. Apart from genetic predisposition, environmental factors seem to play a crucial role in asthma development, including the exposition to indoor and outdoor allergens, such as mites or pollens and irritants like lipopolysaccharides (LPS) [6,7]. Novel techniques and methods, such as next generation sequencing (NGS) have opened a new era in asthma research. Although many genome-wide association studies have been conducted, little replication has been observed [8]. One possible explanation of this phenomenon could be the fact that multiple known phenotypes of human asthma are present and classification depends on clinical phenotypes (severity, treatment-resistance etc.), trigger (allergic/non-allergic asthma, aspirin-induced asthma etc.) and inflammatory phenotype (eosinophilic, neutrophilic, or paucigranulocytic asthma) [9]. Due to the absence of standardized Raw sequencing data quality was assessed with FastQC software [40]. Since truncated adapter sequences were reported by FastQC to be over-represented, a FASTA file was generated with the original adapter sequence plus the sequences corresponding to truncated adapters reported by FastQC for every sample. These multiple adapter sequences, as well as low quality base calls (q < 20), were trimmed with cutadapt (v. 1.8) [41] with the following options: -trim-n -a file:'sample specific adapter sequences in FASTA file' -m 15 -q 20. The resulting read length distribution was further analyzed as next quality control step. The sequencing dataset of 42 controls and 37 cases is available at the European Nucleotide Archive (ENA) [42].

Novel miRNA Identification
For the identification of potential novel equine miRNAs, the quality and adapter trimmed reads were mapped to the reference genome (EquCab2.0) [43] using miRDeep2 (v.0.0.7) [44]. Briefly, miRDeep2 predicts novel miRNAs by aligning the reads to a reference genome and extracting candidate pre-miRNA sequences from the aligned genomic DNA loci and these extracted pre-miRNA sequences are assigned a score based on the ability of the precursor to fold to a pre-miRNA like secondary hairpin structure, the absolute and relative number of reads mapping to the three distinct precursor products that result from dicer processing (5p arm, 3p arm and loop), the possible conservation of the 5 end and the presence of a 3 overhang in the mature sequence [45]. If the algorithm fails to identify a stem-loop product within the candidate, then the miRNA precursor is rejected. The processed reads from all 79 FASTQ files were mapped using the mapper.pl script of the miRDeep2 tool in one run using a configuration file containing a list of all FASTQ files names and a three-letter code for each sample. The module maps the reads to the genome with Bowtie1 [46] keeping only the alignments with 0 mismatches (option -n) in the seed region and only reads that do not map to more than five different loci in the genome are kept (option -m). The module was run with the following parameters: -d -v -u -e -h -o 16 -m -n -q -p. Subsequently, the core algorithm module miRDeep2.pl was used to discover potential novel miRNAs as described previously [29]. MiRDeep2 assigns each novel miRNA a log-odds score (further referenced as miRDeep2 score) that represents the probability that the precursor is a true miRNA precursor based on the theory of miRNA processing by dicer as well as the actual data and the alignment pattern of the reads. Any pre-miRNA candidate with a miRDeep2 score of at least 0 was considered as a potential (predicted) miRNA. To filter for predicted potential novel horse-specific miRNAs with high confidence, a miRDeep2 score cut-off of four was adopted. This cut-off was set based on our initial analysis: the signal-to-noise ratio (total miRNA precursors reported divided by the mean estimated false positive miRNA precursors over 100 rounds of permuted controls) reached a reasonable level of 9.4:1 at that threshold while the value dropped drastically for lower miRDeep2 scores. Additionally, the percentage of the detected novel miRNAs that were estimated to be true positives dropped drastically too at that threshold (82% true positive rate at a miRDeep2 score of 4 versus a rate of 54% true positive rate at a miRDeep2 score of 3; Table S2). Human known miRNAs were used as a guide for the novel equine miRNA identification. Known equine and human miRNAs were obtained from the miRBase database, release 21 [47]. Finally, overlapping novel miRNAs (by at most 1 nucleotide) were merged using BED tools [48].

Identification of a Read Count Threshold to Select Relevant miRNAs
The distribution of miRNA transcript counts in our dataset showed three distinct phases, similar to that described in Koh et al. [49]. The first phase with low count miRNAs likely represents noise in the dataset, due to random expression of these transcripts. The second phase shows similar expression between replicates which may be attributed to steady state expression or consistently expressed miRNA transcripts. Finally, the last phase consists of miRNAs with large count numbers. Kolmogorov-Smirnov test (KS-test) was used to examine the similarity of the empirical distributions between replicates and distinguish biologically relevant miRNAs from the noisy transcripts. The D statistic of a given KS-test near zero indicates similar distributions, while a larger test statistic indicates more severe bias. Similar to Koh et al. we used the D statistic as a cost function, hence we first determined the D statistic of the whole dataset and then determined the D statistic by iteratively excluding miRNAs not achieving a required minimal read count level [49]. As we changed the minimum threshold value for exclusion of low count miRNAs, we arrive at a low D statistic point indicating similar distributions. The D statistic point where this first occurs is designated as the minimum threshold count value of biological significance (code available on GitHub [50]).

Differential Expression Analysis
The dataset used for differential expression analysis consisted of a total of 71 horses: 34 controls (19 females and 15 males) and 37 affected horses with severe equine asthma (one principal component analysis (PCA) outlier and 7 phenotypically changed horses were excluded) (Table S7).
For differential expression analysis, both novel and known equine miRNAs (miRBase v21) were used. The quantification of miRNA expression level was carried out using the quantifier.pl script of the miRDeep2 tool applying the parameters -k -j.
MicroRNAs not meeting the threshold for biological relevance determined by the KS-test approach were filtered out. Furthermore, we removed lowly expressed miRNAs (showing zero counts in more than 10% of samples). The filtered miRNAs were then quality controlled using PCA (top 50 variable miRNAs) based clustering. The PCA plot was produced using variance-stabilized counts [51]. Differential expression analyses were carried out with edgeR (v. 3.16.5) [52] and DESeq2 (v. 1.14.1) [53]. For DESeq2, we used the DESeq function, which estimates size factors and dispersions and finally fits a model in order to perform differential expression tests using negative binomial generalized linear models. For edgeR, the function calcNormFactors was used to normalize the count data, while the function estimateDisp was used to estimate the dispersions and likelihood ratio test for generalized linear models (glmLRT) for the differential expression tests. For both tools, the linear model applied corrected for the levels of hemolysis as well as the underlying population structure (Table S1). MicroRNAs with a false discovery rate (FDR) adjusted p-value (adjusted p) determined by DESeq2 below our FDR threshold of 0.05 were considered as differentially expressed (up/downregulated) and will further be referenced as differentially expressed miRNAs (DEmiRs).

mRNA-miRNA Interaction
We next retrieved potential target genes of the DEmiR, showing the lowest adjusted p-value, eca-miR-128, with TargetScan (v. 6.2) [54]. Human gene symbols of the target genes were converted to horse Ensembl IDs with the Ensembl biomart tool [55]. TargetScan searches for conserved target-sites in the 3 UTR of equine genes. The human homologous gene symbols of the potential target genes were used for Gene Ontology (GO) (GO, release 106) biological process network enrichment analysis with the GeneCodis tool [56][57][58][59]. The list of reported target genes was intersected with a list of previously published differentially expressed mRNA in PBMCs of the same horses in order to investigate potential regulatory miRNA-mRNA networks [33].

Functional Enrichment of DEmiR Target Genes
A high confidence set of experimentally-validated target genes for all DEmiR (human homologues) was obtained from DIANA-TarBase v. 7.0 [60] using a stringent prediction score threshold of 0.9 (Table S8). The resulting set of high-confidence target genes was used with the GeneGo database to perform a functional enrichment analysis. MetaCoreTM v. 6.32 (Thomson Reuters, London, UK) software was used for GeneGo pathway and network analysis. The auto expand algorithm was used for building a functional network.

Code
The scripts that were used for the analyses of this study can be found on GitHub [50].

Small RNA Sequencing and Data Quality Control
As a first step, small RNAs were extracted from 79 horse serum samples ( Figure 1). Few samples were visibly hemolytic and the absorbance values at 414 nm ranged from 0.53 to 4 (median = 1). We used between 12.55 ng and 88.00 ng (mean = 36.60 ng) of the small RNA extracts for the library preparation. Subsequent to an initial quality control step using FASTQC [40], we performed adapter and quality trimming. The mean percentage of total reads remaining after applying cutadapt [41] was 86.6% and truncated adapters were trimmed from an average of 0.34% of total reads. After pre-processing of the raw reads, we obtained 11.7 million reads per library on average. The mean number of reads mapped by Bowtie1 to the horse reference genome per library was 6.99 million.

Figure 1.
Flow chart outlining the pipeline for novel miRNA detection as well as miRNA expression profiling and differential expression analysis. After extracting small RNA from serum samples of 79 horses, the samples were sequenced and analyzed using state of the art bioinformatics tools. While all samples were used for the novel miRNA identification, only 71 horses were used for the differential miRNA expression analysis.
Analysis of the read length distribution of all samples revealed a clear bimodality where on average 4.6% of the total reads constituted the first peak between 21 and 24 nucleotides (nts) while the second and more prominent peak was located between 29 and 33 nts with 95% of the total reads falling into this range ( Figure 2A). Since the percentage of mapped reads to the equine genome showed a rather low value of median 38% using Bowtie1 [46] (mapper.pl module), different mapping algorithms were applied to investigate the impact on the percentage of mapped reads and the impact on the percentage of reads mapped to annotated regions coding for miRNAs. Even though the percentage of reads mapped to the equine genome was increased to a median of 70% for Burrows-Wheeler Aligner (BWA) [61] (parameters -n 1 -o O -e O -k 0 -l 8 -t 4) and 95% for Bowtie2 [62] (parameters -q -very-sensitive-local), the percentage of reads mapping to known miRNA regions relative to the number of total reads was higher in case of Bowtie1 (4% of reads mapped to known miRNAs) when compared to Bowtie2 and BWA (2% mapped to known miRNAs). Thus, the mapper.pl script (implementing Bowtie1) was used to map the small RNA sequencing reads to the equine genome. After mapping, the reads were collapsed into 303000 unique reads per library on average and were used for novel miRNA identification and expression quantification with miRDeep2 [44]. Flow chart outlining the pipeline for novel miRNA detection as well as miRNA expression profiling and differential expression analysis. After extracting small RNA from serum samples of 79 horses, the samples were sequenced and analyzed using state of the art bioinformatics tools. While all samples were used for the novel miRNA identification, only 71 horses were used for the differential miRNA expression analysis.
Analysis of the read length distribution of all samples revealed a clear bimodality where on average 4.6% of the total reads constituted the first peak between 21 and 24 nucleotides (nts) while the second and more prominent peak was located between 29 and 33 nts with 95% of the total reads falling into this range ( Figure 2A). Since the percentage of mapped reads to the equine genome showed a rather low value of median 38% using Bowtie1 [46] (mapper.pl module), different mapping algorithms were applied to investigate the impact on the percentage of mapped reads and the impact on the percentage of reads mapped to annotated regions coding for miRNAs. Even though the percentage of reads mapped to the equine genome was increased to a median of 70% for Burrows-Wheeler Aligner (BWA) [61] (parameters -n 1 -o O -e O -k 0 -l 8 -t 4) and 95% for Bowtie2 [62] (parameters -q -very-sensitive-local), the percentage of reads mapping to known miRNA regions relative to the number of total reads was higher in case of Bowtie1 (4% of reads mapped to known miRNAs) when compared to Bowtie2 and BWA (2% mapped to known miRNAs). Thus, the mapper.pl script (implementing Bowtie1) was used to map the small RNA sequencing reads to the equine genome. After mapping, the reads were collapsed into 303,000 unique reads per library on average and were used for novel miRNA identification and expression quantification with miRDeep2 [44].

Novel miRNA Identification
The core algorithm of miRDeep2 reported a total of 721 putative novel miRNAs with a miRDeep2 score between 0 and 10 (Table S2). After applying our miRDeep2 score cut-off of 4 to filter for high confidence novel miRNAs, a set of 192 novel miRNAs remained. Of this subset of potential novel miRNAs, three mapped to ribosomal or transfer RNA regions (Rfam v. 12.2) [63] and showed non-significant randfold p-value suggesting the secondary structure was unlikely to match the one of a miRNA precursor. This resulting set of high confidence novel miRNAs is listed in Tables S3 and S4. The precursors of 47 of the novel miRNAs overlapped with a set of novel equine miRNA precursors identified in a previous study of our group [33].

Novel miRNA Identification
The core algorithm of miRDeep2 reported a total of 721 putative novel miRNAs with a miRDeep2 score between 0 and 10 (Table S2). After applying our miRDeep2 score cut-off of 4 to filter for high confidence novel miRNAs, a set of 192 novel miRNAs remained. Of this subset of potential novel miRNAs, three mapped to ribosomal or transfer RNA regions (Rfam v. 12.2) [63] and showed non-significant randfold p-value suggesting the secondary structure was unlikely to match the one of a miRNA precursor. This resulting set of high confidence novel miRNAs is listed in Tables S3 and S4. The precursors of 47 of the novel miRNAs overlapped with a set of novel equine miRNA precursors identified in a previous study of our group [33].

Threshold Read Count Value Determination by KS-Test
A threshold value of 13 (DESeq2 mean normalized) read counts was set for biological significance by applying the KS-test approach (see Materials and Methods for details). Therefore, the base mean counts for expressed miRNAs had to be ≥ 13 to be considered biologically relevant.

Threshold Read Count Value Determination by KS-Test
A threshold value of 13 (DESeq2 mean normalized) read counts was set for biological significance by applying the KS-test approach (see Materials and Methods for details). Therefore, the base mean counts for expressed miRNAs had to be ≥13 to be considered biologically relevant.

MiRNA Expression Profile
We detected 515 miRNAs (189 novel predicted miRNAs and 326 known miRNAs) to be expressed in our dataset. After filtering for low counts and mean threshold values of ≥13 normalized counts across all samples, 91 miRNAs remained, which were later used for differential expression analysis. Of the filtered miRNAs 61 (67%) were expressed at a low to moderate level (10-250 DESeq2 normalized counts on average), 16 miRNAs (18%) were expressed at a moderate to high level (251 to 999 DESeq2 normalized counts) while 14 miRNAs (15%) showed high expression levels of more than 1000 DESeq2 normalized counts. Of the total counts 86% were contributed by the top five most highly expressed miRNAs (eca-miR-486-5p, eca-miR-92a, eca-miR-191a, eca-miR-423-5p, eca-miR-148a). The miRNA with the highest expression counts was eca-miR-486-5p showing 63% of the total counts ( Figure 2B). A PCA plot was generated with the 50 most variable miRNAs. Even though no clear clusters of case and control horses were formed, one outlier sample was detected and removed. This outlier was more than three standard deviations away from the mean of the respective principal component loadings on principal component 1 and 3 and was therefore excluded from further analyses ( Figure A1). Thus, the resulting dataset for the following downstream analyses consisted of 71 samples.

Hemolysis Effect on miRNA Expression Profile
Twenty-six miRNAs were significantly affected by the level of hemolysis at a FDR threshold of 0.05 according to DESeq2 (Table S5). Figure A2 depicts two miRNAs representing expression levels of miRNAs that are either positively or negatively affected by the level of hemolysis. The analysis with the tool edgeR showed following differenced in contrast to the analysis with DESeq2: four miRNAs were not significantly affected by the level of hemolysis (eca-miR-744, eca-miR-128, eca-miR-28-3p and eca-miR-125a-5p) and additionally five significantly affected miRNAs were reported: eca-miR-423-5p, eca-let-7g, eca-miR-19b, eca-miR-425, eca-miR-7177b (Table S6). Hence the following linear model was used to determine differentially expressed miRNAs between asthmatic and control horses which accounts for the hemolysis effect and the underlying population structure: The two binary covariates family1 and family2 encode the underlying population structure of three cohorts: family1, family2 and an unrelated cohort.

Functional Enrichment of DEmiR Target Genes
For investigating functional enrichment of DEmiR target genes, we first constructed a highconfidence set of experimentally validated target genes of all the DEmiRs using the DIANA-TarBase database v. 7.0. This resulted in a set of 576 unique target genes. Table S8 illustrates all the target genes for each DEmiR. The top 10 most enriched GeneGO pathway maps associated with the 576 target genes regulated by DEmiRs are shown in Table 2.

Functional Enrichment of DEmiR Target Genes
For investigating functional enrichment of DEmiR target genes, we first constructed a highconfidence set of experimentally validated target genes of all the DEmiRs using the DIANA-TarBase database v. 7.0. This resulted in a set of 576 unique target genes. Table S8 illustrates all the target genes for each DEmiR. The top 10 most enriched GeneGO pathway maps associated with the 576 target genes regulated by DEmiRs are shown in Table 2.
Biological process networks in the GeneGo database analysis are networks of cellular functions and regulations based on functionally interconnected pathways. Process network analysis showed statistically significant enrichment for developmental networks (hedgehog signaling, epithelial-to-mesenchymal transition), signal transduction (NOTCH signaling, WNT signaling), translation (regulation of initiation), cell cycle (G1-S growth factor regulation, G1-S interleukin regulation, G1-S), cell adhesion (amyloid proteins) and proliferation (positive regulation cell proliferation). The enrichment for diseases resulted in an enrichment of respiratory tract related diseases (lung neoplasms, respiratory tract neoplasms, thoracic neoplasms, respiratory tract diseases, lung diseases, rectal neoplasms, genital neoplasms female, rectal disease, glioma, neoplasm neuroepithelial). We then investigated the regulatory miRNA-protein interaction networks using GeneGo database with homologous human miRNAs. This analysis showed nine of the eleven significant DEmiRs are part of an interconnected network containing the following highly interconnected hubs: ATF/CREB (activating transcription factor/cAMP response element modulator) family, polycomb repressive complex 1 (PRC1), Cyclin E, c-Fos and RelA (p65 NF-κB subunit) and histone deacetylase class I (Figure 4). Table 2. Top 10 enriched GeneGo pathway maps (based on MetaCore database). The pathway maps with the number of total genes involved and the number of target genes enriched for the pathway and the false discovery rate adjusted p-value (Adjusted p). Biological process networks in the GeneGo database analysis are networks of cellular functions and regulations based on functionally interconnected pathways. Process network analysis showed statistically significant enrichment for developmental networks (hedgehog signaling, epithelial-tomesenchymal transition), signal transduction (NOTCH signaling, WNT signaling), translation (regulation of initiation), cell cycle (G1-S growth factor regulation, G1-S interleukin regulation, G1-S), cell adhesion (amyloid proteins) and proliferation (positive regulation cell proliferation).
The enrichment for diseases resulted in an enrichment of respiratory tract related diseases (lung neoplasms, respiratory tract neoplasms, thoracic neoplasms, respiratory tract diseases, lung diseases, rectal neoplasms, genital neoplasms female, rectal disease, glioma, neoplasm neuroepithelial). We then investigated the regulatory miRNA-protein interaction networks using GeneGo database with homologous human miRNAs. This analysis showed nine of the eleven significant DEmiRs are part of an interconnected network containing the following highly interconnected hubs: ATF/CREB (activating transcription factor/cAMP response element modulator) family, polycomb repressive complex 1 (PRC1), Cyclin E, c-Fos and RelA (p65 NF-κB subunit) and histone deacetylase class I ( Figure 4).

Discussion
In the light of the constantly increasing number of allergies in humans [1] and also among companion animals [64], there is an urgent need for non-invasive biomarkers to help predict the risk of asthma development. We used a large cohort of 71 mature horses of different age, sex and country of origin, to search for a non-invasive biomarker for asthmatic condition and a better understanding of the pathology of the disease. Using the tool miRDeep2, we were able to identify 142 putative novel equine miRNAs in serum. The high-confidence novel miRNAs were combined with the known equine miRNAs from miRBase and expression profiling followed by differential expression analysis was performed. In agreement with previous reports of plasma expressed miRNAs in horses, one of the top expressed miRNAs across all samples was eca-miR-486-5p [65].

Hemolysis Effect on miRNA Expression Profile
MiRNAs that are affected by the level of hemolysis have been shown to be limited in their applicability as disease biomarkers [32,66] and in line with these results, we report 26 miRNAs with expression levels that are affected by the level of hemolysis (analysis performed with DESeq2). MiR-486 is known to be highly expressed in blood cells, mainly erythroid cells [67,68] and was one of the 26 miRNAs significantly affected by the level of hemolysis. From the four miRNAs commonly reported as showing a hemolysis-affected expression (miR-451, miR-16, miR-92a and miR-486p) [32,39,66] only eca-miR-486-5p (adjusted p = 2.5 × 10 −4 ) and eca-miR-16 (adjusted p = 0.011) showed a fold-change affected by the level of hemolysis. Some miRNAs show a decreased expression level with higher levels of hemolysis which can possibly be explained by increased presence of nucleases or other destabilizing agents released from lysed erythrocytes influencing miRNA turnover and thus leading to the degradation of certain miRNAs [69]. We accounted for hemolysis in our linear model, however lysis of other cell types could impact the results.

Asthma Related miRNAs
Differential miRNA expression analysis between asthmatic and control horses identified eleven statistically significant miRNAs ( Table 1). The differential expression of eca-miR-140-3p is considered to be of low confidence, since the signal was mainly driven by just one sample and the differential expression was not significant after the exclusion of the outlier individual ( Figure A4). However, since the specific outlier individual was not a global outlier (as shown by PCA) but rather only for this one miRNA, this individual was still kept for differential expression analysis. Thus, the remaining ten high-confidence DEmiRs between asthmatic and control horses were used for downstream literature research to investigate their possible biological involvement in the pathophysiology of severe equine asthma.

Biological Implications: Th2/Th17/Th1 Immune Response
The top differentially expressed miRNA (miR-128) has already been implicated in carcinomas [70,71] and serum miR-128 was suggested as potential biomarker for e.g., glioma [72]. A potential role of miR-128 in horse asthma may result from its pro-apoptotic properties [73,74]. This hypothesis is further supported by previous studies on mRNA in asthmatic horses that showed disturbances in the cell cycle [33,[75][76][77].
MiR-128 has recently been shown to be downregulated in asthmatic bronchial epithelial cells and to be part of a regulatory miRNA network that was confirmed to significantly increase the production of interleukin 6 (IL6) and interleukin 8 (IL8) [78]. The increased presence of the pro-inflammatory IL6 and IL8 were previously shown to be associated with the pathophysiology of asthma [79,80]. IL8 is a chemokine that is responsible for the recruitment and activation of neutrophil granulocytes and multiple studies showed its overexpression in bronchoalveolar lavage cells and airway epithelium of horses with severe equine asthma [81][82][83]. Consequently, airway neutrophilia is a prominent feature shown by horses suffering from severe equine asthma [13]. The pleiotropic cytokine IL6 is viewed as a marker of airway inflammation in asthma (showed in humans and animal models) and has been proposed as a therapeutic target in clinical trials [84,85].
The observed differences in the levels of several cytokines in equine asthma can most likely be explained by alterations in the CD4 + T cell development pathway. Downregulated miR-128 has been shown to positively regulate CD4 + T cell differentiation to T helper 2 (Th2) cells and to negatively regulate Th1 cell maturation. Increased levels of IL6 as well as TGFβ are known to positively regulate T cell maturation to T helper 17 (Th17) cells [26]. The differentially expressed miR-744 was shown to target TGFβ, thus its downregulation in horses suffering from severe asthma, as shown in this study, likely leads to increased levels of TGFβ [86].
Furthermore, increased levels of interleukin 17 (IL17A) in bronchoalveolar lavage cells of horses suffering from severe equine asthma were reported, suggest the increased presence of Th17 cells [87]. It is known that Th17 cytokines (IL17A, interleukin 17F (IL17F), and interleukin 22 (IL22)) lead to mucous cell metaplasia as well as increased levels of airway remodeling [88]. Also, differentially regulated miR-197 was shown to be involved in the reciprocal regulation of the IL6/STAT3 pathway [89]. Interestingly, the STAT3 transcription factor is essential for Th17 cell differentiation [90,91].
It has previously been stated that an exaggerated Th2 response as well as a Th17 response is able to explain a large portion of the pathophysiological events underlying severe equine asthma [13].
An enhanced Th17 response is also supported by a finding in a previous RNA sequencing study investigating the transcriptome of unstimulated and stimulated PBMCs collected from the same horses that this study covers. We showed a significant increase in C-X-C motif chemokine ligand 13 (CXCL13) transcript abundance as well as a decrease in interferon gamma (IFNG) expression in horses affected by severe equine asthma when compared to control horses [33]. CXCL13 is a B cell attracting chemokine that was shown to be predominantly produced by Th17 cells but not Th1 or Th2 cells [92]. Increased CXCL13 expression was linked with the formation of ectopic lymphoid structures like inducible bronchus-associated lymphoid tissue (iBALT) [93]. The formation of iBALT was recently shown to be dependent on Th2 as well as Th17 immunity in the course of a fungal infection of the lung in mice [93]. Additionally, iBALT was shown to be present in 90-100% of human asthmatic individuals and the abundance is correlating with asthma severity [94]. On the other hand IFNG is known to be the hallmark cytokine of Th1 cells and its downregulation further supports a decline of Th1 cell abundance [95].
An alternative hypothesis is that horses with severe asthma show dual positive Th2/Th17 cells which were recently discovered and were already reported to be present in an elevated number in the BAL fluid of individuals with severe asthma [96]. These striking concordances with this study further strengthen the hypothesis of a predominant Th2 and Th17 immune response in severe equine asthma.
Accordingly, we also hypothesize that the significant deregulation of the ten miRNAs, as shown in this study, might be a factor influencing the susceptibility of certain individuals to develop asthma due to a deregulation in the T cell maturation pathway leading to polarization of the immune response towards the Th2 and Th17 side and away from the Th1 side ( Figure 5).
Rather than focusing on anti-IL6 agents to control the IL6 pathway (or other cytokines) as a therapeutic approach for asthma and other disease as it has recently been proposed, a valid novel approach could therefore be to consider the reported DEmiRs as potential therapeutic targets. Following this approach, it might be able to prevent an amplified Th2/Th17 response and shift the balance back to an increased Th1 cell differentiation [84]. Additionally, IL6 production is increased by downregulated miR-128 while TGFβ production is enhanced by downregulated miR-744. We hypothesize that this leads to an increased Th2/Th17 immune response upon antigen challenge rendering certain horses susceptible to develop severe equine asthma. This hypothesis is supported by a previous finding of downregulated interferon gamma (IFNG) and upregulated C-X-C motif chemokine ligand 13 (CXCL13) in PBMCs of the same affected horses.

Biological Implications: Asthma and Cell Cycle Regulators
The closely related miRNAs miR-103 and miR-107a are known to pilot cell cycle arrest and their up-regulation in horses suffering from severe equine asthma support previous findings about cell cycle disturbances in equine asthma [33,97,98]. Compelling accordance with our results was provided by a recent study covering the autoimmune disease lupus erythematosus, a chronic inflammatory condition [99]. This study reported a deregulation of the cell cycle characterized by an upregulation of various miR-15/16 members, including miR-103 and miR-107a, as well as a downregulation of miR-744. Hence, upregulation of miR-103 and miR-107 as well as downregulation of miR-744 in asthmatic horses highlight a potential role of this miRNA network in chronic inflammatory conditions.
On the other hand, the significantly differentially expressed miR-361-3p is known to act as a tumor suppressor by interfering with the cell cycle. Interestingly this miRNA was recently shown to be significantly downregulated in patients affected by Sézary Syndrome, an aggressive CD4 + T-cell lymphoma [100]. Thus, the downregulation of eca-miR-361-3p in severe equine asthma might lead to the increased proliferation of CD4 + T-cells, or it might exhibit a yet unknown function in CD4 + T-cell differentiation.

Known Biological Implications of the Remaining DEmiRs
The epidermal growth factor receptor (EGFR)-mediated miRNA miR-7 is a known key player in multiple lung-related diseases and has previously been proposed as a biomarker in serum for chronic obstructive pulmonary disease and is thought to act by suppressing the coupling of SWI/SNF-related matrix-associated actin-dependent regulator of chromatin subfamily D member 1 (SMARCD1) with p53 [101,102]. MiR-148b-3p inhibits the expression of major histocompatibility complex (class I, G; HLA-G). HLA-G is a known susceptibility gene for asthma and therefore it has previously been proposed that miR-148b-3p might overtake a role in asthma susceptibility by interacting with HLA-G [103,104]. The differentially expressed miR-215 has been shown to target IL17RS and IL21 [105,106]. Finally, miR-30d has previously been proposed as biomarker for asthma severity and is known to lead to airway smooth muscle hypertrophy [66,107].
Since all 11 DEmiRs detected in the course of this study showed a small log2 fold change between case and control horses, these miRNAs might not act as a non-invasive diagnostic biomarker for equine asthma due to the low specificity of a potential test. However, the DEmiRs of this study still provide novel insights about the possible underlying pathophysiology of severe equine asthma and thus highlight candidate pathways to target in future therapeutic approaches. Figure 5. DEmiRs affecting CD4 + T cell development. Hypothetical network of DEmiRs affecting T helper (Th) cell maturation. Downregulated miR-128 is known to promote Th cell development towards the Th2 side while negatively regulating Th1 cell maturation. Additionally, IL6 production is increased by downregulated miR-128 while TGFβ production is enhanced by downregulated miR-744. We hypothesize that this leads to an increased Th2/Th17 immune response upon antigen challenge rendering certain horses susceptible to develop severe equine asthma. This hypothesis is supported by a previous finding of downregulated interferon gamma (IFNG) and upregulated C-X-C motif chemokine ligand 13 (CXCL13) in PBMCs of the same affected horses.

Biological Implications: Asthma and Cell Cycle Regulators
The closely related miRNAs miR-103 and miR-107a are known to pilot cell cycle arrest and their up-regulation in horses suffering from severe equine asthma support previous findings about cell cycle disturbances in equine asthma [33,97,98]. Compelling accordance with our results was provided by a recent study covering the autoimmune disease lupus erythematosus, a chronic inflammatory condition [99]. This study reported a deregulation of the cell cycle characterized by an upregulation of various miR-15/16 members, including miR-103 and miR-107a, as well as a downregulation of miR-744. Hence, upregulation of miR-103 and miR-107 as well as downregulation of miR-744 in asthmatic horses highlight a potential role of this miRNA network in chronic inflammatory conditions.
On the other hand, the significantly differentially expressed miR-361-3p is known to act as a tumor suppressor by interfering with the cell cycle. Interestingly this miRNA was recently shown to be significantly downregulated in patients affected by Sézary Syndrome, an aggressive CD4 + T-cell lymphoma [100]. Thus, the downregulation of eca-miR-361-3p in severe equine asthma might lead to the increased proliferation of CD4 + T-cells, or it might exhibit a yet unknown function in CD4 + T-cell differentiation.

Known Biological Implications of the Remaining DEmiRs
The epidermal growth factor receptor (EGFR)-mediated miRNA miR-7 is a known key player in multiple lung-related diseases and has previously been proposed as a biomarker in serum for chronic obstructive pulmonary disease and is thought to act by suppressing the coupling of SWI/SNF-related matrix-associated actin-dependent regulator of chromatin subfamily D member 1 (SMARCD1) with p53 [101,102]. MiR-148b-3p inhibits the expression of major histocompatibility complex (class I, G; HLA-G). HLA-G is a known susceptibility gene for asthma and therefore it has previously been proposed that miR-148b-3p might overtake a role in asthma susceptibility by interacting with HLA-G [103,104]. The differentially expressed miR-215 has been shown to target IL17RS and IL21 [105,106]. Finally, miR-30d has previously been proposed as biomarker for asthma severity and is known to lead to airway smooth muscle hypertrophy [66,107].
Since all 11 DEmiRs detected in the course of this study showed a small log2 fold change between case and control horses, these miRNAs might not act as a non-invasive diagnostic biomarker for equine asthma due to the low specificity of a potential test. However, the DEmiRs of this study still provide novel insights about the possible underlying pathophysiology of severe equine asthma and thus highlight candidate pathways to target in future therapeutic approaches.

mRNA-miRNA Interactions
For further in-silico downstream analysis we focused on the most significant DEmiR eca-miR-128. Potential target genes of miR-128 are involved in signal transduction, an indisputable part of immune response. Signal transduction was also the most enriched biological process by the differentially expressed Genes (DEGs) related to asthma in equine PBMCs [33].
All three up-regulated DEGs (RAB20, member RAS oncogene family (RAB20), bromodomain adjacent to zinc finger domain 2B (BAZ2B), H2.0 like homeobox (HLX)) in mock stimulated control and asthmatic PBMCs derived from the same blood sample collection as the serum samples used in this study are potential targets of miR-128 and were previously implicated in human asthma-like diseases [108,109]. A single nucleotide polymorphism (SNP) in RAB20 has been associated with childhood asthma in European-American and Hispanic-American populations using genome wide association study (GWAS) [110]. Using the same technique, an SNP within BAZ2B was associated with longitudinal changes in lung function and mean rates of decline by smoking pattern [111,112]. However, little is known about the BAZ2B function.
The third gene, HLX, is a Th1 specific transcription factor (TF) that interacts with another Th1 specific TF, T-box 21 (TBX21) required for the Th1 cells maturation and Th1-specific cytokine expression, including high expression of IFNG and repression of IL4 expression [113][114][115]. Variation in the genes of both the TFs have been associated with childhood asthma and variation within HLX significantly decreased its activity [116][117][118]. It is likely that an increased expression of HLX in asthmatic horses is maintained due to the decreased activity of the TF enforced by silencing of the HLX transcript by miR-128. Whereas childhood asthma is thought to be a Th2-related type of asthma [119], severe asthma in humans is characterized by high levels of Th1 cells [120]. Therefore, the decreased expression of miR-128 and increased expression of its target, the HLX serve as potential therapeutic targets for equine as well as human asthma.

Functional Enrichment of DEmiR Target Genes and Network Analysis
Functional pathway enrichment using MetaCore (Table 2) revealed that the most significantly enriched pathway map was the regulation of epithelial-to-mesenchymal transition (EMT) (adjusted p = 1.02 × 10 −4 ). EMT has been proposed to be a key player in airway remodeling in asthma [121,122]. The second most significant pathway map, PIP3 signaling, takes over a prominent function in airway inflammation. PIP3 is produced by PI3K and leads to the stimulation of several downstream targets including the Akt kinase. The PI3K pathway overtakes a major role in CD4 + T cell differentiation and activation, supporting our hypothesis of an altered CD4 + T cell development in severe equine asthma [123]. Blocking the PI3K/Akt signaling pathway has previously been proposed as a therapeutic approach for early stages of airway remodeling induced by the abnormal epithelial-to-mesenchymal transition [124].
Phosphorylation of EIF2 function pathway (adjusted p = 4.29 × 10 −4 ) is required for the activation of NF-κB which in turn leads to increased lung inflammation in response to an allergen challenge [125]. The regulation of GSK3β functional pathway (adjusted p = 0.29 × 10 −4 ) is essential for regulatory T cell function and has been proposed as a therapeutic target against allergic airway inflammation [126,127].
The biological process networks enrichment analysis resulted in a significant enrichment in the NOTCH signaling pathway, which is of special interest, since it was shown that the NOTCH pathway is crucial for Th17 cell differentiation [128].
Several of these most prominent hubs were previously implicated in asthma pathophysiology. The autoimmune disease modulating ATF/CREB pathway has been shown to promote Th17 differentiation together with CRTC2 [129]. Increased c-Fos expression in T lymphocytes has been shown to be involved in corticosteroid-resistant bronchial asthma [130]. The well-characterized NF-κB is a key regulator of adaptive and innate immune response that plays a pivotal role in allergic airway diseases [131]. Besides, the histone deacetylase class I hub confirms recent findings in human and murine model asthma research where histone deacetylases were reported to play an important role in asthma pathogenesis and histone deacetylate inhibitors showed promising results in asthma treatment studies in animal models [132]. PRC1 is responsible for gene silencing by post-translational modification of histones and exerts important functions in T cell differentiation. PRC1 recognizes H3K27me3 epigenetic modifications and condenses chromatin, thus leading to stabilized Th2 cell function and restricting the plasticity of the cells towards the Th1 side [133].

Conclusions
Using small RNA sequencing data from 71 individual horses, we identified 11 significantly differentially expressed miRNAs in the serum of asthmatic horses compared to controls. Several of the DEmiRs have previously been implicated in cell cycle control and some of them are known modulators of CD4 + cell differentiation and airway remodeling. This confirms a previous finding of our group reporting an impaired cell cycle control in RAO horses. We hypothesize that the immune response underlying the pathophysiology of severe equine asthma follows a Th2/Th17 driven manner rather than one of type Th1. This is backed by the fact that downregulated miR-128 was shown to negatively regulate T cell maturation towards Th1 but to also positively regulate the maturation towards Th2 cells. Additionally, a modulated cytokine profile towards the IL6 and TGFβ side caused by decreased levels of miR-128 and miR-197, as well as increased levels of miR-744 positively affect the maturation of T cells towards the Th17 side.
Therefore, we propose that the decreased levels of serum miR-128 might yield insights into the molecular mechanisms underlying asthma and might also provide room for novel therapeutic strategies. Rather than focusing on anti-interleukin agents for therapeutic approaches we propose that the identified differentially expressed miRNAs might act as novel therapeutic target to tackle severe equine asthma pathophysiology. However, to evaluate the potential of the DEmiRS to act as a non-invasive biomarker further research is needed. The 11 DEmiRs constitute future targets for future asthma research, helping to shed light on the still unknown molecular mechanisms underlying the disease. Since miRNAs are highly conserved and equine asthma shares striking similarities to human asthma, these new insights about the miRNA profile in the serum of asthmatic horses offers novel avenues for investigating the molecular pathology of human asthma [134,135].
Supplementary Materials: The following are available online at www.mdpi.com/2073-4425/8/12/383/s1. Supplementary file 1: Figure S1: Flowchart highlighting the changing sample size; Table S1: Table listing the horses of this study; Table S2: MiRDeep2 output table; Table S3: Identified putative novel equine miRNAs; Table S4: Identified putative novel equine miRNAs;  Table S7: Table outlining the statistics of the sample set used for differential miRNA expression analysis (n = 71); Figure S2: Bioanalyzer report. Supplementary file 2:        Figure A4. Expression of significantly differentially expressed miRNAs between asthmatic and control horses. Boxplots showing the DESeq2 normalized read counts for the significant DEmiRs. The respective FDR-adjusted p-value (adjusted p) is indicated. The first plot indicates the borderline significant eca-miR-215 followed by the boxplot for eca-miR-140-3p which clearly indicates that the signal for that miRNA is predominantly driven by one extreme outlier individual. After the exclusion of that outlier, eca-miR-140-3p did not show significant differential expression anymore (last plot; adjusted p = 0.106).  Figure A4. Expression of significantly differentially expressed miRNAs between asthmatic and control horses. Boxplots showing the DESeq2 normalized read counts for the significant DEmiRs. The respective FDR-adjusted p-value (adjusted p) is indicated. The first plot indicates the borderline significant eca-miR-215 followed by the boxplot for eca-miR-140-3p which clearly indicates that the signal for that miRNA is predominantly driven by one extreme outlier individual. After the exclusion of that outlier, eca-miR-140-3p did not show significant differential expression anymore (last plot; adjusted p = 0.106). Figure A4. Expression of significantly differentially expressed miRNAs between asthmatic and control horses. Boxplots showing the DESeq2 normalized read counts for the significant DEmiRs. The respective FDR-adjusted p-value (adjusted p) is indicated. The first plot indicates the borderline significant eca-miR-215 followed by the boxplot for eca-miR-140-3p which clearly indicates that the signal for that miRNA is predominantly driven by one extreme outlier individual. After the exclusion of that outlier, eca-miR-140-3p did not show significant differential expression anymore (last plot; adjusted p = 0.106).