Distinct miRNA Profile of Cellular and Extracellular Vesicles Released from Chicken Tracheal Cells Following Avian Influenza Virus Infection

Innate responses provide the first line of defense against viral infections, including the influenza virus at mucosal surfaces. Communication and interaction between different host cells at the early stage of viral infections determine the quality and magnitude of immune responses against the invading virus. The release of membrane-encapsulated extracellular vesicles (EVs), from host cells, is defined as a refined system of cell-to-cell communication. EVs contain a diverse array of biomolecules, including microRNAs (miRNAs). We hypothesized that the activation of the tracheal cells with different stimuli impacts the cellular and EV miRNA profiles. Chicken tracheal rings were stimulated with polyI:C and LPS from Escherichia coli 026:B6 or infected with low pathogenic avian influenza virus H4N6. Subsequently, miRNAs were isolated from chicken tracheal cells or from EVs released from chicken tracheal cells. Differentially expressed (DE) miRNAs were identified in treated groups when compared to the control group. Our results demonstrated that there were 67 up-regulated miRNAs, 157 down-regulated miRNAs across all cellular and EV samples. In the next step, several genes or pathways targeted by DE miRNAs were predicted. Overall, this study presented a global miRNA expression profile in chicken tracheas in response to avian influenza viruses (AIV) and toll-like receptor (TLR) ligands. The results presented predicted the possible roles of some DE miRNAs in the induction of antiviral responses. The DE candidate miRNAs, including miR-146a, miR-146b, miR-205a, miR-205b and miR-449, can be investigated further for functional validation studies and to be used as novel prophylactic and therapeutic targets in tailoring or enhancing antiviral responses against AIV.

on previous studies in chickens) [2,36,38]. The control groups received complete FBS-free Medium 199. The final volume for each well was 500 µL. After 2 h of stimulation, tracheal rings were washed twice and incubated at 37 • C in fresh complete FBS-free Medium 199. The tracheal rings designated for cellular miRNA extraction (TOC samples), consisting of six biological replicates (six individual embryos) per treatment group, were collected 3-and 18-h post-infection or -stimulation and stored at −80 • C in RNAlater (Invitrogen, Burlington, ON, Canada). The tracheal supernatants of the rings that were designated for EV isolation and subsequent EV miRNA extraction (EV samples) were collected 24-h post-infection or -stimulation. For EV samples, there were three biological replicates (pools of five individual embryos per replicate) for each treatment group.

Extracellular Vesicle Isolation
All of the centrifugations and ultracentrifugations were performed at 4 • C. Tracheal supernatants from EV samples were centrifuged at 300× g for 10 min. to remove cellular debris. Subsequently, supernatants were recovered and centrifuged at 2000× g for 20 min. Supernatants were again recovered and ultracentrifuged at 10,000× g for 30 min. (Optima L-100XP, Beckman Coulter, Mississauga, ON, Canada). Supernatants were recovered and filtered with 0.2 µm filters (VWR, Montreal, QC, Canada). Supernatants were then ultracentrifuged at 100,000× g for 60 min. (Optima L-100XP, Beckman Coulter, Mississauga, ON, Canada). Supernatants were discarded and pellets were resuspended in FBS-free complete Medium 199 and ultracentrifuged at 100,000× g for a final 60 min. (Optima L-100XP, Beckman Coulter, Mississauga, ON, Canada). The supernatants were discarded, and the pellets were resuspended in 200 µL of phosphate-buffered saline (PBS, Gibco, Burlington, ON, Canada). Finally, for miRNA isolation from EV samples, 700 µL of QIAzol reagent (QIAGEN, Toronto, ON, Canada) were added to samples that were designated for RNA isolation and then stored at −80 • C. Samples designated for Western Blot and transmission electron microscopy analyses were stored at −80 • C in PBS.

Western Blot
Western Blot analysis was performed to confirm biomarkers compatible with EVs. The protein concentrations of the isolated EVs were determined using Micro BCA Protein Assay Kit according to manufacturer's instructions (Thermo Fisher Scientific, Burlington, ON, Canada). EV samples were lysed with RIPA Lysis Buffer (EMD Millipore, Oakville, ON, Canada) and incubated on ice for 5 min. The samples were then treated with NuPAGE LDS Sample Buffer (4X) containing 5% 2-mercaptoethanol (Invitrogen, Burlington, ON, Canada) and denatured for 5 min. at 95 • C. Equal amounts of proteins (about 30 µg) were separated by SDS-PAGE using pre-cast Mini-PROTEAN TGX electrophoresis gels (Bio-Rad, Montreal, QC, Canada). The separated proteins were transferred onto activated polyvinylidene difluoride (PVDF) membranes (Sigma-Aldrich, Oakville, ON, Canada). The membranes were washed in Tris-buffered saline with 0.05% Tween 20 (TBS-T), blocked in 5% skimmed milk TBS-T solution, and then washed again in TBS-T. The membranes were incubated with the primary antibody overnight at 4 • C, followed by the appropriate secondary antibody diluted in TBS-T for 1h at room temperature. Proteins were detected using Clarity Max enhanced chemiluminescence (ECL) Substrate (Bio-Rad, Montreal, QC, Canada), according to the manufacturer's instructions.

