MicroRNA transcriptome analysis of chicken embryo fibroblast cells infected with Newcastle disease virus variants

Abstract

Variations in the pathogenicity of Newcastle disease virus (NDV), the agent causing Newcastle disease, are associated with variants of different virulence. A few studies have characterized the expression of microRNAs (miRNAs) in NDV-infected avian cells. Here, the expression of miRNAs in chicken embryo fibroblasts (CEFs) infected with Herts/33 and LaSota NDV strains (highly virulent and nonvirulent, respectively) was determined using RNA sequencing. miRNAs involved in NDV infection included 562 previously documented and 184 novel miRNAs. miRNA target genes involved transcription factors, cell apoptosis, ubiquitin-mediated proteolysis, and protein processing in the endoplasmic reticulum. Potential target genes associated with autophagy were verified by qRT-PCR. No studies have documented the miRNA profiles of CEFs infected with NDVs variants. This study adds to our knowledge of the cellular miRNAs involved in NDV infection and the complex molecular mechanisms mediating virus-host interactions. The results of this study will aid the development of strategies against the chicken virus.

Introduction

Newcastle disease (ND) is a major disease that induces economic losses to the poultry industry (Lancaster 1976; Suarez et al. 2020). ND was first documented in Newcastle-upon-Tyne, England, and Java, Indonesia, in the mid-1920s (Dortmans et al. 2011). Infection with Newcastle disease virus (NDV) causes ND, which can be fatal in poultry. ND has been documented in various countries (Alexander 2000).

NDV belongs to Avulavirinae, a Paramyxoviridae subfamily (Alexander 2000). NDV has an RNA genome with the following structural proteins: fusion protein, phosphoprotein, hemagglutinin-neuraminidase with a large polymerase protein, matrix protein, and nucleoprotein (Chambers et al. 1986). NDV can infect more than 200 avian species through natural or experimental routes (Ganar et al. 2014). NDV strains can be lentogenic, mezogenic, and velogenic. Lentogenic and mesogenic NDV strains are often used to develop anti-NDV vaccines (Dortmans et al. 2011; Alexander 2000; Suarez et al. 2020).

MicroRNAs (miRNAs) are involved in cell proliferation, immune defense, cell differentiation, and tumorigenesis (Krol et al. 2010). miRNAs are often components of complexes with proteins (e.g., Argonaute and Dicer) and can form RNA-induced silencing complexes (Kim 2005). miRNA-mediated gene regulation in animals involves the degradation or inhibition of mRNA by targeting a small complementary sequence in mRNA (Bushati and Cohen 2007).

Cellular miRNAs regulate the replication and translation of viruses by either binding directly to the genomes of viruses or mediating changes within the transcriptomes of hosts (Trobaugh and Klimstra 2017). miR-451 was recently identified to modulate inflammatory cytokine expression during influenza infection (Rosenberger et al. 2012). miRNA-302a suppresses the induction of the cytokine storm and expression of influenza A virus-stimulated interferon regulatory factor-5 (Chen et al. 2017). In SARS-CoV-2 infection, virus-encoded miRNA prevents an innate immune response (Singh et al. 2021). In mosquito cells, small RNA autoregulates the replication of Dengue virus 2 DENV-vsRNA-5 (Hussain and Asgari 2014). Cellular miRNAs thus mediate the response of hosts to viral infections. However, the changes in miRNAs during NDV infection in the cells of chickens remain unclear.

Here, a miRNA map was constructed by RNA sequencing (RNA-seq), and the miRNAs associated with NDV infection in chicken embryo fibroblasts (CEFs) treated with the Herts/33 and LaSota NDV strains (which are highly virulent and nonvirulent, respectively) was identified. This study adds to our knowledge of the function of miRNAs in NDV infection and host-virus interactions. The information generated by this research will also facilitate the development of vaccines and enhance our understanding of host–NDV interactions.

Results

Coding transcript identification in NDV-infected CEFs

The infected cells and RNAs used in this study were the same as those in Liu et al. (2018). Eight million raw reads (average) were generated in each sequencing library using the Illumina HiSeq2500 platform. Clean reads accounted for approximately 98% raw reads. Most clean reads ranged in length from 20 to 24 nucleotides (peak at 22 nt). Reads of this size are typical for miRNAs obtained from Dicer-derived products. RNA sequences 18–35 nt in size were matched to the chicken genome sequence for distribution analysis. Approximately 90% of the sequences in each library were mapped to the chicken genome (Table 1).

