Integrated analysis of miRNA and mRNA expression profiles in the bursa of Fabricius of specific pathogen-free chickens infected with avian reticuloendotheliosis virus strain SNV
Yubo Zhao, Qing Zhang, Meng Wang, Bingrong Wu, Saisai Zhao, Xinhui Wei, Youxiang Diao, Yi Tang, Jingdong Hu
TL;DR
This study explores how miRNA and mRNA expression changes in chickens infected with a virus, revealing key immune-related pathways involved in the infection process.
Contribution
The study identifies a dynamic miRNA-mRNA regulatory network and immune-related pathways in chickens infected with avian reticuloendotheliosis virus.
Findings
213 differentially expressed miRNAs and 3311 differentially expressed genes were identified in infected chickens.
1376 negatively correlated miRNA-mRNA pairs were found, with immune-related pathways being significantly enriched.
qRT-PCR confirmed key immune-related miRNAs and their target genes involved in the infection.
Abstract
Reticuloendotheliosis virus (REV) is a gamma retrovirus that can cause immunosuppression, dwarf syndrome and acute reticulocytoma in poultry. The molecular mechanism by which REV infection leads to immunosuppression and tumorigenesis is poorly understood. In this study, we elucidated the regulatory network of miRNA-mRNA and the major signaling pathways involved in REV-SNV infection. Therefore, we used the spleen necrosis virus (SNV) model of REV to inoculate one-day-old specific pathogen-free (SPF) chickens and then performed global miRNA and mRNA expression profiling by conducting high-throughput sequencing of 18 bursa of Fabricius samples collected at 7, 14, and 21 dpi. In total, 213 differentially expressed miRNAs (DEMs) and 3311 differentially expressed genes (DEGs) were identified. In the miRNA-mRNA network constructed based on the association analysis of these DEMs and DEGs, 1376…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsAnimal Virus Infections Studies · interferon and immune responses · Viral gastroenteritis research and epidemiology
Introduction
Avian reticuloendotheliosis virus (REV) is a gamma retrovirus that can cause immunosuppression, dwarf syndrome and acute reticular cell tumors in poultry (Wang et al., 2012; Witter and Glass, 1984). Spleen necrosis virus (SNV) is a type C retrovirus isolated from ducks. It can kill ducklings within one week of inoculation. Phylogenetically, SNVs are grouped among the REV subspecies. The REV strain was first isolated from turkey tumors in the United States in 1958 (Robinson and Twiehaus, 1974). REV is widely distributed in different regions and is reported in countries such as the United States, China, Australia, and Canada (MacDonald et al., 2019; Niewiadomska and Gifford, 2013; Singh et al., 2003). In flocks infected by REV, the vaccinations for Avian influenza, Newcastle disease, Fowl pox and Marek's disease have failed due to immunosuppression (Sun et al., 2009), which has caused significant economic losses in the poultry industry and may threaten human health (Cheng et al., 2007).
MicroRNAs (miRNAs) belong to a class of small noncoding RNAs with a size of 20–25 nucleotides. They play important roles in gene regulation by regulating post-transcriptional silencing of target genes (Lu and Rothenberg, 2018). They are major regulators of many cellular processes, including cell proliferation, differentiation, and apoptosis (Kargutkar et al., 2023). Changes in the expression of specific miRNA genes can lead to inflammation, cancer, and other diseases (Tsuchiya et al., 2006).
Infection by REV damages immune organs, including the spleen, bursa of Fabricius, and thymus. Our previous studies revealed that REV infection induces changes in the spleen transcriptome, including changes in metabolic pathways, immune response processes, and genes involved in tumorigenesis (Gao et al., 2019a; Gao et al., 2019b). The bursa of Fabricius is an important central immune organ of birds; it is the site where the target B cells are reached during virus infection. The REV-T strain can differentially express 88 miRNAs in the bursa of Fabricius (Yu et al., 2017). Therefore, we studied the bursa of Fabricius to assess the effects of REV-SNV infection on different immune organs.
Sequencing the whole transcriptome via next-generation sequencing (NGS) has become efficient, sensitive, and inexpensive over the years (Chang et al., 2011). Therefore, we used specific pathogen-free (SPF) chickens artificially infected with REV-SNV, collected bursa of Fabricius tissues from SPF chickens infected with REV for NGS, and conducted a comprehensive analysis of the overall changes in mRNA and miRNA expression. By screening the miRNA-mRNA interaction network to determine the functional relationship between miRNAs and mRNAs, our study might help elucidate the molecular mechanism by which REV infection leads to immunosuppression and tumor disease.
Materials and methods
Ethical statement
The animal tests involved were carried out in accordance with national and institutional guidelines. This study was conducted in strict accordance with the regulations of the Animal Protection and Utilization Committee of Shandong Agricultural University (No. SDAUA-2018-165).
Animal test design and sample collection
The REV strain SNV (GenBank: DQ003591.1) was provided by Professor Shuhong Sun and Zhizhong Cui of Shandong Agricultural University. The SPF chickens used in this study were purchased from Jinan SAIS Poultry Co. Ltd., China. SPF chickens (n = 40; one day old) were randomly divided into a control group (group C) and a REV-SNV-infected group (group I). The chickens in group I were intraperitoneally injected with 0.2 mL containing 3000 TCID50, and the chickens in group C were intraperitoneally injected with 0.2 mL of DMEM. All chickens were raised in negative pressure filtered air isolators. At 7, 14, and 21 dpi, five chickens in each group were randomly captured and euthanized, and the bursa of Fabricius was collected quickly and aseptically, immediately frozen with liquid nitrogen, and stored at -80 °C the following day.
Isolation and detection of total RNA
Total RNA was extracted from the bursa of Fabricius samples collected at 7, 14, and 21 dpi after inoculation using the TRIzol kit (Invitrogen, Carlsbad, CA, USA) following the instructions provided with the reagent. Five bursa of Fabricius samples from each group were mixed and divided into three parts for extraction. Total RNA integrity and gDNA contamination were evaluated using 1 % agarose gel electrophoresis, and no degradation or contamination of total RNA was detected. The concentration and purity of total RNA were measured using a spectrophotometer (DS 11, DeNovix, Wilmington, DE, USA). High-purity RNA samples were subsequently used (i.e., OD260/OD280 between 1.9–2.1, OD260/OD230≥2.0, and RNA integrity number ≥7).
Reverse transcription (RT)-polymerase chain reaction (RT-PCR) was performed to test whether SPF chickens were infected with REV-SNV, and the Oligo 7.0 software was used to design long terminal repeat (LTR) fragment primers. The sequences of primers used were as follows: 5′-TCGCTGATATCATTTCTCGGA-3′ LTR forward sequence and 5′-CAGCCAACACCACGAACAAA-3′ LTR reverse sequence. The PCR products were sequenced by TSINGKE Biological Technology (Beijing, China).
Analysis of sRNA sequencing data
Small RNA (sRNA) libraries were constructed using 2 µg of total RNA from each bursa of Fabricius sample. The sRNA libraries were sequenced by Gene Denovo Biotechnology Co. (Guangzhou, China) using Illumina HiSeq Xten on the Illumina sequencing platform. By screening the raw sequencing data, we removed the reads without inserted fragments and reads with inserted fragment lengths less than 18 nt. Reads with more than one base with a quality value lower than 20 in the data or those that contained N were removed, and reads that contained polyA were removed to obtain sRNA clean tag sequences for the subsequent analysis. Each sample was subsequently analyzed for clean reads with sequence lengths of 18–35 nt, and the selected clean reads were mapped using MiRdeep2 to the gallus reference genome (ftp://ftp.ensembl.org/pub/release-81/fasta/gallus_gallus/dna/). The tag sequences obtained by sequencing were annotated to noncoding RNAs from the GenBank and Rfam databases. Then, the Bowtie (version 1.1.2) software was used to compare the miRNA sequences of this species in the authoritative miRNA database miRbase. Finally, the content of existing miRNAs in the sample and the distribution of bases were obtained. Then, the MiRdeep2 software was used to predict the special secondary structure of the miRNAs and identify new miRNAs.
Differential miRNA expression analysis prediction of miRNA targets
The miRNAs obtained from each sample were summarized, and the expression level of the TPM (tags per million) of each miRNA is calculated (Zhou et al., 2010). The formula used was TPM = T*10^6^/N (T represents the number of miRNA tags, and N represents the total number of miRNA tags). The expression profiles of all miRNAs in all samples were obtained. The edgeR software was used to analyze differences in miRNA expression between group C and group I. The screening criteria for differentially expressed miRNAs were a more than two-fold change in expression and p < 0.05. The target genes of different miRNAs were predicted via Miranda (v3.3a) and TargetScan (version 7.0). The target genes of the differentially expressed miRNAs in each group were significantly enriched according to GO function and KEGG pathway, and the differentially expressed protein genes were identified.
Transcriptome sequencing and data analysis
The total RNA used for transcriptome sequencing was the same as the total RNA sample used for sRNA sequencing. The total RNA was removed with a Ribo-Zero^TM^ Magnetic Kit (Epicenter, Madison, WI, USA), the mRNA was enriched with Oligo (dT) magnetic beads, the mRNA was cut into short fragments with buffer solution, and random oligonucleotides were used as primers. Double-stranded cDNA was synthesized using a reverse transcriptase system. The double ends of the cDNA were purified and repaired, A base was added, the sequence was connected, and PCR amplification was performed. The library was sequenced by Gene Denovo Biotechnology Co. (Guangzhou, China) using the Illumina NovaSeq 6000 platform.
Reads were filtered using fastp (version 0.18.0) (Chen et al., 2018). Reads containing adapters were removed, reads containing more than 10 % unknown nucleotides (N) were removed, and low-quality reads containing more than 50 % low-quality (Q value≤20) bases were removed. Finally, high-quality clean reads were obtained. The reads were subsequently analyzed against the gallus reference genome using the HISAT2 software. The expression of each gene was calculated based on an FPKM (fragment per kilobase of transcript per million mapped reads) value. Genes with a false discovery rate (FDR) parameter < 0.05 and an absolute fold change ≥2 were considered to be differentially expressed genes.
Construction of the miRNA-mRNA interaction network
The correlation of expression between miRNAs and mRNAs was evaluated via DEGs and DEMs using the Pearson correlation coefficient (PCC). The negatively correlated miRNA-mRNA pairs were screened by a PCC < –0.7 and p < 0.05. A miRNA-mRNA interaction network was constructed, and the mRNAs in the network were analyzed for significant GO functional enrichment and significant KEGG pathway enrichment.
The expression of differential miRNAs and their target genes was verified by qRT-PCR
We tested the reliability of the high-throughput sequencing results via qRT-PCR detection of common differential miRNAs and their target gene expression levels at the three selected time points. A miRNA 1st Strand cDNA Synthesis Kit (Tail A) (Vazyme, Nanjing, China) was used to synthesize cDNA using the miRNA tailing method with U6 snRNA as the internal reference gene, and the primer sequences of the differential miRNAs are listed in Table 1. For qRT-PCR detection of differential target mRNAs, highly homogeneous cDNA was synthesized for qPCR quantification using a HiScript III 1st Strand cDNA Synthesis Kit (+gDNA wiper) (Vazyme, Nanjing, China). With GAPDH as the internal reference gene, the primer sequences of the target mRNAs are listed in Table 2. The relative levels of expression were normalized to the internal control. Real-time quantitative PCR analysis was performed using a LightCycler®96 (Roche Diagnostics GmbH, Germany) and ChamQ Blue Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China). The reactions were as follows: predenaturation at 95 °C for 30 s, followed by 40 cycles at 95 °C for 10 s and 60 °C for 30 s. Each cDNA sample was divided into three groups, and the qPCR data obtained were analyzed using the 2^–∆∆CT^ method to calculate the relative expression levels of the differentially expressed miRNAs and mRNAs. The primers were designed using the Oligo 7.0 software. All primers used in this study were synthesized by TSINGKE Biological Technology (Beijing, China).Table 1. The miRNAs primers used in this study.Table 1. PrimersSequence (5’–3’)Size (bp)gga-miR-215-5pATGACCTATGAATTGACAGAC21gga-miR-194TGTAACAGCAACTCCATGTGGA22gga-miR-1329-3pCCTCGTAGCTTGATCACGATAT22gga-miR-12224-5pAACAGTGGAATGCACCTAGCTGA23miR-125-xTCCCGGAGACCCTAACTTGTGT22miR-194-xTGTAACAGCAACTCCATGTGGACA24miR-12135-zTAAGGGTTCGTTTGTAAA18miR-215-xATGACCTATGAATTGACAGACTT23miR-194-zTGTAACAGCAACTCCATGTGGACC24miR-21-zATAGCTTATCAGACTGATGTTGACTTT27novel-m0035-5pTCCCTGAGACCCTAACTTGTGATGT25novel-m0124-3pAACAGGCACGACCGAGGCTGAAA23U6TGGAACGCTTCACGAATTTGCG22Table 2The mRNAs primers used in this study.Table 2. PrimersSequence (5’–3’)Size (bp)TLR3Forward: CAAGAGCCTGAGAACACTAGAGAT24Reverse: GCCAACTAGCAAACCAAGCAAT22CCL19Forward: GAACAGCAGCAGAACGCAGCAA22Reverse: ACCAGGAGACTAAGGCAGAGAACG24CCL1Forward: GTGCTGCTTCAACTTCCATACCA23Reverse: ATCTCCTTGCCATTCTTCACCTTG24CCL4Forward: GGCAGACTACTACGAGACCAACA23Reverse: TGATGAACACAACACCAGCATGAG24STAT1Forward: CCTCTGGAACGATGGCTGTATC22Reverse: TCCTGGTCTCTGGTCCTTCAATA23TRIM25Forward: TAACCACCACCCTCAGCGTTTC22Reverse: CAGATGCCAATGCCACAGAAGTT23C1RForward: TTATGGCTTCTACACCAAGATCCTC25Reverse: GGAACACTTAGAGCGTATCCTCC23C1SForward: CCAATTATCCTCAGGCATACCCAA24Reverse: CTTCACAGAGTCATATTCGCAGT23IFIH1Forward: CTTCCTCTACATGATCTCCTGCTTC25Reverse: CCGCACCTTGTCCTTCTCCT20IRF1Forward: CGTCCATCACCTTAGACCTCTCGT24Reverse: GATTTCCATCGGCATCCTCCAGTC24IFNGForward: AGCTCCCGATGAACGACTTG20Reverse: TCCTCTGAGACTGGCTCCTT20Vpreb2Forward: TCTCTGGGTATCGGCAGACC20Reverse: AGGTAAATTTCCAGTTGGATGGTG24YF5Forward: CGCTACTTCCTGACCGGGAT20Reverse: CCACGTACCCGACGGC16CDKN1AForward: GAGCGGTGGAACTTCGACTT20Reverse: CGGGCAAGTAGACCTTAGGC20CSF2RBForward: AAGGTGGCATTCATCCCAGG20Reverse: GTCCGCACAGACACTTGGTA20GADPHForward: GCCATCACAGCCACACAGAAGA22Reverse: GGCAGGTCAGGTCAACAACAGA22
Statistical analysis
All data were expressed as the mean ± standard deviation. Correlation analysis was performed between the NGS results and the qRT-PCR results using the GraphPad Prism 8.0 software.
Results
Verification of REV-SNV infection in SPF chickens
We performed RT-PCR analysis of the RNA extracted from the bursa of Fabricius samples collected from infected chickens and control group SPF chickens. For all samples in the infected group, PCR amplification bands corresponding to the 354 bp region of the TLR region were detected by conducting agarose gel electrophoresis (Fig. 1). After sequencing, we found that the PCR product sequences were identical to those of the REV-SNV (GenBank: DQ003591.1), which had the same LTR region; this finding confirmed that REV-SNV infection occurred in the infected group.Fig. 1. Detection of Reticuloendotheliosis virus-spleen necrosis virus (REV-SNV) RNA from the bursa of Fabricius of infected chickens. The DNA marker (M) used DL2000 (Takara, Dalian China). Lanes 1–3 are the PCR products of the IFN group at 7, 14, and 21 dpi, respectively. Lanes 4–6 are the PCR products of the CON group at 7, 14, and 21 dpi, respectively.Fig 1
mRNA-seq analysis of 18 bursa of Fabricius samples
We performed mRNA sequencing of 18 bursa of Fabricius samples, including six 7 dpi samples (con1-1, con1-2, con1-3, inf1-1, inf1-2, and inf1-3), six 14 dpi samples (con2-1, con2-2, con2-3, inf2-1, inf2-2, and inf2-3), and six 21 dpi samples (con3-1, con3-2, con3-3, inf3-1, inf3-2, and inf3-3). The sequencing data for each library were obtained (Table 3). According to the Q30 value of each library, the sequencing quality was satisfactory.Table 3. Summary of the mRNA sequencing data.Table 3. SampleRaw ReadsClean ReadsQ30 (%)Mapped Readscon1-139,627,35439,455,62091.6136,578,916con1-241,827,66041,644,66691.1738,622,712con1-339,583,01639,418,10091.4836,562,124con2-136,500,85036,333,76691.6733,251,063con2-240,097,82839,916,34091.9236,641,599con2-337,783,92437,601,55690.7234,610,802con3-141,622,85241,430,13491.7737,802,318con3-240,301,08040,096,48492.0336,178,414con3-336,973,05436,792,91291.6833,190,554inf1-136,944,85436,764,40491.4333,743,420inf1-240,045,67239,875,06291.5936,568,349inf1-339,562,71839,393,16891.9236,191,496inf2-139,395,67439,200,53490.5935,439,994inf2-240,186,51840,021,04891.9536,439,797inf2-346,801,53046,629,75494.3042,382,317inf3-140,118,62439,918,50890.9636,241,902inf3-240,162,77040,004,82092.4736,586,563inf3-339,918,22439,719,96491.4035,956,491
DEG analysis
The genes identified by an FDR < 0.05 and an absolute fold change ≥2 were considered to be DEGs. We identified 3311 DEGs in the CON and INF groups at 7, 14, and 21 dpi, and 144 DEGs at the intersection of the three time points (Fig. 2a). In the 7 dpi samples, 801 genes were identified, including 76 genes whose expression was upregulated and 725 genes whose expression was downregulated (Fig. 2b). In the samples collected at 14 dpi, 557 genes were identified, including 120 upregulated genes and 437 downregulated genes (Fig. 2c). Additionally, 2662 genes were identified in the 21 dpi samples, including 199 genes whose expression was upregulated and 73 downregulated (Fig. 2d). There were considerably more downregulated genes than upregulated genes at the three time points. Overall, the REV-SNV infection significantly affected the overall gene expression profile.Fig. 2. Venn diagrams. (a) Number and overlap of differential mRNAs during REV infection of the bursa of Fabricius at 7, 14, and 21 dpi. Venn diagram analysis revealed 144 mRNAs at all three tested time points. Volcano plot of differential mRNA expression at 7 (b), 14 (c), and 21 dpi (d). The x-axis represents log2 (FC), and the y-axis represents -log10 (P-Value). The red points represent the upregulated mRNAs, whereas the yellow points represent the downregulated mRNAs.Fig 2
miRNA-seq analysis of 18 bursa of Fabricius samples
We performed miRNA sequencing of 18 bursa of Fabricius samples, including six 7 dpi samples (micon1-1, micon1-2, micon1-3, miinf1-1, miinf1-2, and miinf1-3), six 14 dpi samples (micon2-1, micon2-2, micon2-3, miinf2-1, miinf2-2, and miinf2-3), and six 21 dpi samples (micon3-1, micon3-2, micon3-3, miinf3-1, miinf3-2, and miinf3-3). Sequencing data for each library were obtained (Table 4). According to the high-quality values of each library, the sequencing quality was satisfactory. Most tags were 21–24 nt long with only one peak and were most abundant at 22 nt (Fig. 3), which was similar to our previous results (Gao et al., 2019a). These results indicated that miRNAs were continuously enriched from 18 libraries.Table 4. Summary of the small RNA sequencing data after filtering and mapping.Table 4. SampleRaw readsClean readsHigh quality (%)Select ReadsSelect match readsMatch rate (%)micon1-112,526,39912,423,66798.1810,817,1358,680,79380.25micon1-210,569,8018,783,63483.107,776,6226,249,66380.36micon1-310,264,2968,521,32683.027,555,6926,062,76480.24micon2-112,510,97910,394,62383.088,901,0467,114,24179.93micon2-215,971,43415,765,47398.7112,709,2639,931,63378.14micon2-310,872,34810,775,94999.118,732,6366,899,96279.01micon3-111,234,58311,137,53199.148,691,3136,557,43775.45micon3-214,608,27614,515,87299.3711,269,4958,665,52976.89micon3-310,553,84910,443,30998.958,069,9846,132,83376.00miinf1-116,494,81016,302,76598.8413,708,36210,743,24778.37miinf1-29,660,5838,031,43783.147,024,4945,465,58477.81miinf1-310,354,73210,255,68199.048,714,9766,851,87878.62miinf2-110,691,05610,601,09399.169,663,6327,528,95277.91miinf2-210,285,6338,559,32183.227,925,2546,284,54879.30miinf2-39,881,8268,208,12483.067,447,7395,806,63477.97miinf3-19,439,3289,354,11499.107,827,9745,895,04475.31miinf3-216,767,43516,570,40598.8213,391,69810,148,61975.78miinf3-312,322,47712,210,47699.0910,401,4897,892,61675.88Fig. 3Length distribution results for the miRNA sequences. Only the micon1-1 example results are shown. The x-axis represents the tag length (nt), and the y-axis represents the tag count.Fig 3
DEM analysis
Using the edgeR software, normalization and statistical analysis were performed based on miRNA sequences. We identified 213 differentially expressed miRNAs (DEMs) in the CON and INF groups at 7, 14, and 21 dpi and 12 DEMs at the intersection of the three time points (Fig. 4a). Among the samples collected at 7 dpi, 69 DEMs were identified, including 35 upregulated and 34 downregulated miRNAs (Fig. 4b). Among the samples collected at 14 dpi, 105 were identified, including 34 upregulated and 71 downregulated miRNAs (Fig. 4c). Additionally, 108 miRNAs were identified in the 21 dpi samples, including 35 upregulated and 73 downregulated miRNAs (Fig. 4d).Fig. 4. Venn diagrams. (a) Number and overlap of differential miRNAs during REV infection of the bursa of Fabricius at 7, 14, and 21 dpi. Venn diagram analysis revealed 12 miRNAs at all three tested time points. Volcano plot of differential miRNA expression at 7 (b), 14 (c), and 21 dpi (d). The x-axis represents log2 (FC), and the y-axis represents -log10 (P-Value). The red points represent the upregulated miRNAs and yellow points represent the downregulated miRNAs.Fig 4
Functional enrichment analysis of differential miRNA expression
To determine the biological functions of different miRNAs in host chickens infected with REV-SNV, two methods, Miranda (v3.3a) and TargetScan (version 7.0), were used to predict candidate target genes. GO functional significance enrichment analysis was performed to identify the major biological functions of the target genes of the differentially expressed miRNAs. The results of the GO analysis (Fig. 5a-c) revealed that the genes associated with the GO biological process (GO-BP) category were enriched mainly in cellular terms “cellular process”, “metabolic process”, “regulation of biological process”, “response to stimulus”, “cellular component organization or biogenesis”, and “immune system process”. The GO cell component (GO-CC) category mainly included the terms “cell part”, “organelle part”, “membrane part”, “protein-containing complex”, “membrane-enclosed lumen”, and “extracellular region part”. The GO molecular function (GO-MF) category included the terms “binding”, “catalytic activity”, “molecular function”, “regulator, transcription regulator activity”, “transporter activity”, and “molecular transducer activity”.Fig. 5. Functional enrichment analysis of candidate target genes at 7 (a), 14 (b) and 21 dpi (c). The x-axis represents the second-order GO term, and the y-axis represents the number of miRNA target genes in this term. Different colors indicate different types of GO terms.Fig 5
The most significant biochemical metabolic pathways and signal transduction pathways of the miRNA candidate genes at 7, 14, and 21 dpi were determined by conducting the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. The results of the KEGG analysis (Fig. 6a-c) revealed that many of the genes enriched by the KEGG pathway were associated with Metabolic pathways, Viral protein interaction with cytokine and cytokine receptors, Lysosome, SNARE interactions in vesicular transport, Hippo signaling pathway-fly, Apoptosis and Protein processing in endoplasmic reticulum were significantly correlated. These results suggested that DEMs may play important roles in virus-host interactions when chickens are infected with REV-SNV.Fig. 6. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of candidate target genes at 7 (a), 14 (b), and 21 dpi (c). The first 20 pathways with the smallest Q values were mapped. The y-axis represents the pathway, and the x-axis represents the enrichment factor. The size represents the number of miRNA target genes; the redder the color, the smaller the Q value.Fig 6
miRNA-mRNA integrated analysis
The negative correlations and interaction networks of differential miRNAs and differential target mRNAs were comprehensively analyzed to provide valuable insights into the role of DEMs in the process of REV infection. Based on the 213 DEMs and 3311 DEGs screened in the previous steps, 1376 miRNA-mRNA interaction pairs were obtained, involving 84 DEMs and 610 DEGs, according to the criteria of negative correlation of expression levels and targeted relationships. The results revealed that most miRNAs have multiple target genes.
The functional enrichment analysis of miRNA-mRNA interactions with target mRNAs provided insights into their functional role during chicken infection with REV-SNV. The results of the GO analysis (Fig. 7) revealed that the GO biological process (GO-BP) category was enriched in mainly included the terms “cellular process”, “metabolic process”, “regulation of biological process”, “response to stimulus”, “cellular component organization or biogenesis” and “immune system process”. The GO cell component (GO-CC) category mainly included the terms “cell part”, “organelle”, “membrane part”, “macromolecular complex”, “membrane-enclosed lumen”, and “extracellular region part”. The GO molecular function (GO-MF) category included “binding”, “catalytic activity”, “catalytic activity”, “molecular function regulator”, “nucleic acid binding transcription factor activity”, “transporter activity”, and “transcription factor activity, protein binding”.Fig. 7. Functional enrichment analysis of target genes. The x-axis is the second-order GO term, and the y-axis is the number of miRNA target genes in this term. Different colors indicate different types of GO terms.Fig 7
The results of the KEGG analysis (Fig. 8) revealed that a large number of genes enriched by KEGG were associated with Viral protein interaction with cytokine and cytokine receptor, Proteoglycans in cancer, Cell adhesion molecules, Cytokine-cytokine receptor interaction, Chemokine signaling pathway, PD-L1 expression and PD-1 checkpoint pathway in cancer, NOD-like receptor signaling pathway, Rap1 signaling pathway, T cell receptor signaling pathway, Pathways in cancer, and JAK-STAT signaling pathway were significantly correlated. Using the Cytoscape network construction software, known miRNAs and mRNAs related to signal transduction, signaling molecules and interaction, transport and catabolic, cell growth and death, immune system and other aspects with negative correlations were visualized (Fig. 9).Fig. 8KEGG pathway analysis of target genes. The first 20 pathways with the smallest Q values were mapped. The y-axis represents the pathway, and the x-axis represents the enrichment factor (the number of miRNA target genes in this pathway divided by all the numbers in this pathway). The size represents the number of miRNA target genes; the redder the color, the smaller the Q value.Fig 8. Fig. 9Network diagram of miRNAs and negatively correlated target genes involved in immunity. The v shapes represent genes, and the circles represent miRNAs that target those genes.Fig 9
qRT-PCR validated significant DEMs and DEGs
Quantitative RT-PCR was performed to validate the differentially expressed miRNAs and mRNAs in the NGS data. The 12 DEMs at the intersection of 7, 14, and 21 dpi were selected (gga-miR-215-5p, gga-miR-194, gga-miR-1329-3p, gga-miR-1329-3p, miR-125-x, miR-125-x, miR-12135-z, miR-215-x, miR-194-z, and miR-21-z). The gene expression changes in the 12 DEMs were comparable among the three time points, with correlation coefficients of R^2^ = 0.9450 (p < 0.0001), 0.9037 (p < 0.0001), and 0.9176 (p < 0.0001) (Fig. 10a-c).Fig. 10. The DEMs and DEGs were confirmed by conducting by quantitative RT-PCR (qRT-PCR). (a-c) Correlation analyses of DEMs expression at 7, 14, and 21 dpi, respectively. (d-f) The level of expression of the DEGs at 7, 14, and 21 dpi, respectively. The x-axis represents the log2 (fold change) by qRT-PCR, and the y-axis represents the log2 (fold change) by NGS.Fig 10
The 15 DEGs at the intersection of 7, 14, and 21 dpi were selected (TLR3, CCL19, CCL1, CCL4, STAT1, TRIM25, C1R, C1S, IFIH1, IRF1, IFNG, Vpreb2, YF5, CDKN1A, CSF2RB). The gene expression changes of the 15 DEGs were comparable among the three time points, with correlation coefficients of R^2^ = 0.9125 (p < 0.0001), 0.9233 (p < 0.0001), and 0.9204 (p < 0.0001) (Fig. 10d-f). The results of qRT-PCR assays were consistent with those of high-throughput sequencing, which indicated that the sequencing data used in this study were reliable.
Discussion
High-throughput sequencing technology can be used for diagnosing diseases in plants and animals and for detailed studies of host-pathogen interactions at the genomic and transcriptome levels (Head et al., 2014; Van Borm et al., 2015). Through integrated analysis, the expression profiles of miRNAs and the genome in plant and animal tissues can be revealed to investigate the relationships between miRNAs and mRNAs, which has important reference significance for the disease and biology of animals and plants (Jia et al., 2024; Ortega-Acosta et al., 2024; Wu et al., 2020). REV-infected chickens attack mainly B-cell-mediated and T-cell-mediated immune organs and severely affect the immune functions of the spleen, bursa of Fabricius, and thymus (Xu et al., 2020). We previously studied changes in the miRNA and mRNA expression profiles of the spleen of SPF chickens infected with REV-SNV (Gao et al., 2019a; Gao et al., 2019b; Jiang et al., 2021). In chickens, the spleen is a peripheral immune organ, but the bursa of Fabricius is an important central immune organ and the main site of B lymphocyte development and differentiation (Zhang et al., 2023). To assess the effects of REV-SNV infection on different immune organs, in this study, we analyzed the expression profiles of miRNAs and mRNAs in the bursa of Fabricius of SPF chickens at 7, 14, and 21 dpi post REV-SNV infection to identify miRNAs and target genes that respond to REV-SNV. We constructed miRNA-mRNA interaction networks based on the screened DEMs and DEGs. The differential expression of 213 miRNAs and 3311 mRNAs in the bursa of Fabricius in REV-SNV-infected SPF chickens was identified, and 1376 negatively correlated miRNA-mRNA pairs were constructed, of which 82 pairs were identified at 7 dpi, 203 pairs were identified at 14 dpi, and 873 pairs were identified at 21 dpi. These results suggested that as the age of chickens increases, their immune system continues to develop and improve, and infected chickens exhibit increasingly strong immune protection responses against virus invasion. Some single miRNAs can target multiple genes. For example, gga-miR-214b-3p can simultaneously target PARD3 (par-3 family cell polarity regulator), TLR1A (toll-like receptor 1 family member), TRIM25 (tripartite motif containing 25), NOD1 (nucleotide binding oligomerization domain containing 1), and CD164 (CD164 molecule). A single gene can be targeted by multiple miRNAs, such as TNFSF8 (tumor necrosis factor superfamily member 8), which can be simultaneously targeted by gga-miR-2131-5p, gga-miR-301b-3p, and gga-miR-2131-3p.
By regulating the expression of target genes, miRNAs participate in various regulatory pathways, including those involved in viral replication, cell proliferation, differentiation, apoptosis, the antiviral immune response, and tumorigenesis (Fani et al., 2018; Gottwein and Cullen, 2008; Trobaugh and Klimstra, 2017). Many researchers have studied miRNAs at a greater depth and provided new insights into the study of host-virus interactions (Chen et al., 2020; Zhuo et al., 2013). Hongmei Li performed miRNA microarray analysis of the livers of 10-week-old chickens infected with avian leukosis virus (ALV) and identified 12 miRNAs with significant expression (Li et al., 2012). Mohammad Heidari performed miRNA sequencing on 21-day MD-resistant and MD-susceptible chickens in the bursa of Fabricius tissues infected with Marek's disease virus (MDV); they identified 54 new miRNAs and 10 known miRNAs. Several immune-related pathways, including the Toll-like receptor signaling pathway and other pathways, were identified through enrichment analysis (Heidari et al., 2020). To highlight the miRNAs associated with REV-SNV interactions with SPF chickens, we enriched mRNAs in the miRNA-mRNA negatively correlated networks to study their molecular mechanisms and biological functions. The analysis of the KEGG pathways of 610 mRNAs in the negatively correlated network of miRNA-mRNA revealed that the genes associated with this network were enriched mainly with Viral protein interaction with cytokine and cytokine receptor Proteoglycans in cancer and Cell adhesion molecules, Cytokine-cytokine receptor interaction, Chemokine signaling pathway, PD-L1 expression and PD-1 checkpoint pathway in cancer, Rap1 signaling pathway, T cell receptor signaling pathway, Pathways in cancer, JAK-STAT signaling pathway. Previous studies have reported that these pathways are involved in host antiviral processes (Goenka et al., 2023; Tzeng et al., 2021; Zheng, 2021). Overall, miRNAs influence the regulation of the host immune response, which may be related to the replication of REV-SNV during early infection in SPF chickens.
In the interaction network of miRNA-mRNA, TLR1A is negatively correlated with gga-miR-214b-3p and miR-1655-y. It is involved in the recognition of virus invasion by Toll-like receptor (TLR) and initiates the innate immune response through a series of signal transduction pathways (Bovenhuis et al., 2022; Kawai and Akira, 2007). CCR8 is negatively correlated with gga-miR-18b-5p, gga-miR-18a-5p, gga-miR-6542-3p, miR-2954-x, and miR-18-x and is involved in viral protein interaction with cytokine and cytokine receptor as chemokine receptors that are highly expressed on regulatory T cells (Tregs) (Suzuki et al., 2022). CCR8 often acts synergistically with other chemokine systems, as well as immune regulatory cytokines and growth and differentiation factors, to inhibit/promote tumors (Moser, 2022). NOD1 is negatively correlated with gga-miR-16c-5p, gga-miR-155, gga-miR-15c-5p, gga-miR-214b-3p, and miR-155-x, which participate in the NOD-like receptor signaling pathway. This pathway plays an important role in host innate immunity, regulating cell proliferation and apoptosis, autophagy, and the inflammatory response (Kuss-Duerkop and Keestra-Gounder, 2020; Trindade and Chen, 2020). CDKN1A is negatively correlated with gga-miR-15c-5p, gga-miR-12211-5p, and miR-2954-x. CDKN1A (p21 gene) is a member of the cell cycle protein-dependent kinase inhibitor family, which is involved in cell cycle regulation, cell migration, differentiation, and apoptosis. It is also associated with the occurrence and metastasis of tumors (Kreis et al., 2019; Xiao et al., 2020). EGFR is negatively correlated with gga-miR-16c-5p, gga-miR-155, gga-miR-20b-5p, miR-155-x, and miR-2995-x and is a key regulator of cell proliferation, growth, differentiation, and cancer development (Sabbah et al., 2020; Talukdar et al., 2020). EGFR plays a key role in regulating signaling processes during the cell cycle and is an essential component, considering that it can promote mitotic and carcinogenic effects (Lui and Grandis, 2002). LCK is negatively correlated with gga-miR-2188-5p, gga-miR-15c-5p, gga-miR-1456-5p, gga-miR-214-3p, and miR-347-z and plays a key role in T lymphocyte differentiation and activation. It can regulate the initiation of the T-cell receptor signaling pathway (TCR), T-cell development, and T-cell homeostasis (Bommhardt et al., 2019; Isakov and Biesinger, 2000). IRF1 is negatively correlated with gga-miR-106-5p, gga-miR-2131-5p, gga-miR-15c-5p, and gga-miR-12211-5p and plays an important role in protecting the host from pathogen invasion and maintaining homeostasis in the host, acting as the main regulator of innate immunity. Several studies have shown that IRF1 is associated with tumor suppression, the development of natural killer (NK) cells and T-cells, and B-cell biology (Feng et al., 2021; Mboko et al., 2015; Yamada et al., 1991).
To summarize, these results suggested that these miRNAs can regulate the function of target genes and promote or inhibit the replication of REV-SNV by targeting important regulators of immune-related signaling pathways.
Author contributions
Conceived and designed the experiments and reviewed the article: YXD YT JDH YBZ. Performed the experiments and wrote the paper: YBZ. Analyzed the data: YBZ and QZ. QZ, MW, BRW, SSZ, and XHW helped YBZ with animal husbandry and sample collection.
Declaration of competing interest
The authors declare that the research was conducted without any commercial or financial relationships that could be construed as potential conflicts of interest.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Bommhardt U.Schraven B.Simeoni L.Beyond TCR signaling: emerging functions of lck in cancer and immunotherapy Int. J. Mol. Sci.20201910.3390/ijms 20143500 PMC 667922831315298 · doi ↗ · pubmed ↗
- 2Bovenhuis H.Berghof T.V.L.Visker M.Arts J.A.J.Visscher J.van der Poel J.J.Parmentier H.K.Divergent selection for natural antibodies in poultry in the presence of a major gene Genet. Sel. Evol.542022243531379810.1186/s 12711-022-00715-9PMC 8939063 · doi ↗ · pubmed ↗
- 3Chang S.T.Sova P.Peng X.Weiss J.Law G.L.Palermo R.E.Katze M.G.Next-generation sequencing reveals HIV-1-mediated suppression of T cell activation and RNA processing and regulation of noncoding RNA expression in a CD 4+ T cell linem Bio 2201110.1128/m Bio.00134-11PMC 317562521933919 · doi ↗ · pubmed ↗
- 4Chen S.Zhou Y.Chen Y.Gu J.fastp: an ultra-fast all-in-one FASTQ preprocessor Bioinformatics 342018 i 884i 8903042308610.1093/bioinformatics/bty 560PMC 6129281 · doi ↗ · pubmed ↗
- 5Chen X.Ali Abdalla B.Li Z.Nie Q.Epigenetic regulation by non-coding RN As in the avian immune system Life (Basel)20201010.3390/life 10080148 PMC 745977932806547 · doi ↗ · pubmed ↗
- 6Cheng Z.Shi Y.Zhang L.Zhu G.Diao X.Cui Z.Occurrence of reticuloendotheliosis in Chinese partridge J. Vet. Med. Sci.692007129512981817602910.1292/jvms.69.1295 · doi ↗ · pubmed ↗
- 7Fani M.Zandi M.Rezayi M.Khodadad N.Langari H.Amiri I.The role of micro RN As in the viral infections Curr. Pharm. Des.242018465946673063658510.2174/1381612825666190110161034 · doi ↗ · pubmed ↗
- 8Feng H.Zhang Y.B.Gui J.F.Lemon S.M.Yamane D.Interferon regulatory factor 1 (IRF 1) and anti-pathogen innate immune responses P Lo S Pathog.172021 e 100922010.1371/journal.ppat.1009220 PMC 781961233476326 · doi ↗ · pubmed ↗