Negative Staining and Transmission Electron Microscopy (TEM)
TEM analysis was performed at the Facility for Electron Microscopy Research (FEMR), McGill University, to confirm and validate the purity of isolated EVs. Briefly, a 1:1 dilution of frozen EV samples was thawed and fixed with 2.5% glutaraldehyde in 0.1M sodium cacodylate buffer. The samples were then allowed to equilibrate at room temperature for 30 min. A 200-mesh copper TEM grid with carbon support film (Agar Scientific Ltd., Stansted, UK) was negatively glow discharged at 25 mA for 30 s (Electron Microscopy Sciences 100 Glow Discharge System, Hatfield, PA, USA) and loaded with a 5 µL droplet of sample for 3 min. Excess solution was carefully blotted off with Whatman Grade 1 filter paper and the grid was washed twice with a droplet of glycine for 1 and 2 min., respectively. The grid was then washed three times with a droplet of MilliQ water for 1 min. Finally, for negative staining, the grid was first floated on a droplet of 2% uranyl acetate (Electron Microscopy Sciences, Hatfield, PA, USA) to remove excess water and then on a second droplet for 1 min. Excess liquid was blotted off with Whatman Grade 1 filter paper. The samples were imaged by the FEI Tecnai G2 Spirit 120 kV TEM (Thermo Fisher Scientific, Hillsboro, OR, USA) equipped with a Gatan Ultrascan 4000 CCD camera Model 895 (Gatan, Inc., Warrendale, PA, USA). The micrographs were taken at appropriate magnifications to record the fine structure of EVs. The proprietary Digital Micrograph 16-bit images (DM3) were converted to unsigned 8-bit TIFF images.

MiRNA Isolation
Small RNAs were isolated while using the miRNeasy Mini Kit (QIAGEN, Toronto, ON, Canada) following the Quick-Start protocol of miRNeasy Mini Kit manufacturer's instructions. Briefly, tracheal rings from TOC samples (stored in RNAlater) were collected, lysed in 700 µL of QIAzol reagent (QIAGEN, Toronto, ON, Canada), and homogenized for two minutes using 0.5 mm glass beads (Biospec Products Inc., Bartlesville, OK, USA) and a tissue homogenizer (MP FastPrep-24 Classic Instrument, MP Biomedicals, Solon, OH, USA). The isolated EVs from EV samples were likewise lysed in 700 µL of QIAzol reagent (QIAGEN, Toronto, ON, Canada). The samples were then deproteinized in chloroform and centrifuged for 15 min. at 12,000× g at 4 • C. The upper aqueous portion of the samples was precipitated in 95% ethanol. The samples were sequentially washed and centrifuged with RWT and RPE buffers using RNeasy Mini columns in 2 mL collection tubes (QIAGEN, Toronto, ON, Canada). Finally, the purified RNA was eluted in 27 µL RNase-free water and quality control of RNA was performed using the RNA ScreenTape Analysis kit (Agilent Technologies, Santa Clara, CA, USA) and the Agilent 4200 TapeStation Analysis Software A.02.01 SR1 (Agilent Technologies, Santa Clara, CA, USA), according to the manufacturer's instruction. The samples were then stored at −80 • C.

Small RNA Library Preparation and Sequencing
MiRNA sequencing was performed at the Research Center of the CHU de Québec-Université Laval, Quebec, Canada. Twenty-four libraries for TOC samples (pooling of two replicates from each treatment group) to give three replicates per group and twelve libraries for EV samples to have three replicates within each group were prepared using the NEBNext Multiplex Small RNA Library Prep Set for Illumina (New England Biolabs, Ipswich, MA, USA), according to manufacturer's instructions. The total RNA from each sample was sequentially ligated to 3 and 5 small RNA adapters using T4 RNA ligase (New England Biolabs, Ipswich, MA, USA). Next, cDNAs were synthesized through reverse transcription using ProtoScript II Reverse Transcriptase (New England Biolabs, Ipswich, MA, USA) and amplified by PCR. Clean up and size selection of fragments was made using the Monarch PCR & DNA Cleanup Kit (5 µg) (New England Biolabs, Ipswich, MA, USA) and 6% polyacrylamide gel. Finally, the RNA libraries were sequenced (50 bp single-end) on an Illumina HiSeq 2500 platform (rapid mode) (Illumina, San Diego, CA, USA).