Table 1 Statistics of small RNA library sequences

Identification of miRNAs related to NDV infection

The clean sequences were compared with public databases, and tRNA, snRNA, and rRNA sequences were removed. Totally 29,612 reads showed a complete match to known chicken miRNAs, including more than 3,000 unique sequences in each sample. In that, 562 mature miRNAs were identified, corresponding to 844 precursor sequences (Table S2). A total of 1,643 reads were predicted to be novel miRNAs, including approximately 200 unique sequences. All in all, 368 mature novel miRNAs, including 184 precursor sequences, were detected in the nine small RNA libraries using mirdeep2 and miREvo. Detailed information on these new miRNAs, including mature sequences and precursor sequences, is provided in Table S3.

Sequence similarity was used to group all miRNAs identified into specific families. A total of 224 miRNAs from 121 miRNA gene families were identified. Twelve families of miRNAs were only present in chicken, which might represent chicken-specific or NDV infection-specific miRNAs (Table S4). These miRNA families were compared miRNA families in 220 species, including model animal species such as Drosophila melanogaster, Caenorhabditis elegans, zebrafish, Xenopus tropicalis, chickens, mice, and humans.

miRNA expression profiles under NDV infection

Totally 562 previously documented and 184 novel miRNAs were identified in the nine NDV-infected samples and normalized using TPM values. A total of 27 downregulated and 12 upregulated miRNAs were identified in the Herts/33 vs. control comparison (Table S4). A total of 30 downregulated and 19 upregulated miRNAs were identified in the LaSota vs. control comparison (Table S5). Venn diagrams and volcanic plots were used to visualize differences in the abundance of miRNAs (Fig. 1).

Fig. 1

Differentially expressed miRNAs among samples. A Venn diagram of differentially expressed miRNAs between the two comparison groups (blank vs. LaSota; blank vs. Herts/33). B Hierarchical heatmap showing the converted expression values of miRNAs. Red: upregulation; blue: downregulation. B represents uninfected samples, H represents infected with Herts/33, and L represents infected with LaSota. C Volcano plot showing differentially expressed miRNAs between the two comparison groups (blank vs. LaSota and blank vs. Herts/33) (adjusted P < 0.05). miRNAs without significant differential expression are denoted by blue circles, and upregulated miRNAs are denoted by red circles, and downregulated miRNAs are denoted by orange circles

The number of miRNAs identified during the Herts/33 infection was similar to those identified during LaSota infection. Twenty two miRNAs were shared in both strains (Fig. 1A). Thus, miRNA expression profiles did not significantly differ among cells infected with different strains.

Functional analysis of differentially expressed miRNA targets

The functions of the miRNA targets were divided into three main categories: molecular function, cellular component, and biological process. Comparison of the miRNA targets of CEFs infected with each virus strain with controls revealed different gene function distributions. GO terms related to cellular components were significantly enriched in the Herts/33 vs. control comparison, in addition to binding and protein binding. GO terms related to metabolism were significantly enriched in the LaSota vs. control comparison, in addition to cellular part, cellular component, and organelles. GO terms related to molecular functions, including binding to proteins and ions, were significantly enriched (Fig. 2 and Tables S6 and S7).

Fig. 2

Gene Ontology (GO) analysis of differentially expressed miRNAs. A GO analysis of differentially expressed miRNAs in the Herts/33 vs. control comparison. B GO analyses of differentially expressed miRNAs in the LaSota vs. control comparison. The P-value was calculated using the hypergeometric distribution and corrected using the false discovery rate. Significantly enriched GO terms were determined according to an adjusted P < 0.05

The pathways related to the differentially expressed miRNA targets in Herts/33, and LaSota cells and control cells were determined using KEGG analysis to enhance our understanding of differences in differentially expressed miRNA targets between Herts/33, LaSota cells, and control cells. The functions of miRNA targets were mainly related to protein processing in the endoplasmic reticulum, transcription factors, cell apoptosis, and ubiquitin-mediated proteolysis (Fig. 3).

