Genome-wide identification of chicken bursae of Fabricius miRNAs in response to very virulent infectious bursal disease virus

Infectious bursal disease virus (IBDV) can cause a highly contagious immunosuppressive disease in young chickens. MicroRNAs (miRNAs) are crucial regulators of gene expression and are involved in the pathogenesis of IBDV infection. To investigate the roles of miRNA in chicken bursae of Fabricius in response to very virulent IBDV (vvIBDV) infection, RNA sequencing was performed to compare the small RNA libraries from uninfected and vvIBDV-infected group which was infected for 3 days. A total of 77 differentially expressed (DE) miRNAs were identified in BF, of which 42 DE miRNAs were upregulated and 35 DE miRNAs were downregulated. A gene ontology analysis showed that genes associated with cellular processes, cells, and binding were enriched. Moreover, pathway analyses suggested that apoptosis, T cell receptor signaling pathways, and chemokine signaling pathways may be activated following vvIBDV infection. In addition, we predicted the target genes of DE miRNAs and constructed an miRNA-mRNA regulatory network. In total, 189 pairs of miRNA-target genes were identified, comprising 67 DE miRNAs and 73 mRNAs. In this network, gga-miR-1684b-3p was identified with the highest fold change, as well as gga-miR-1788-3p and gga-miR-3530-5p showed a high degree of change. The above three miRNAs were considered to play vital roles in vvIBDV-host interactions. This study was the first to perform a comprehensive analysis of DE miRNAs in the bursa of Fabricius in response to vvIBDV infection, and it provided new insights into molecular mechanisms underlying vvIBDV infection and pathogenesis.


Introduction
Infectious bursal disease virus (IBDV), an important member of the family Birnaviridae, can cause infectious bursal disease which is highly contagious and immunosuppressive disease in chickens. There were two known IBDV serotypes, of which serotype-I strains include very virulent, classically virulent, and attenuated IBDV variants which show varying degrees of pathogenicity in chickens [1,2]. Very virulent IBDV (vvIBDV) can attack early B cells in the bursa of Fabricius (BF) [3,4], and cause depletion in the number of B lymphocytes and serious disruption or even necrosis of bursal tissue [5], and the release of various cytokines also caused severe damage to the BF [6]. Severe immunosuppression caused by vvIBDV increases the susceptibility of chickens to other pathogens such as those causing Marek's disease and Newcastle disease, which were responsible for serious economic losses to the global poultry industry. Therefore, exploring the underlying regulatory mechanisms is essential to prevent and control vvIBDV infections.
With the development of transcriptome analysis and high-throughput sequencing, the molecular mechanisms underlying host responses to virus infection regarding microRNAs (miRNAs) have been studied extensively [7][8][9][10]. miRNAs were endogenous small non-coding RNAs of approximately 20-25 nt that mainly regulate degradation or translation inhibition of mRNAs by binding to the 3′-untranslated regions [11,12]; moreover, they exert various functions associated with inflammation [13], cancer [14], and immune responses [15][16][17][18]. Accumulating evidence suggested that miRNAs were universal regulators of gene expression in animals, and plants and play a key role in the pathogenesis of virus infections. miRNAs of several avian viruses including avian leucosis virus subgroup J [19], Newcastle disease virus [20], infectious bronchitis virus [21], and duck enteritis virus [22] have been identified and examined by RNA sequencing on a genomic scale. Moreover, aberrantly expressed miRNAs in cells and tissues play a crucial role in virus-host interactions [23]. It has been reported that classical IBDV infection affected miRNA expression in DF-1 cells [24]. The discovery of miRNAs provides new insights into the mechanisms of gene regulation [25][26][27]. However, few reports were available regarding the effect of IBDV infection on chicken BF miRNAs, particularly with respect to vvIBDV infection. The objectives of this study were to assess genome-wide expression profiles of BF miRNA in response to vvIBDV strain LJ-5 infection using RNA sequencing; to screen and analyze differentially expressed (DE) miRNAs; and to reproduce a miRNA-mRNA regulatory network. In this study, to identify key miRNAs associated with IBDV infection, which would provide clues for better understanding the pathogenesis and immune mechanisms of how vvIBDV effects BF functions at an miRNA regulation level.

Experimental animals and sample collection
Three-week-old SPF chickens provided by the Harbin Institute of Veterinary Medicine (Harbin, China) were randomly assigned to a vvIBDV infection group (n = 30) and a control group (n = 30). The chickens in the IBDV infection group were inoculated with vvIBDV strain LJ-5 through eye-nose drops at a dosage of 10 3 ELD 50 /0.2 mL [28], and the chickens in the control group were treated with an equal volume of sterile phosphate-buffered saline (PBS). All chickens were housed in an animal facility under negative-pressure pathogen-free conditions and were provided with a standard diet and water. After a three-day period of infection, the chickens were euthanized under anesthesia as previously described [28]. BF tissues of each chicken were collected immediately following dissection and were separated into two subsamples, one of which was fixed in 4% paraformaldehyde and 2.5% glutaraldehyde solution, and the other was immediately frozen in liquid nitrogen and was stored at -80 °C. All instruments were treated with DEPC prior to use so as to inactivate RNases.