MiRNA Expression Analysis
Raw sequencing reads were processed by removing adapters and low-quality sequences using Trimmomatic, a read trimming tool for Illumina NGS data [39]. We used FastQC, a quality control tool for high throughput sequencing data, and MultiQC v1.7, a tool that creates a single report for different metrics and alignment statistics across many samples, in order to assess the quality of the generated reads [40,41]. We aligned the clean reads obtained from each library to mature miRNAs and to the reference genome of Gallus gallus in miRBase version 22.1 database to identify known miRNA expression levels in each group [42]. The MirDeep2 package, which maps reads against a library of known miRNAs from miRBase, was used for miRNA quantification from the reads, coordinated by a Nextflow workflow [43,44]. Small RNA sequencing libraries were normalized to counts per million (CPM) and evaluated for expression. MiRNA differential expression modeling and calculation were then done using R and R packages edgeR, tidyverse, magrittr and ComplexHeatmap [45][46][47][48][49]. Furthermore, the identified miRNAs for each treatment group were considered differentially expressed (DE) if their normalized expression fold changes (FC) relative to the control group was greater than or equal to 1.5-fold (log2FC ≥ 0.58 or log2FC ≤ −0.58) and if their False Discovery Rate (FDR) was less than 0.05 (FDR < 0.05). Venn diagram analysis for the DE miRNAs among the different treatment groups was performed using the online tool http://bioinformatics.pbs.ugent.be/webtools/Venn/. Furthermore, the distribution and intersection of the DE miRNAs were visualized using the UpSet software (http://vcglab.org/upset) [50].