Fig. 3

KEGG annotations of differentially expressed miRNAs. KEGG analyses based on differentially expressed miRNAs in the A Herts/33 vs. control comparison and B LaSota vs. control comparison. The number of genes (circles); richness factor (colors). Significantly enriched KEGG pathways were determined according to adjusted P < 0.05

Validation of miRNA expression and target genes by qRT-PCR

The expression patterns of selected miRNAs were validated by stem-loop qRT-PCR. miRNAs that showed significantly downregulated expression during Herts/33 and LaSota infection included novel-483, gga-miR-133a-5p, and gga-miR-140-3p (P < 0.05*). The expression of gga-miR-1434 was significantly upregulated during both Herts/33 and LaSota infection (P < 0.05*), and gga-miR-6696-5p was significantly upregulated after Herts/33 infection (P < 0.05*). The Pearson correlation coefficients between the RNA-seq and qRT‒PCR data were high for both Herts/33 infection (0.8613) and LaSota infection (0.8138) (Fig. 4). gga-let-7b, gga-miR-103-3p, gga-miR-140-3p, and gga-miR-26a-3p are candidate binding sites of a circular RNA associated with NDV replication (data not shown).

Fig. 4

qRT-PCR analysis of differentially expressed miRNAs and target genes in CEF cells. A Expression patterns of differentially expressed miRNAs related to NDV infection. The expression of miRNA novel-483, gga-miR-133a-5p, and gga-miR-140-3p was significantly downregulated in CEFs infected with the Herts/33 and LaSota strains (P < 0.05*). The expression of gga-miR-1434 was significantly upregulated in CEFs infected with both the Herts/33 and LaSota strains (P < 0.05*). The expression of gga-miR-6696-5p was significantly upregulated following Herts/33 infection (P < 0.05*). B Expression patterns of target genes related to autophagy after NDV infection as determined by qPCR. The expression of ULK2 and RUBCN was significantly upregulated, and the expression of ATG3, DRAM2, RICTOR, ATG4a, ATG7 and ATG13 was significantly downregulated in CEFs infected with both the Herts/33 and LaSota strains (P < 0.05*). Data (mean ± SD) are from three samples. t test for gga-mir-6696-5p and one-way ANOVA for other miRNAs and target genes

The target genes of these differentially expressed miRNAs were predicted using miRDB. The target genes associated with autophagy (and predicted by miRDB) were verified using qRT‒PCR. The expression of ULK2 and RUBCN was significantly upregulated, and the expression of ATG3, DRAM2, RICTOR, ATG4a, ATG7 and ATG13 was significantly downregulated during Herts/33 and LaSota infection (P < 0.05*). The expression of gga-let-7b and gga-miR-103-3p was significantly downregulated after NDV infection (P < 0.05*). Thus, based on the general rules between the miRNA and target genes, we infer that ULK2 and RUBCN might be the target genes of gga-let-7b and gga-miR-103–3, respectively. ATG7 might be the target gene of gga-miR-6696-5p.

miRNA–mRNA network analysis

The expression of gga-miR-6696-5p and gga-miR-1434 was significantly upregulated after Herts/33 infection (|Log₂Foldchange|> 2; P < 0.01) and after LaSota infection (|Log₂Foldchange|> 2; P < 0.01), respectively. The expression of gga-miR-133a-5p was significantly downregulated after Herts/33 and LaSota infection (|Log₂Foldchange|> 2; P < 0.01). The miRNA–mRNA regulatory network is shown in Fig. 5.

Fig. 5

miRNA-mRNA interaction network associated with NDV infection in CEFs. Green circles indicate miRNAs. Red circles indicate target genes that have interacted with specific miRNAs. The green circle represents gga-miR-6696-5p. The yellow circle represents gga-miR-133a-5p, and the blue circle represents miRNA gga-miR-1434. The red circles represent different miRNA target genes

Discussion