Histopathology and transmission electron microscopy (TEM)
Histopathology and TEM observation were performed as previously described [6,29]. Briefly, for histopathology observation, bursa tissues fixed in 4% paraformaldehyde over 24 h were processed by paraffin-wax cutting into 5-μm thick slices which were stained using hematoxylin and eosin (Beyotime Biotechnology, Shanghai, China) and were observed using a light microscope (Nikon E100; Nikon, Tokyo, Japan). For TEM observation, bursa tissues fixed with 2.5% glutaraldehyde were rinsed twice with PBS and were then fixed with 1% buffered osmium tetroxide. After this, the tissues were dehydrated with a graded alcohol series and were embedded in epoxy resin. Ultrathin sections were then stained with uranyl acetate and lead citrate for transmission electron microscope observation (H-7650; Hitachi, Tokyo, Japan).

miRNA sequencing
Three samples of each group were randomly selected for RNA sequencing. Total RNA was extracted with TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. RNA purity and integrity were examined using a NanoDrop 2000 device (Thermo Fisher Scientific, USA) and an Agilent 2100 device (Agilent, Santa Clara, CA, USA), respectively. Obtained total RNA was used to determine miRNA sequencing using Illumina HiSeq 2500 platform (Illumina, San Diego, USA). Raw sequencing reads were filtered to obtain clean sRNA sequences by removing (1) low quality reads, (2) reads lacking 3′ adapters but containing 5′ adapters, (3) reads without an insertion fragment, (4) reads containing poly-As, and (5) reads shorter than 18 nt. The remaining clean reads were aligned with sRNAs in the GenBank database (Release 209.0) and Rfam database (11.0) to identify and remove small nucleolar RNAs (snoRNAs), small nuclear RNA (snRNA), small cytosol RNA (scRNA), ribosomal RNA (rRNA), and transfer RNA (tRNA) reads. The remaining reads were then mapped to the reference genome Gallus gallus GRCg6a using TopHat2 (v2.1.1), and genes mapped to exons, introns, and repeat sequences were removed. The remaining clean reads were used for analysis.
The clean reads were subjected to a search using the chicken miRBase database (Release 21), were identified as existing miRNAs using Bowtie (v1.1.2), and were matched to the miRBase database of other animal species which were thought to be known miRNAs in chicken BF. Unannotated reads were mapped to the reference genome, and novel candidate miRNAs were identified based on their genome location and on hairpin structures as predicted using Mireap_v0.2 software. Only perfect alignments were considered novel miRNAs. miRNAs with a fold change ˃ 2 and a P-value < 0.05 in a comparison were considered significant DE miRNAs between the vvIBDV infection group and the control group.

KEGG and GO enrichment analysis of target genes
Bioinformatic analyses were used to predict DE miRNA targets as described previously [30]. Three software packages, including RNAhybrid (v2.1.2) [31], Miranda (v3.3a) [32], and TargetScan (v7.0) [33], were used to predict potential target genes of DE miRNAs. Results which were consistently identified were considered putative target genes of miRNAs. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases were used to analyze the biological functions of miRNA target genes [34]. Enrichment values of GO terms and KEGG pathways were produced using a P hypergeometric test, and P-values < 0.05 was considered to indicate significant enrichment.

qRT-PCR validation
To assess accuracy and reliability of DE miRNAs based on sequencing results, nine upregulated and nine downregulated miRNAs were selected for qRT-PCR confirmation. Primers were designed using Oligo 6 software (Table 1). Total RNA was extracted from IBDV-infected and uninfected BF tissues, and qRT-PCR was performed using a miRcute Enhanced miRNA Fluorescence Quantitative Kit (Tiangen, Beijing, China) and an ABI 7500 Real-Time PCR system (Applied Biosystems). U6 was used as an internal control for miRNA detection. The 2 −ΔΔCt method was used to calculate relative expression levels of target genes [35][36][37]. Three technical replicates were performed.

Statistical analyses
GraphPad Prism8 software was used to process qRT-PCR results, and data were statistically analyzed using a oneway ANOVA. All data were expressed as means ± standard error of the mean. Statistical significance was reported at P < 0.05.

Histopathology and TEM
Histopathological changes in the BF were shown in Fig. 1.
The control group showed intact tissue structures, a clear boundary between the BF cortex and medulla, and no pathological changes (Fig. 1A), whereas vvIBDV-infected BF tissues exhibited severe structural alterations. Massive infiltration of inflammatory cells, nuclear debris, necrosis, and exudation occurred in the follicles, in addition to widening of the interfollicular space (Fig. 1B). TEM results were shown in Fig. 1C and D. The ultrastructure of BF tissues of the control group was normal (Fig. 1C), whereas in vvIBDV-infected BF tissues, the number of lymphocytes was decreased, mitochondria were swollen, cristae were broken, and many vacuoles and damaged organelles occurred in the cytoplasm (Fig. 1D).