In Silico Target Gene Prediction and Pathway Analysis
We downloaded targets of Gallus gallus' miRNAs from release 7.2 of TargetScan and version 6.0 of miRDB in order to perform functional enrichment analysis on the targets of DE miRNAs [51,52]. Next, we excluded targets with a low confidence score (below 99 'context++ score percentile' for TargetScan and below 95 'target prediction score' for miRDB), resulting in 4987 unique target genes for 982 miRNAs from TargetScan and 4887 unique target genes for 973 miRNAs from miRDB. Subsequently, we formed gene sets corresponding to targets of differentially expressed, up-regulated, and down-regulated miRNAs, according to TargetScan and miRDB (separately). We excluded genes that were targets of both up-regulated and down-regulated miRNAs in order to calculate the exact number of host target genes or host target pathways.
We used the KnowEnG analytical platform (https://knoweng.org/analyze) to perform pathway and gene ontology (GO) enrichment analysis on these gene sets [53,54]. KnowEnG is a computational system for analysis of 'omic' datasets in light of prior knowledge in the form of various biological networks. We used the KnowEnG's gene set characterization (GSC) pipeline in the standard mode (no knowledge network) and its knowledge-guided mode with the STRING co-expression network [55]. The knowledge-guided mode of this pipeline implements DRaWR, a method that utilizes random walk with restarts (RWR) in order to incorporate gene-level biological networks in the enrichment analysis to improve identification of important pathways and GO terms [56]. For the GSC pipeline, we did not use the bootstrapping option, selected Gallus gallus as 'species', and used default values for all other parameters. In addition, we mapped the Gallus gallus genes to Homo sapiens genes and used them in the same manner in the GSC pipeline with 'species' selected as Homo sapiens. The p-values of the Fisher's exact test corresponding to the standard mode of GSC pipeline were corrected for multiple hypotheses (Benjamini-Hochberg method) and FDR < 0.05 was considered to be significant. For the network-guided mode of GSC pipeline, the results were filtered based on the 'Difference Score' and those with a value larger than 0.5 were considered. The difference score is the normalized difference between the query probabilities and the baseline probabilities in the RWR algorithm, with the best score observed as one [56]. Finally, we used GraphPad Prism 8.4.3 for the illustration of the targeted pathways by the treatment group [57].

Chicken Tracheal Cells Release Extracellular Vesicles
While there is heterogeneity among EV cargo, depending on cell type and origin, the protein contents of extracellular vesicles and exosomes have been extensively studied [6].
The average protein concentration across all EV samples was 506 µg/mL. Based on the International Society for Extracellular Vesicles (ISEV)' recommendation, our results highlighted the presence of "EV-enriched" markers and the absence of non-exosomal markers in EV samples by Western Blot analysis [63]. Our results confirmed the presence of lysosomal LAMP1, a type I transmembrane protein, in both EV and cell lysate samples (Figure 1a). The LAMP1 protein bands identified in both samples were of slightly different sizes. This was expected, as the antibody detects a band of approximately 90-130 kDa. The variability in molecular weight is observed as a result of different levels of glycosylation of the target in different cell and tissue types [64]. In addition, endosomal TSG101, a cytosolic protein, and cytoskeletal β-actin proteins were detected in EV samples ( Figure 1a). Moreover, as expected, the EV samples lacked endoplasmic reticulum protein, GRP94 (Figure 1a) [63]. Finally, EVs were examined for morphological characteristics by transmission electron microscopy (TEM), which revealed vesicles of 50-150 nm in diameter with morphological characteristics of EVs ( Figure 1b). Biotechnology Information [58,59]. Two platforms, miRanda and RNAhybrid, were used to scan the eight segments of the AIV viral genome for potential target sites of the identified DE miRNA. The source code for miRanda, (written in C and downloaded from http://www.microrna.org/microrna/getDownloads.do) was used with default parameters for scaling parameter (4.0), strict 5′ seed pairing (off), gap-opening penalty (−4.0), and gap-extend penalty (−9.0), and adjusted parameters for score, set to greater than or equal to 160 (sc ≥ 160), and the minimum free energy (mfe), set to less than or equal to −16kcal/mol (en ≤ −16 kcal/mol), as previously suggested [60,61].

Chicken Tracheal Cells Release Extracellular Vesicles
While there is heterogeneity among EV cargo, depending on cell type and origin, the protein contents of extracellular vesicles and exosomes have been extensively studied [6].
The average protein concentration across all EV samples was 506 μg/mL. Based on the International Society for Extracellular Vesicles (ISEV)' recommendation, our results highlighted the presence of "EV-enriched" markers and the absence of non-exosomal markers in EV samples by Western Blot analysis [63]. Our results confirmed the presence of lysosomal LAMP1, a type I transmembrane protein, in both EV and cell lysate samples (Figure 1a). The LAMP1 protein bands identified in both samples were of slightly different sizes. This was expected, as the antibody detects a band of approximately 90-130 kDa. The variability in molecular weight is observed as a result of different levels of glycosylation of the target in different cell and tissue types [64]. In addition, endosomal TSG101, a cytosolic protein, and cytoskeletal β-actin proteins were detected in EV samples ( Figure 1a). Moreover, as expected, the EV samples lacked endoplasmic reticulum protein, GRP94 ( Figure 1a) [63]. Finally, EVs were examined for morphological characteristics by transmission electron microscopy (TEM), which revealed vesicles of 50-150 nm in diameter with morphological characteristics of EVs (Figure 1b).

Cellular and EV Treatment Groups Have Distinct miRNAs Expression Profiles
RNA quality control of the extracted RNA determined an average RNA integrity number (RIN) of 6.0 and an average concentration of 406 ng/µL. Following small-RNA sequencing, low-quality reads were filtered and adaptor sequences were trimmed. A total of 235,810,088 clean reads were obtained from TOC samples (three libraries per treatment group). A total of 27,298,247 clean reads were obtained From EV samples (three libraries per treatment group) (Table S1). Across all samples, the distribution of the small RNA sequencing length was consistent with the known length of mature miRNA, an average of 22 nucleotides [65]. The fragment lengths were primarily concentrated at 22 nucleotides ( Figure S1).
The miRNA expression profiles were evaluated in order to determine the ability of AIV infection, and LPS and polyI:C stimulation to influence the expression of cellular and EV miRNAs. A total of 692 known mature chicken miRNAs were detected. Following differential expression filtering (FC ≥ 1.5 and FDR < 0.05) and mean difference plot analysis, 228 DE unique miRNAs were identified ( Figures S2-S4). In addition, some DE miRNAs were present in several treatment groups ( Figure 2, Figure S5).  Within the TOC 3 h groups, a total of 18 unique DE miRNAs, of which two miRNAs were found in more than one group (Figure 3a-c). Our results demonstrated that 11, three, and two miRNAs were up-regulated, while four miRNAs were down-regulated in the TOC 3 h, AIV group (Tables 1 and 2). There were two miRNAs, gga-miR-1608 and gga-miR-6705-5p, which were DE in both the TOC 3 h AIV and TOC 3 h LPS groups ( Figure 2a, Table S2). Furthermore, within the TOC 18 h groups, our results showed a total of 88 unique DE miRNAs, of which six miRNAs were found in more than one group (Figure 4a-c). In the TOC 18 h groups, three, 32, and two miRNAs were up-regulated, and two, 15, and 40 miRNAs were down-regulated following treatment with AIV, LPS and polyI:C, respectively (Tables 1 and 2). There were three miRNAs, gga-miR-1451-5p, gga-miR-1563, and gga-miR-12234-5p, which were DE in both TOC 18 h AIV and TOC 18 h LPS groups and three miRNAs, gga-miR-6569-5p, gga-miR-12228-3p, and gga-miR-2184a-3p, which were DE in both the TOC 18 h LPS and TOC 18 polyI:C groups ( Figure 2b, Table S3). Table 1. Up-regulated miRNAs in TOC 3 h and TOC 18 h following treatment with avian influenza viruses (AIV), lipopolysaccharide (LPS), and polyI:C. Following differential expression filtering (FC ≥ 1.5 and FDR < 0.05), 11, three, and two miRNAs were found to be up-regulated in TOC 3 h AIV, LPS, and polyI:C groups, respectively. For TOC 18 h AIV, LPS, and polyI:C groups, three, 32, and two miRNAs were found to be up-regulated in each group, respectively.      With respect to the EV groups, a total of 145 unique DE miRNAs were identified, in which 59 miRNAs were found in more than one group (Figure 5a-c). Our results exhibited that 21, five, and 14 miRNAs were up-regulated in AIV, LPS, and polyI:C groups, respectively, while 57, 17, and 90 miRNAs were down-regulated in AIV, LPS, and polyI:C groups, respectively (Tables 3 and 4). There were 10 miRNAs, gga-miR-107-5p, gga-miR-1784b-5p, gga-miR-449b-5p, gga-miR-205b, gga-miR-210a-5p, gga-miR-1727, gga-miR-1464, gga-miR-6665-5p, gga-miR-7482-5p, and gga-miR-12284-3p, which were DE in all EV groups. In addition, there was one miRNA, gga-miR-383-5p, which was DE in both the EV AIV and EV LPS groups. Furthermore, 34 miRNAs were DE in both EV AIV and EV polyI:C groups. Moreover, four miRNAs, gga-miR-211, gga-miR-132b-5p, gga-miR-1597-5p, and gga-miR-12223-3p were DE in both EV LPS and EV polyI:C groups (Figure 2c, Table S4). Across all TOC 3 h, TOC 18 h, and EV groups, there were 67 up-regulated miRNAs, 157 down-regulated miRNAs, and four miRNAs, gga-miR12244-5p, gga-miR210a-5p, gga-miR30b-5p, and gga-miR-383-5p, which were found to be up-regulated in some groups and down-regulated in others. The miRNA gga-miR12244-5p was up-regulated in the TOC 18 h AIV group and down-regulated in the EV AIV group and gga-miR210a-5p was up-regulated in the TOC 3 h AIV and TOC 18 h LPS groups and down-regulated in all three EV groups. The miRNA gga-miR30b-5p was up-regulated in the TOC 18 h LPS group and down-regulated in the EV polyI:C group and, finally, gga-miR-383-5p was up-regulated in the EV AIV and EV polyI:C groups and down-regulated in the TOC 3 h AIV group (Tables 1-4). Table 3. Up-regulated miRNAs in EVs following treatment with AIV, LPS, and polyI:C. Following differential expression filtering (FC ≥ 1.5 and FDR < 0.05), 21, five, and 14 miRNAs were found to be up-regulated in EV AIV, LPS, and polyI:C groups, respectively.

Target Gene Prediction and Functional Annotation Reveals DE miRNAs Target Multiple Pathways
The miRDB and TargetScan databases were used to predict the possible gene set or pathway targets of DE miRNAs to characterize molecular and immunoregulatory functions of DE miRNAs presented in the current study. Target genes were predicted based on Homo sapiens (human), Mus musculus (mouse), and Gallus gallus (chicken) targets. However, there was limited information on Gallus gallus targets; therefore, the results based on the chicken database are provided in the Supplementary Materials (Tables S5-S9). The results presented in this study are based on human and mice databases. A total of 105, 34, and 151 projected genes and three, four, and 14 projected pathways were targeted by 15, three, and two DE cellular miRNAs in 3 h TOC AIV, LPS, and polyI:C groups, respectively. Moreover, 25, 578, and 384 projected genes and 17, 14, and 17 projected pathways were targeted by five, 47, and 42 DE miRNAs in 18 h TOC AIV, LPS, and polyI:C treatment groups, respectively (Tables 5-7). Finally, a total of 487, 147, and 817 genes and 16, 11, and 23 pathways were identified as targets for 78, 22, and 104 DE EV miRNAs following AIV, LPS, and polyI:C treatments (EV samples), respectively (Tables 5, 8, and 9). The targeted pathways related to immune responses and intracellular signaling are illustrated in Figure 6, while a complete list of targeted pathways and the respectively responsible miRNAs are provided in Tables 6-9.
The results presented here demonstrated that EV up-regulated miRNAs following AIV infection potentially could target gene sets related to mRNA processing and c-myc, while other gene sets or pathways could be targeted by EV down-regulated miRNA following AIV infection, such as the TGF-beta pathway and caspase cascade in apoptosis. Some gene sets, such as c-myc or TGF-beta signaling related genes related genes, can be regulated by both up-regulated or down-regulated EV miRNAs following AIV infection.
It was projected that some DE miRNAs, such as gga-miR-449b-5p or gga-miR-205b may target several pathways. For example, gga-miRNA-205b can potentially regulate mRNA processing, p53 transcription factor and Ras-related C3 botulinum toxin substrate 1 (RAC1) pathway (Table 9). Furthermore, gga-miR-12284-3p was also a miRNA down-regulated in all EV groups. Our data from the target gene set prediction within all three different treatments showed that this miRNA could regulate p53 transcription factor ( Table 9). The gga-miR-383-5p, which was down-regulated in the TOC 3 h treatment group, was predicted to target and potentially regulate the c-myc pathway (Table 8). In addition, this miRNA was up-regulated following infection with AIV and stimulation with LPS in the EV treatment groups and it may also target the c-myc pathway based on the EV data (Table 9).

The Functional Annotation Reveals DE miRNAs Target Multiple Segments of the AIV Viral Genome
We also investigated the target of miRNAs within the AIV viral genome in order to determine the possible roles of DE miRNAs in AIV infection. A total of 26 miRNAs were determined to target at least one segment of the viral genome (Table 10). One miRNA, gga-miR-1784b-5p, was found to potentially target two viral segments, segment 5 (NP protein) and segment 7 (M1 protein). Figure 7 illustrates the target sites within the viral segments. Our results demonstrated that some down-regulated EV miRNAs such as gga-miR-107-5p, gga-miR-6665-5p and gga-miR-1784b-5p, possibly can target AIV segments (Table 10). Among the DE miRNAs targeting the AIV viral genome in both the TOC and EV treatment groups, 4 miRNAs were up-regulated, and 22 miRNAs were down-regulated. In addition, the secondary structures for the miRNA-RNA interactions were predicted using RNAhybrid, which was also used to predict and confirm the mfe values predicted by miRanda ( Figure S6). MiRNA target sites within AIV viral genome. The target sites for miRNAs within the AIV viral genome were predicted using the miRanda and RNAhybrid algorithms. All of the segments were found to be targeted by at least one miRNA, except segment 8 (NS1, NEP). One miRNA, gga-miR-1784b-5p, was found to target both segments 5 (NP protein) and 7 (M1 protein). Specific target positions for each miRNA can be found in Table 10.

Discussion
Studies investigating cellular miRNAs following viral infections are essential for providing insight into the role of miRNAs in intracellular communication and the induction of antiviral responses. Identifying the mechanisms that can be regulated by miRNAs following infections or Segment 1: PB2

Discussion
Studies investigating cellular miRNAs following viral infections are essential for providing insight into the role of miRNAs in intracellular communication and the induction of antiviral responses. Identifying the mechanisms that can be regulated by miRNAs following infections or treatments will expand the current knowledge of host-pathogen interactions. This regulation is complex and, while host miRNAs can positively regulate antiviral responses, viruses have also been shown to impact miRNA expression to favor viral infection [66]. Here, we report for the first time that chicken tracheal cells secrete EVs. This work revealed that EV miRNAs have the potential to regulate antiviral responses in chickens. Furthermore, we observed that the miRNA profiles of chicken tracheal cells depend on the source (cellular versus EV) and the treatment (AIV, LPS, or polyI:C). These differences in the profile of miRNAs can be associated with active machinery of cargo loading during EV maturation and release. Previous studies have demonstrated that multiple mechanisms, including the endosomal sorting complex required for transport (ESCRT)-dependent and ESCRT-independent pathways, regulate the content of the EVs [67,68]. In addition, viral infections or stimulation via TLRs can change the EV cargo composition by targeting these pathways [8].
Additionally, we found a group of miRNAs that are in common among the different treatment groups in EV samples, while they are significantly affected by treatments. This suggests that, while the treatment groups have distinct miRNA expression profiles, certain miRNAs may have a fundamental role in which their loading into the EVs is independent of treatment or independent of the active cargo loading process.
EVs and their contents, including miRNAs, play critical roles in the regulation of the immune responses in the recipient cells by either activating or suppressing immune response genes [69,70]. In the context of influenza virus infections, previous studies demonstrated that miRNAs regulate antiviral responses by suppressing intracellular signaling pathways or activating pathways downstream of pattern recognition receptors (PRRs). For example, miR-92-5p, which was up-regulated in EVs following AIV infection, is able to enhance the activity of the NF-κB pathway [71]. Meanwhile, the expression of miR-449b-5p was suppressed following virus infection. The miR-449 family mainly interfere with influenza virus infection by enhancing type I IFNs [18]. Therefore, as highlighted in this study, the type of stimuli affected the EV cargo and, subsequently, can have potential impact on the recipient cells. It is likely to expect that the treatment of host cells with specific miRNAs in EVs that activate antiviral responses in the recipient cells may limit the replication of the virus in the neighboring cells. On the other hand, the virus may down-regulate the expression of some EV miRNAs and interfere with antiviral responses.
Network-guided gene set characterization analysis using KnowEnG's analytical platform enabled the identification of important GO terms and pathways, while incorporating co-expression relationships among targets of DE miRNAs [53]. These pathway analyses demonstrated that DE EV miRNAs can target several pathways that are critical during influenza virus infection, such as RAC1 and TGF-beta signaling pathways. The RAC1 pathway has been shown to enhance the replication of the influenza virus [72]. Therefore, targeting this pathway, through specific miRNAs, such as miR-205a, can be a strategy to limit the replication of the influenza virus either in the target cells or neighboring cells through EVs.
Analysis of the miRNA profile in chicken tracheal cells at different time points showed that the type of stimuli impacts the profile of DE miRNAs. For instance, the majority of the differentially expressed miRNAs were down-regulated in the polyI:C treatment group, while the majority of the differentially expressed miRNAs were up-regulated in AIV infection group at an early time point. Previous studies in chickens highlighted the possible role of miRNAs in the regulation of immune responses and the modulation of Marek's disease virus or avian influenza virus. As highlighted in the study by Wang et al., the profile of expressed miRNAs or miRNA-related mechanisms that are involved in the regulation host immune response could vary depending on the type of virus, cell stimuli, and the time point [33,73].
Moreover, as previously suggested, our analyses showed that multiple miRNAs can target the same gene and that, equally, a single miRNA can have multiple gene targets [74]. This suggests that it is potentially a combination of miRNA activities that modulate target gene expression. Among the cellular TOC samples, the distinct miRNA expression characteristics at different post-infection or -stimulation time points (3 h and 18 h) suggests that the expression of certain miRNAs is time-dependent. There were not any common DE miRNAs within each treatment group among different time points.
While the functions of most of the DE miRNAs identified in the various treatment groups remain unknown, some have been reported to be involved in the regulation of host immune responses and host-pathogen interactions. Previously, we showed that, in the same tissues, LPS is a potent TLR ligand for inducing pro-inflammatory cytokines, such as IL-6 and IL-1β cells, which limit the replication of AIV [1,2,36]. The current study showed that LPS up-regulated the expression of miR-146a in tracheal cells at a later time point. A previous study demonstrated that the activation of NF-κB signaling pathway leads to the up-regulation of miR-146a, which regulates and controls pro-inflammatory responses [75]. Moreover, miR-181 family, including miR-181a and -b, regulates pro-inflammatory responses by regulating NF-κB signaling pathway [76][77][78]. In this study, the expression of both miRNA-181a and -b was increased following LPS treatment. This result is in alignment with previous reports indicating increased expression of miR181b following LPS treatment in chicken macrophages [79]. Therefore, expressed miRNAs, such as miR-181 family or 146a following LPS treatment, can control pro-inflammatory responses through the NF-κB signaling pathway.
In addition, the results of the current study demonstrated that some miRNAs have the potential to target both innate and adaptive immune responses. For example, some DE miRNAs, such as gga-miR-205b and gga-miR-124a-5p, target p53 signaling pathway (Figure 6c,e, Tables 7,9). It has been demonstrated that p53 serves as a host antiviral factor by enhancing cytokine and antiviral gene responses in the lung, increasing the activity of dendritic cells (DC) and influenza-specific CD8+ T cells [80]. Therefore, these miRNAs could be employed to develop anti-influenza strategies and vaccine adjuvants for the control of AIV in chickens.
In addition, previous studies have also implied that miR-146 is involved in the regulation of immune responses to Salmonella infection in mice and it has direct targets within TLR signaling [81][82][83]. Our results indicated that gga-miR-12244-5p may target TLR signaling pathways. The miRNA gga-miR-146b-5p was previously suggested to be involved in the general process of inflammation during Avian pathogenic Escherichia coli (APEC) infection in chickens [28].
The DE miRNAs play a role in the regulation of antiviral responses, but they also might target the AIV genome. We were interested in miRNA target sites in the AIV genome, as recent studies have indicated that miRNAs may be capable of regulating virus translation and replication by directly binding to the virus genome [22][23][24][25][26][27]. While investigating the potential targets in the AIV viral genome, we determined that seven out of eight segments were targeted by at least one miRNA, suggesting that it is possibly a combination or cooperation of miRNA activity that can control viral activity. Our results demonstrated that both gga-miR-146a-5p and gga-miR-146b-5p, which are part of the same gene family, target segment 1 (PB2 protein) of the AIV genome. This suggests that, while gga-miR-146a-5p and gga-miR-146b-5p may play a role in modulating TLR signaling, they potentially have a dual function in regulating the PB1 protein. In addition, gga-miR-1784b-5p was the only miRNA identified that targeted more than one viral segment of AIV, both segment 5 (NP protein) and segment 7 (M1 protein). Our results suggest that this miRNA has multiple targets, both in the host and viral genomes. Several other miRNAs, gga-miR-129-5p, gga-miR-6641-5p, and gga-miR-1663-5p, target host mRNA processing and NF-kB signaling pathways, while also potentially targeting segments 2 (PB1 protein) and 3 (PA protein) of AIV.

Conclusions
In summary, we investigated the global miRNA expression profile in chicken tracheas in response to AIV infection and TLR ligand stimulation and the potential roles of specific DE miRNAs. Our results highlighted the possible role of EV contents in regulating antiviral responses. This accumulated evidence presented here elucidated the potential function of the candidate miRNAs on the host responses to AIV infection. We inferred that miRNAs play a major role in AIV infection and they can coordinate host defense against viral infections. In this study, the miRNA expression profiles were significantly regulated by the treatments with AIV, LPS, and polyI:C. These DE miRNAs were found to potentially target both host and viral genes. This study identified several miRNAs, such as miR-146a, miR-146b, miR-205a, miR205b, and miR-449, which could be employed in alternative strategies for the control of AIV in chickens, such as miRNA-based antiviral agents or vaccine adjuvants. Further functional studies are required to confirm the effect of treatment with candidate miRNA to inhibit AIV infection and induce stronger immune responses in in vitro and in vivo models.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2076-393X/8/3/438/s1, Figure S1: Sequence length distribution of miRNAs, Figure S2: Up-and down-regulated DE miRNAs in the (a) TOC 3 h AIV, (b) TOC 3 h LPS and (c) TOC 3 h polyI:C treatment groups, Figure S3: Up-and down-regulated DE miRNAs in the (a) TOC 18 h AIV, (b) TOC 18 h LPS and (c) TOC 18 h polyI:C treatment group, Figure S4: Up-and down-regulated DE miRNAs in the (a) EV AIV, (b) EV LPS and (c) EV polyI:C treatment groups, Figure S5: Visualization of intersecting sets of DE miRNAs in (a) TOC 3h groups treated with AIV, LPS and polyI:C, (b) TOC 18h groups treated with AIV, LPS and polyI:C and (c) EXO groups treated with AIV, LPS and polyI:C, Figure S6: Predicted secondary structures for miRNA targeting viral genome, Table S1: Normalized reads from small RNA sequencing, Table S2: Intersecting sets of DE miRNAs among TOC 3 h, Table S3: Intersecting sets of DE miRNAs among TOC 18 h, Table S4: Intersecting sets of DE miRNAs among EV groups, Table S5: Pathway analysis summary using the miRDB Gallus gallus database, Table S6: Prediction of pathways targeted by up-regulated cellular DE miRNAs using the miRDB Gallus gallus database, Table S7: Prediction of pathways targeted by down-regulated cellular DE miRNAs using the miRDB Gallus gallus database, Table S8: Prediction of pathways targeted by up-regulated EV DE miRNAs using the miRDB Gallus gallus database, Table S9: Prediction of pathways targeted by down-regulated EV DE miRNAs using the miRDB Gallus gallus database