Host miRNAs interfere with viral replication and infection, and often disrupts the expression of host miRNAs (Poltronieri, et al. 2015; Trobaugh and Klimstra 2017). Cellular miRNAs protect hosts from infection various viruses, such as the Nipah virus, hepatitis E virus, and respiratory syncytial virus (Foo et al. 2016; Haldipur et al. 2018; Othumpangat et al. 2012). However, the expression profiles of cellular miRNAs and the molecular mechanisms mediating the functions of miRNAs following NDV infection in chickens remain unclear. Here, we characterized the miRNA expression profiles of CEF cells infected with two NDV strains, the nonvirulent LaSota strain, and the highly virulent Herts/33 strain, using RNA-seq combined with bioinformatics analysis. In this study, we analyzed the abundance of miRNA in cells infected with two different NDV strains. Global miRNAs and target genes associated with NDV infection in CEFs were identified.

Previous studies have identified 73 and 64 differentially expressed miRNAs at 6 and 12 h after NDV infection using DF-1 cells and the JS 5/05 strain, respectively (Chen et al. 2019b) and 112 and 115 differentially expressed miRNAs at 24 and 48 h after NDV infection in chicken HD11 cells with LaSota strain infection, respectively (Mu et al. 2019). In our study, 39 and 49 differentially expressed miRNAs were identified following Herts/33 infection and LaSota infection, respectively. This variation in the number of differentially expressed miRNAs among studies might stem from differences in the NDV strains used or the multiplicity of infection of the infected cells.

Patterns of enrichment of the differentially expressed miRNA targets in biological processes differed in CEFs under Herts/33 infection and LaSota infection, suggesting that infections by different strains can affect biological processes. Several studies have indicated that NDV enters macrophage cells by caveolae-mediated endocytosis (Cantín et al. 2007; Sánchez-Felipe et al. 2014; Tan et al. 2018; Zhao et al. 2021). This study’s KEGG pathway analysis indicated that the miRNA targets were significantly enriched in endocytosis under Herts/33 infection and LaSota infection. These results suggest that NDV may infect cells through the endocytic pathway and that the endocytic pathway is independent of cell type. Virions constantly alter the actin cytoskeleton as they move through cells, from entry to exit. Viral infection may be recognized by cells detecting these cytoskeletal changes (Delorme-Axford and Coyne 2011). A focal adhesion localized tyrosine kinase called focal adhesion kinase (FAK) sends signals to the cytoplasm from the extracellular matrix (Bozym et al. 2012). FAK phosphorylation plays a pivotal role in herpes simplex virus entry, and FAK is involved in rabies virus infection through its interaction with viral phosphoprotein P (Cheshenko et al. 2005). FAK has been shown to play a critical role in the endosomal trafficking of influenza A viruses, and FAK kinase activity is important in promoting IAV replication by regulating the polymerase activity of multiple influenza A virus subtypes (Elbahesh et al. 2016). This study enriched miRNA targets in focal adhesion under Herts/33 infection and LaSota infection. This result suggests that FAK could play an important role in NDV entry. However, additional research is needed to investigate the mechanisms underlying these observations.

Upregulation of the expression of miR-181 in Marc-145 cells has been observed following infection by porcine reproductive and respiratory syndrome virus (PRRSV), and miR-181 inhibits PRRSV replication (Gao et al. 2013; Guo et al. 2013). miR-23 can induce the expression of type I interferon during PRRSV infection, suppressing PRRSV infection (Zhang et al. 2014). miR-491, miR-323, and miR-654 inhibit the replication of influenza A H1N1 virus by binding to PB1 (Song et al. 2010). miR-324-5p is downregulated in A549 cells under H5N1 and NDV infection and decreases H5N1 virus replication through its effect on CUEDC2 in the host and PB1 in the virus (Kumar et al. 2018). In this study, gga-miR-6696-5p and gga-miR-1434 expressions were upregulated following NDV infection. gga-miR-133a-5p mediates breast muscle development in chickens (Fu et al. 2018). gga-miR-133a-5p expression was significantly downregulated following Herts/33 and LaSota infection, suggesting that gga-miR-133a-5p mediates the response to NDV infection and important function during NDV infection. However, the specific mechanism by which gga-miR-133a-5p contributes to the response to NDV infection requires further investigation. miR-1756 and miR-6608 were detected only in chickens in this study. This finding suggests that these two miRNAs are chicken specific; the functions of these two miRNAs require clarification. miRNA expression profiles following NDV infection have been documented in various types of cells, including DF-1 cells, HeLa cells, and macrophage HD11 cells (Chen et al. 2019b, 2021; Mu et al. 2019). For the first time, our study uses primary cell CEFs, the natural hosts of NDV, to investigate NDV infection; however, the DF-1 cell line is a CEF-derived continuous cell line. Thus, the results of this study shed light on the interactions between NDVs and hosts. miRNAs associated with NDV replication have been identified in several studies. For example, both gga-miR-1794 and gga-miR-1603 directly target the viral L gene and inhibit NDV replication; the expression of gga-miR-375 is significantly upregulated after NDV infection; miR-375 markedly reduces NDV replication by acting on the NDV M gene and enhances the growth of NDV by decreasing the expression of ELAVL4; gga-miR-19b-3p inhibits NDV replication by targeting RNF11 and ZMYND11, which inhibits the inflammatory response; and NDV replication is inhibited by gga-miR-455-5p by targeting factors that suppress cytokine signaling (Chen et al. 2019a, 2021; Wang et al. 2019a, b). We characterized the miRNA expression profile in CEFs after NDV infection; however, the specific miRNAs involved in NDV replication and the mechanisms underlying these miRNAs remain unclear and require further investigation.