Deep sequencing analysis of small RNAs
To identify miRNAs involved in vvIBDV infection in chickens, six small RNA (sRNA) libraries were constructed. After sequencing and assembling, 48,284,949 and 42,663,699 raw reads were obtained from control and IBDV-infected BF tissues, respectively, ( Fig. S1; Table S1). Low-quality reads and adapter sequences were filtered, and 42,354,180 and 36,118,954 (83.67%-89.04%) clean reads were remained, respectively (Fig. S1, Table S1). A similar length distribution of sRNA sequences was observed in the six libraries (Fig. S2). sRNA sequences were mainly in the range of 18-25 nt (> 90%), of which 22 nt miRNAs had the highest compared with other miRNA lengths [38]. Clean reads of each library (including CK-1, CK-2, CK-3, LJ-1, LJ-2, and LJ-3) were aligned using the GenBank database, the Rfam database, and reference genome exons, introns, and repeat sequences, after which rRNA, scRNA, snoRNA, snRNA, tRNA, exons, introns, and repeat sequences were removed ( Fig. S3; Table S2). The remaining clean reads were used for further miRNA analysis.

Identification of existing, known, and novel miRNAs in BF tissues
In total, 1717 miRNAs were identified in the control and vvIBDV-infected groups; the numbers of miRNAs overlapping between the infection treatment and the controls were visualized in a Venn diagram, with 411 miRNAs belonging to the control group and 180 to the vvIBDV-infected group (Fig. 2). Among these, 1126 were expressed in both groups, and 571 existing miRNAs matched chicken miRNAs (Supporting information file 1), and 870 known miRNAs matched miRNAs from other animal species (Supporting

DE miRNA profiles
miRNAs were key factors regulating antiviral responses; therefore, genome-wide expression changes of BF miR-NAs in response to vvIBDV infection were studied. DE miRNAs associated with vvIBDV infection were screened between the control group and the vvIBDV infection group. We identified 77 DE miRNAs between the two groups, including 20 existing miRNAs, 36 known miR-NAs, and 21 novel miRNAs, according to the criteria of fold change ˃ 2 and a P-value < 0.05 (Supporting information file 4); 42 DE miRNAs were upregulated and 35 DE miRNAs were downregulated (Fig. 3). The top 20 miRNAs with the most significant differential expression were listed in Table 2. The identified novel miRNA sequences were listed in Table S3. These results confirmed that vvIBDV can significantly affect expression patterns of miRNAs in chicken BF tissues.

Functional annotation of DE genes
To assess the functions of DE miRNAs involved in vvIBDV infection, the target genes were annotated using GO and KEGG analyses. GO includes three ontologies: molecular functions, cellular components, and biological processes, and Fig. 4A showed that cellular processes, single-organism processes, and metabolic processes were significantly enriched in the molecular function category.
Additionally, the terms cell, cell part and organelle, and binding and catalytic activity were significantly enriched in the cellular components and molecular functions ontologies, respectively. Pathway-based analysis helps elucidate biological functions of genes. KEGG was used to identify significantly enriched pathways of target genes. Fig. 4B showed that 20 pathways were significantly enriched, including apoptosis, T cell receptor signaling pathways, and chemokine signaling pathways.

Target prediction and analysis of the miRNA-mRNA regulatory network
miRNAs play a role in gene regulation by directly silencing or indirectly reducing the expression of their target genes. To further investigate functions and potential regulatory roles of DE miRNAs, we constructed miRNA-mRNA co-expression networks (Fig. 5), and 189 pairs of miRNA-target genes were identified including 67 DE miRNAs and 73 DE mRNAs.

Confirmation of DE genes by quantitative reverse-transcription PCR (qRT-PCR)
To verify the reliability and accuracy of the sequencing results, 18 DE miRNAs (nine upregulated and nine downregulated genes) were selected and subjected to qRT-PCR. Fig. 6 showed that qRT-PCR results displayed similar trends as RNA sequencing results, suggesting validity of the sequencing data.

Discussion
IBDV is a pathogen which causes highly contagious and immunosuppressive disease, primarily in chicken BF tissues [39,40]. Transcriptional analysis regulation in BF tissues after IBDV infection can be used as a tool to obtain valuable insights into virus-host interactions and to elucidate the role of miRNAs in IBDV responses. Here, we reported for the first time miRNA alterations of vvIBDV-infected chicken BF tissues using RNA sequencing, showing that vvIBDV can affect miRNA expression profiles.
In the present study, vvIBDV-infected BF tissues were used for RNA sequencing, showing 1,710 miRNAs in the control and vvIBDV infection groups. Among them, 77 miRNAs, including 20 existing miRNAs, 36 known miR-NAs, and 21 novel miRNAs, were identified as DE miRNAs. These DE miRNAs were considered to play regulatory role during IBDV infection. Additionally, KEGG analysis predicted that the target genes of DE miRNAs were enriched regarding including apoptosis, T cell receptor signaling pathways, and chemokine signaling pathways, and miRNAs may be involved in regulating the expression of target genes associated with vvIBDV infection through these pathways.