Autophagy is a cell survival pathway in development, immunity, and cell death (Deretic and Levine 2009; Levine et al. 2011; Maiuri et al. 2007; Yue et al. 2003). Long-lived proteins, protein aggregates, and organelles are among the cytoplasmic cargo that autophagosomes surround during the conserved cellular process of autophagy. These cargoes are then sent to lysosomes for destruction (Morishita and Mizushima 2019). Autophagy is mediated by evolutionarily conserved autophagy-related (ATG) genes (Levine and Kroemer 2019). ATG3 encodes a ubiquitin-like-conjugating enzyme and acts as an E2 ubiquitin-like conjugating enzyme in the ATG8 conjugation systems involved in autophagy, the process of degradation, turnover, and recycling of cytoplasmic constituents in eukaryotic cells (Fang et al. 2021). ATG5 and ATG12 are conjugated by the E1-like ligase ATG7, which is essential for the formation of autophagosomes. ATG7 is a necessary catalyst in both conjugation systems and is essential for autophagy (Tanida et al. 1999, 2001). In yeast, ATG13 is a crucial component of the ATG1 autophagy initiation complex, while ATG13 in mammals is required for the ULK1 complex to induce autophagy (Popelka and Klionsky 2017). Autophagy facilitates the degradation of intracellular components and can play a role in viral infection with either provirus or antiviral function depending on the virus and the stage of the viral replication cycle (Dong et al. 2021; Mao et al. 2019). Recent studies have suggested roles for autophagy in ZIKV replication, HPV, influenza virus, and HIV infection (Chiramel and Best 2018; Mattoscio et al. 2018; Mehrbod et al. 2019). NDV infection induces autophagy, which facilitates the replication of NDV in chickens (Sun et al. 2014). Several autophagy-related genes that might be target genes of miRNAs were identified, such as ATG3, ATG7 and ULK2. This result suggested that miRNA could regulate NDV infection through the autophagy pathway. However, the relationships between the target genes and miRNAs, and the specific mechanisms by which miRNAs regulate NDV replication through autophagy, require investigation.

Conclusions

In summary, we identified miRNAs and used deep sequencing methods to identify the cellular miRNAs in chicken CEF cells infected with the Herts/33 and LaSota NDV strains. Several previously documented and novel differentially expressed miRNAs related to NDV infection were identified. The potential miRNA target genes and relevant pathways were verified. Our findings reveal the miRNA expression profiles in CEFs infected with Herts/33 and LaSota. Additional studies of the role of specific miRNAs are needed to understand host-virus interactions during NDV infection.

Methods

Cells, viral propagation, viral infection, and RNA preparation

The primary CEF cell, virus, and RNA preparation protocols were the same as those described. The infection samples from Liu et al. were used to perform RNA-seq for both coding and noncoding RNA (Liu et al. 2018).

RNA sequencing

The program for constructing the transcriptome library was obtained from China Beijing Novogene Biotechnology Co., Ltd. Each sample consisted of 3 μg of RNA. The sequencing libraries were used to generate the NEBNext Multiplex Small RNA Library Prep Set for Illumina (NEB, Ipswich, MA, USA). Experiments were conducted in three replicates. Fifty-bp single-ended reads were generated using an Illumina HiSeq 2500/2000 platform.

miRNA identification and sequence analysis

The small RNA read length distribution was determined, and low-quality reads, and adapter sequences were removed. A specific set of lengths from the clean reads was used in the subsequent analysis. Bowtie was used to map the clean reads to the chicken genome (Langmead et al. 2009). Tags derived from repeat sequences, protein-coding genes, snRNA, tRNA, rRNA and snoRNA were removed.

We searched small RNA tags against the chicken miRNAs in miRBase (Release 21.0), and known miRNAs were identified. Potentially novel miRNAs were identified by aligning the unannotated small RNA (sRNA) sequences with the chicken genome. miREvo and mirdeep2 software, coupled with genomic sequences containing sRNA, were used to predict the folding energy and hairpin structures (Friedländer et al. 2012; Wen et al. 2012). Candidate novel miRNAs were identified as sequences with a typical hairpin structure with only a stem loop. Samples from other species were used to determine the presence of different miRNA families. Known miRNAs defined in miFam.dat (http://www.mirbase.org/ftp.shtml) were used to search for miRNA families; Rfam families were searched by submitting novel miRNA precursors to Rfam (http://rfam.sanger.ac.uk/search/).

Expression analysis

Levels of miRNA expression were expressed in transcripts per million (TPM) (Zhou et al. 2010). Normalized expression values were determined using the DESeq R package (V. 1.8.3) as follows (Wang et al. 2010): (number of reads × 1,000,000)/total number of miRNA reads in the library. The threshold for significance was an adjusted P-value of 0.05 (Benjamin-Hochberg method).

GO and KEGG analysis

The target gene candidates of differentially expressed miRNAs (referred to as “target gene candidates” in the following) were subjected to a Gene Ontology (GO) enrichment analysis. For GO enrichment analysis, a Wallenius noncentral hypergeometric distribution based on GOseq that could correct for gene length bias was used (Young et al. 2010). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was carried out using KOBAS (http://kobas.cbi.pku.edu.cn/) (Mao et al. 2005). A corrected P-value of 0.05 was required to consider enrichment significant.

Verification of differentially expressed miRNAs and target genes

Total RNA from cells infected with NDV strains was used for qPCR analysis. Uninfected cells served as a negative control, and the experiment was repeated with three biological replicates. After assessing the integrity and concentration of RNA, reverse transcription PCR was conducted. One reverse transcription reaction of one μg of total RNA was conducted using a mixture of eight miRNA primers to sequence miRNAs; another reverse transcription reaction was conducted using random primers to sequence chicken U6 and target genes. The procedure was run as follows: 37℃ for 10 min, 50℃ for 30 min, 85℃ for 5 min, followed by holding at 4℃. The expression of the eight miRNAs and target genes was validated using a Bio-Rad CFX96 Touch™ Real-Time PCR Detection System (Bio-Rad, CA, USA) and a standard SYBR Green PCR kit (Dongsheng Biotech, Guangzhou, Guangdong, China). miRNAs were verified using the stem-loop qPCR method. The miRNA primer sequences of miRNA, target genes, and mature miRNAs are shown in Table S1. Chicken U6 snRNA or beta-actin was an internal control for miRNA and target genes. Triplicate reactions were conducted as follows: 95℃ for 5 min, 95℃ for 5 s, and 59℃ for 10 s (40 cycles). Expression was quantified using the 2^−ΔΔCT method. Statistical analyses were conducted using the t-test for gga-mir-6696-5p and one-way ANOVA for other miRNAs and target genes. P-values < 0.05 were considered to be statistically significant. The correlation coefficient between the RNA-seq and qRT‒PCR results was calculated using GraphPad Prism 6 (V. 6.02, GraphPad Software, La Jolla, CA, USA).

miRNA-mRNA integrative network

We used the prediction tools at miRDB (http://mirdb.org/) for predicting the target mRNAs of differentially expressed miRNAs; upregulated or downregulated miRNAs were determined by |log₂Foldchange|> 2 and P < 0.01. The miRNA–mRNA regulatory network was constructed using different miRNAs and target genes through Cytoscape software (V. 3.7.1, https://cytoscape.org).