Intense Innate Immune Responses and Severe Metabolic Disorders in Chicken Embryonic Visceral Tissues Caused by Infection with Highly Virulent Newcastle Disease Virus Compared to the Avirulent Virus: A Bioinformatics Analysis

The highly virulent Newcastle disease virus (NDV) isolates typically result in severe systemic pathological changes and high mortality in Newcastle disease (ND) illness, whereas avirulent or low-virulence NDV strains can cause subclinical disease with no morbidity and even asymptomatic infections in birds. However, understanding the host’s innate immune responses to infection with either a highly virulent strain or an avirulent strain, and how this response may contribute to severe pathological damages and even mortality upon infection with the highly virulent strain, remain limited. Therefore, the differences in epigenetic and pathogenesis mechanisms between the highly virulent and avirulent strains were explored, by transcriptional profiling of chicken embryonic visceral tissues (CEVT), infected with either the highly virulent NA-1 strain or the avirulent vaccine LaSota strain using RNA-seq. In our current paper, severe systemic pathological changes and high mortality were only observed in chicken embryos infected with the highly virulent NA-1 strains, although the propagation of viruses exhibited no differences between NA-1 and LaSota. Furthermore, virulent NA-1 infection caused intense innate immune responses and severe metabolic disorders in chicken EVT at 36 h post-infection (hpi), instead of 24 hpi, based on the bioinformatics analysis results for the differentially expressed genes (DEGs) between NA-1 and LaSota groups. Notably, an acute hyperinflammatory response, characterized by upregulated inflammatory cytokines, an uncontrolled host immune defense with dysregulated innate immune response-related signaling pathways, as well as severe metabolic disorders with the reorganization of host–cell metabolism were involved in the host defense response to the CEVT infected with the highly virulent NA-1 strain compared to the avirulent vaccine LaSota strain. Taken together, these results indicate that not only the host’s uncontrolled immune response itself, but also the metabolic disorders with viruses hijacking host cell metabolism, may contribute to the pathogenesis of the highly virulent strain in ovo.


Introduction
Newcastle disease virus (NDV), also known as avian orthoavulavirus 1 (AOAV-1) or avian paramyxovirus 1 (APMV-1), is an enveloped single-stranded negative-sense RNA virus of birds from the genus Orthoavulavirus belonging to the subfamily Avulavirinae of family Paramyxoviridae, with a nonsegmented genome of 15 kb containing coding sequences for six genes [1][2][3]. The virus has infected at least 236 bird species tested to date, including domestic and free-living species [1]. Isolates are of a single serotype, but have a broad virulence from avirulent (lentogenic) to intermediate virulent (mesogenic), and highly virulent (velogenic), depending on the severity of the disease that they cause in birds [1]. Highly virulent viruses that cause diseases with severe morbidity are termed Newcastle disease (ND), which poses a major threat to many species of birds and causes severe economic losses for commercial poultry and pet birds [4]. On the contrary, avirulent viruses cause asymptomatic infections in adult birds or mild respiratory suffering and are widely used as live ND vaccines in commercial poultry [1]. Therefore, vaccines with live, avirulent viruses are the most popular and accepted control and prevention approaches for combating ND in the poultry industry worldwide [5]. However, research on the host's innate immune responses to infection with either a highly virulent strain or an avirulent strain, and how this response may contribute to severe pathological damage (and even mortality) in response to the infection of a highly virulent strain, remains limited.
Specific-pathogen-free (SPF) chicken embryos are widely used as an ideal and broad experimental model for the infection and propagation of multiple viruses, such as NDV, infectious bronchitis virus (IBV), and avian influenza virus (AIV), because of their accessibility, cost-effectiveness, and fast growth [6][7][8][9][10]. Meanwhile, numerous studies suggested that the immune system in birds starts to grow early during embryogenesis, as well as the fact that numerous immune responses are generated in the latter stage of chicken embryos [11,12]. Furthermore, severe inflammatory pathological lesions and high viral titers were presented in 9-to 10-day-old SPF CEVT when the embryonated chicken eggs (ECE) were challenged with various viruses [13][14][15]. Therefore, to better explore the exact differences in epigenetic and pathogenesis mechanisms between the highly virulent and avirulent strains, we established an SPF chicken embryo model infected with either a highly virulent strain NA-1 (a predominant genotype VII NDV strain presented in parts of most of Asia, including China) or an avirulent widely used vaccine strain LaSota (a genotype II NDV strain) for 24 and 36 h post-infection (hpi), and performed RNA-sequencing (RNA-seq) on viral infected CEVT, to obtain their transcriptional profiles. These findings will not only bring on an understanding of differences in infectivity and pathogenesis between highly virulent and avirulent NDV strains, but also be helpful to novel vaccine development and other control strategies.

Embryonated Eggs and Viruses
The SPF chicken embryos were purchased from Jinan SAIS Poultry Company (Shandong, China). The eggs were incubated at 37.9 • C, with 60% ± 5% humidity for 10 days, and then were infected with viruses or mock. Genotype VII NDV virulent strain NA-1 (GenBank No. DQ659677) [16] and genotype II NDV vaccine strain LaSota (GenBank No. AF077761) were propagated in the allantoic cavity of 9-10-day-old SPF ECE and purified directly from the allantoic fluid; then, the HA titers of allantoic fluid were measured by HA standard assays [17] and stored at −80 • C until used again [18].

Virus Infection, Sample Collection and Preparation
In the current work, 88 10-day-old chicken embryos were randomly divided into three groups, including NA-1 (n = 32), LaSota (n = 32), and mock (n = 24). The allantoic cavity of eggs was injected with either NA-1 or LaSota at a single dose of 7 × 10 3 TCID50 per egg, and eggs were injected with equal volumes of phosphate-buffered saline (PBS) as the mock-infected group. The inoculated chicken embryos continued to incubate at 37.9 • C and were candled every 6 h to detect dead embryos. Every eight chicken embryos from the viral infected group and six chicken embryos from the mock group were randomly selected at 12, 24, 36 and 48 hpi; then, the allantoic fluid and visceral tissues were collected and stored at −80 • C after quick freezing in liquid nitrogen. The hemagglutination assay tests were performed on allantoic fluid harvested from the ECE [17], and their visceral tissues were sent for RNA-seq.

RNA Extraction, cDNA Library Construction and High Throughput Sequencing
Total RNA was extracted using Total RNA Extractor (Trizol) according to the manufacturer's instructions (Sangon Biotech, Shanghai, China). Subsequently, RNA concentration was detected by Qubit2.0 (ThermoFisher, Waltham, MA, USA; Carlsbad, CA, USA), RNA purity was evaluated by NanoDrop 2000 (ThermoFisher, Waltham, MA, USA; Carlsbad, CA, USA), and the agarose gel electrophoresis was used to detect RNA integrity and genome contamination, and then qualified RNA samples were selected for cDNA library construction.
For RNA-seq, the mRNA in total RNA was enriched by specific binding of oligo (dT) magnetic beads. Then, the mRNA with poly (A) was eluted from the magnetic beads and fragmented into 200-300 bps fragments. Additionally, cDNA was synthesized by reverse transcription of the enriched mRNA fragments with random primers. After adding the adapter to the two segments of the DNA fragment, the obtained fragment was enriched by PCR to obtain the cDNA library. RNA-seq was performed by Sangon Biotech (Shanghai, China) based on an Illumina Hiseq system (Illumina, San Diego, CA, USA).

Bioinformatic Analysis of Sequencing Data
Transcriptome assembly and annotation protocols were provided by Sangon Biotech in China. At first, low-quality reads, readings containing adaptors, and reads containing poly-N > 10% were removed from the raw data according to the manual of Sangon Biotech, then the sequences of clean reads were aligned with the Gallus gallus reference genome using HISAT2 software, followed by gene expression levels quantified as transcripts per million (TPM) values. The differentially expressed (DE) genes (DEGs) in the groups between NA-1 and LaSota were identified by DESeq2, when corrected to p < 0.01 and |log2 (Fold change, FC)| > 1. DEGs were subjected to functional enrichment analysis using clusterProfiler, including gene ontology (GO) enrichment analysis [19] and Kyoto Encyclopedia of Genes Genomes (KEGG) pathway analysis [20]. Studying the distribution of DEGs in the annotation function will clarify the expression of sample differences in gene function. Identified KEGG pathways and GO terms with a corrected p < 0.05 were regarded as significantly enriched.

Quantitative Real-Time PCR (qPCR)
To validate the results of RNA-seq, qPCR was carried out using the ABI StepOne Real-Time PCR system (Applied Biosystems, Foster City, CA, USA), and the expression levels of six randomly selected DEGs (OASL, CYP3A5, LDHA, MYD88, TGM2, and CCL4) with annotations from statistical analysis of RNA-seq were measured using the Fast Start Universal SYBR Green Master kit (Roche, Basel, Switzerland). The gene GAPDH was selected as the reference gene. The relative expression level of each target gene was calculated based on the 2 −∆∆Ct relative expression formula [21]. To investigate the biological properties of different virulence NDV strains, 10-day-old SPF ECE were challenged with a single dose of 7 × 10 3 TCID50 of either virulent strain NA-1 or avirulent vaccine strain LaSota, and the chicken embryos' mortality, HA titers of the harvested allantoic fluid and pathological changes in the chicken embryos were evaluated at 12, 24, 36 and 48 hpi, respectively. HA titers of the harvested chicken embryonic allantoic fluid revealed no differences between NA-1 and LaSota groups, although the HA titers at different stages of infection exhibited significant time-dependent increases after 12 hpi which propagated fastest from 12 to 36 hpi ( Figure 1A). However, the 10-day-old ECE infected with the virulent strain NA-1 started to appear visible typical hemorrhagic pathological changes in the partly brain and body surface at 24 hpi, and exhibited ecchymosis hemorrhages mainly in the brain and body surface at 36 hpi, while no typical hemorrhagic pathological changes were observed in the brain and body surface of avirulent vaccine strain LaSota group ( Figure S1). In addition, the mortality rate reached 50% at 36 hpi and 100% at 43 hpi when the ECE was infected with virulent strain NA-1, while no mortalities were observed following infection with avirulent vaccine strain LaSota and mock ( Figure 1B). Collectively, these chicken embryo infection data demonstrated that distinct biological characterizations of different virulence NDV were present in ovo. Therefore, based on these results, the CEVT at 24 and 36 hpi were selected for further RNA-seq.

Quantitative Real-Time PCR (qPCR)
To validate the results of RNA-seq, qPCR was carried out using the ABI StepOne Real-Time PCR system (Applied Biosystems, Foster City, CA, USA), and the expression levels of six randomly selected DEGs (OASL, CYP3A5, LDHA, MYD88, TGM2, and CCL4) with annotations from statistical analysis of RNA-seq were measured using the Fast Start Universal SYBR Green Master kit (Roche, Basel, Switzerland). The gene GAPDH was selected as the reference gene. The relative expression level of each target gene was calculated based on the 2 −ΔΔCt relative expression formula [21].

Distinct Biological Characterization of NDV Highly Virulent Strain NA-1 and Avirulent Vaccine Strain LaSota In Ovo
To investigate the biological properties of different virulence NDV strains, 10-dayold SPF ECE were challenged with a single dose of 7 × 10 3 TCID50 of either virulent strain NA-1 or avirulent vaccine strain LaSota, and the chicken embryos' mortality, HA titers of the harvested allantoic fluid and pathological changes in the chicken embryos were evaluated at 12, 24, 36 and 48 hpi, respectively. HA titers of the harvested chicken embryonic allantoic fluid revealed no differences between NA-1 and LaSota groups, although the HA titers at different stages of infection exhibited significant time-dependent increases after 12 hpi which propagated fastest from 12 to 36 hpi ( Figure 1A). However, the 10-day-old ECE infected with the virulent strain NA-1 started to appear visible typical hemorrhagic pathological changes in the partly brain and body surface at 24 hpi, and exhibited ecchymosis hemorrhages mainly in the brain and body surface at 36 hpi, while no typical hemorrhagic pathological changes were observed in the brain and body surface of avirulent vaccine strain LaSota group ( Figure S1). In addition, the mortality rate reached 50% at 36 hpi and 100% at 43 hpi when the ECE was infected with virulent strain NA-1, while no mortalities were observed following infection with avirulent vaccine strain LaSota and mock ( Figure 1B). Collectively, these chicken embryo infection data demonstrated that distinct biological characterizations of different virulence NDV were present in ovo. Therefore, based on these results, the CEVT at 24 and 36 hpi were selected for further RNA-seq.

Global mRNA Expression Patterns in Chicken Embryonic Visceral Tissues after NA-1 and LaSota Infection
To make a thorough inquiry about the exact mechanistic differences behind the above observations between infection by virulent and avirulent NDV strains, a transcriptome sequencing for the CEVT samples was performed. When comparing NA-1 to the LaSota,

Global mRNA Expression Patterns in Chicken Embryonic Visceral Tissues after NA-1 and LaSota Infection
To make a thorough inquiry about the exact mechanistic differences behind the above observations between infection by virulent and avirulent NDV strains, a transcriptome sequencing for the CEVT samples was performed. When comparing NA-1 to the LaSota, a total of 2209 DEGs were identified, with 586 (including 290 upregulated and 296 downregulated) and 1623 (including 949 upregulated and 674 downregulated) at 24 and 36 hpi, respectively (Figure 2A-C and Table S1). Among these 2, 209 DEGs, 55 genes were commonly DE at the two-time points, although 531 and 1568 DEGs were uniquely DE at 24 and 36 hpi, respectively ( Figure 2D). a total of 2209 DEGs were identified, with 586 (including 290 upregulated and 296 downregulated) and 1623 (including 949 upregulated and 674 downregulated) at 24 and 36 hpi, respectively (Figure 2A-C and Table S1). Among these 2, 209 DEGs, 55 genes were commonly DE at the two-time points, although 531 and 1568 DEGs were uniquely DE at 24 and 36 hpi, respectively ( Figure 2D).  To investigate the potential biological functionalities of the above DEGs, the functional annotation of these DEGs was performed through a GO enrichment analysis ( Figure 3A,B and Table S2). When comparing NA-1 to LaSota, the top ten terms within the GO category of biological process that were associated with metabolic processes, such as organic substance metabolic process, primary metabolic process, organonitrogen compound metabolic process, and regulation of metabolic processes, were significantly enriched at 24 and 36 hpi, although two terms named with metabolic process and cellular metabolic process were significantly enriched only at 36 hpi (p < 0.05). Meanwhile, common GO terms of cellular components, including organelles, membrane-bounded organelles, cytoplasms, intracel-lular membrane-bounded organelles, and cytoplasmic parts, were significantly enriched at both time points. Conversely, unique GO terms at 24 and 36 hpi were involved in intracellular, intracellular parts, and intracellular organelles, as well as extracellular regions, extracellular region parts, and extracellular spaces, respectively. Interestingly, molecular functions related to bindings and activities were enriched in completely different ways, as shown in binding, ion binding, protein binding, cation binding, metal ion binding, enzyme binding, zinc ion binding, transcription regulator activity, and enzyme regulator activity. Furthermore, protein binding, small-molecular binding, anion binding, nucleotide binding, nucleoside phosphate binding, carbohydrate-derivative binding, ribonucleotide binding, catalytic activity, hydrolase activity, and transferase activity were significantly enriched at 24 and 36 hpi. Therefore, based on the results from the top ten GO categories of DEGs at these two time points, different enriched GO terms were presented between infection by virulent and avirulent strains, with some common occurrences within three main GO categories (including cellular components, molecular functions, and biological processes), regardless of virulence. Moreover, the unique significantly enriched GO terms during an NA-1 infection relative to those that were identified during LaSota infection could be attributed to differences in epigenetic and pathogenesis mechanisms between the virulent and avirulent strains.
Next, KEGG pathway enrichment analysis was applied as an additional way to elucidate the functions of these DEGs, regardless of upregulation and downregulation. According to the KEGG mapping of the DEGs, 37 significantly enriched KEGG pathways (p < 0.05) were obtained during the process of NDV infection at 24 and 36 hpi when comparing NA-1 to LaSota; there were ten at 24 hpi and twenty-seven at 36 hpi ( Figure 3C,D and Table S3). Among the 37 KEGG pathways, 10 of the 27 pathways were related to immune response and 13 had metabolism at 36 hpi, but only 1 of the 10 pathways was associated with metabolism and 1 with the immune response at 24 hpi. At the virus infection later stage of 36 hpi, DEGs were enriched in the innate immune response, and inflammatory responses were found to be mainly associated with cytokine-cytokine receptor interaction, the NF-kappa B signaling pathway, Jak-STAT signaling pathway, TNF signaling pathway, Toll-like receptor signaling pathway, HIF-1 signaling pathway, complement and coagulation cascades, RIG-I-like receptor signaling pathway, NOD-like receptor signaling pathway, and Toll and Imd signaling pathway. On the other hand, the enriched pathways were associated with metabolism, including arginine and proline metabolism, metabolism of xenobiotics by cytochrome, ascorbate and aldarate metabolism, drug metabolism-cytochrome P450, D-Glutamine and D-glutamate metabolism, fructose and mannose metabolism, nitrotoluene degradation and five other related pathways. These results, in conjunction with the severe pathological damage observed at a later stage post-infection in ovo, demonstrate that significantly increased enriched pathways related to metabolism and host innate immune responses could contribute to the pathogenesis of the virulent strain in comparison with the avirulent strain.
To further identify the data from RNA-Seq, a subset of six randomly selected genes (OASL, CYP3A5, LDHA, MYD88, TGM2, and CCL4) with annotations from the statistical analysis of RNA-seq was evaluated by qPCR analysis ( Figure 2E and Table S4). As shown in Figure 2E, six genes exhibited a concordant direction both in RNA-seq and qPCR analysis. The comparable results indicated that the data obtained by RNA-seq were valid.  Next, KEGG pathway enrichment analysis was applied as an additional way to elucidate the functions of these DEGs, regardless of upregulation and downregulation.

Highly Virulent Virus Induces a Robust Innate Immune Response in Chicken Embryonic Visceral Tissues Compared to the Avirulent Vaccine Virus
Previous results demonstrated that infections of virulent NDV elicit strong host innate immune responses in vivo, in vitro, and in ovo, evidenced by the increased expression of inflammatory cytokines, type I and II interferons, IFN effectors, and other host innate immune-related genes and proteins [13,[22][23][24]. In Table 1, high expression levels of 25 kinds  of inflammatory cytokines and their receptors were induced by the highly virulent strain  NA-1 at 36 hpi, including IFNAR1, FLT1, IL10RB, VEGFA, IL8L2, IL13RA1, CSF1, CX3CL1,  IL8L1, CRLF2, TNFSF15, IL10RA, TNFRSF9, C-C motif chemokine 26-like, TNFRSF23,  TGFBR2, IL17RA, IL1B, TNFRSF10B, CCL4, CSF3, IL15RA, EDA2R, IL6, and CCL19. These findings are consistent with a strong innate immune response to highly virulent NDV after respiratory exposure [22,25]. Among these host factors, the capability of IL-6 to stimulate cytokine production, the migration of macrophages, lymphocytes, dendritic cells, and neutrophils into inflammatory sites [26], as well as the action of IL8L1 and IL8L2 (equivalent to IL8 in humans) to stimulate keratinocytes, endothelial and epithelial cells with the resultant movement of lymphocytes into tissues, indicate their potential roles in viral pathogenesis [27]. Under normal circumstances, these cytokines help coordinate the response of the host immune system against infectious agents [28], such as viruses and bacteria. However, sometimes, the host produces excessive amounts of pro-inflammatory cytokines that trigger inflammation and not enough feedback from the anti-inflammatory cytokines that modulate inflammation, in response to infections of highly virulent influenza virus subtypes, H1N1 and H5N1 [29,30]. Therefore, elevations of several inflammatory cytokines seem to be involved in the development of an acute cytokine storm syndrome, which could be the partial cause of death in chicken embryos infected with highly virulent NDV.
Of the upregulated or downregulated DEGs and their enriched KEGG pathways obtained, pathways specifically involved in innate antiviral response, host inflammatory response, and well-defined cytokines were selected for further analysis in this work. According to the KEGG mapping for the upregulated DEGs, a total number of 27 pathways were significantly enriched (p < 0.05) with five at 24 hpi and twenty-two at 36 hpi when comparing NA-1 to LaSota ( Figure 4A,B), while 43 significantly enriched pathways for the downregulated DEGs were obtained with 13 at 24 hpi and 30 at 36 hpi ( Figure 4C,D). Interestingly, when comparing NA-1 to LaSota, 81.82% (18/22) KEGG pathways (enriched by the upregulated DEGs) associated with host innate immune response were induced at 36 hpi, including cytokine-cytokine receptor interaction, the MAPK signaling pathway, NF-kappa B signaling pathway, TNF signaling pathway, and Jak-STAT signaling pathway, whereas no innate immune response-related pathway was enriched at 24 hpi ( Figure 4A,B and Table 1). However, 15.15% (5/33) of KEGG pathways related to the host's innate immune response were enriched based on downregulated DEGs in NA-1 relative to LaSota ( Figure 4C,D and Table 1). In addition, numerous upregulated DEGs were specifically involved in the host immune system, including antigen processing and presentation, complement and coagulation cascades, intestinal immune network for IgA production, and the cytosolic DNA-sensing pathway.
Three classes of pattern-recognition receptors (PRRs), designated NOD-like receptors (NLRs), retinoic acid-inducible gene I (RIG-I)-like receptors (RLRs), and Toll-like receptors (TLRs), were demonstrated to be involved in the recognition of pathogen-associated molecular patterns (PAMPs) in various innate immune cells [31,32]. Among these receptor types, RLRs and TLRs are essential for the production of proinflammatory cytokines and type I interferons (IFNs), whereas NLRs are important to regulate interleukin-1β (IL-1β) production [31][32][33]. In the host, the antiviral innate immune response should be of a suitable magnitude and duration to competently eradicate the invading viruses; however, an overactive innate immune response can cause immune pathology and subsequent undesirable damage [34]. Similarly, in our work, when comparing NA-1 to LaSota at 36 hpi, a total of 23 KEGG pathways related to host innate immune response-including three classes of PRRs signaling pathways, the NF-kappa B signaling pathway, and others-were significantly enriched by DEGs (Figure 4). Therefore, we speculate that the infection of highly virulent NDV causes the hyperactivation of various host innate immunity and promotes the pathogenesis of ND illness.  , TRIM25, IL8L2, IL8L1, MAP3K14, TNFAIP3, PTGS2, VCAM1, NFKB2, IL1B, CCL4, E3  ubiquitin/ISG15 ligase TRIM25-like, PCASP3, TRAF6, BIRC3, MYD88,

Severe Metabolic Disorders in Chicken Embryonic Visceral Tissues Caused by Highly Virulent Virus Compared to the Avirulent Vaccine Virus
Metabolism consists of numerous interconnected cellular pathways within the cells of living organisms to ultimately provide cells with the energy required to sustain life [35]. Over the last decades, it has been confirmed that viruses dramatically promote cellular anabolism for the generation of the macromolecules needed for virion assembly, viral genome replication, and envelopment of the enveloped viral particles, since they do not innately have their own metabolism [36,37]. Therefore, the recognition of how viruses reprogram cellular metabolism, and where in the virus life cycle these metabolic alterations are required, will present a comprehensive interpretation of virus propagation needs and potentially deliver defense strategies against viral invasions. Several host central cellular metabolic pathways, including lipid metabolism, vitamin metabolism, thymidine metabolism, and amino acid metabolism, are significantly reprogrammed by NDV infection [13,38,39]. However, there is very limited knowledge regarding the recognition of how

Severe Metabolic Disorders in Chicken Embryonic Visceral Tissues Caused by Highly Virulent Virus Compared to the Avirulent Vaccine Virus
Metabolism consists of numerous interconnected cellular pathways within the cells of living organisms to ultimately provide cells with the energy required to sustain life [35]. Over the last decades, it has been confirmed that viruses dramatically promote cellular anabolism for the generation of the macromolecules needed for virion assembly, viral genome replication, and envelopment of the enveloped viral particles, since they do not innately have their own metabolism [36,37]. Therefore, the recognition of how viruses reprogram cellular metabolism, and where in the virus life cycle these metabolic alterations are required, will present a comprehensive interpretation of virus propagation needs and potentially deliver defense strategies against viral invasions. Several host central cellular metabolic pathways, including lipid metabolism, vitamin metabolism, thymidine metabolism, and amino acid metabolism, are significantly reprogrammed by NDV infection [13,38,39]. However, there is very limited knowledge regarding the recognition of how different virulent NDV reprogram cellular metabolism and whether highly virulent NDV-induced metabolic pathways alterations could be the main cause of severe systemic pathological changes and high-mortality ND illness. According to the KEGG mapping, in total, 36 metabolic pathways were significantly enriched (p < 0.05) with 7 at 24 hpi and  Figure 4C,D). With the progress of a highly virulent NDV propagation, the numbers of the reprogrammed cellular metabolic pathways increased; therefore, the results presented here are in line with previous studies [38,40].
Many viruses have been indicated to alter various aspects of host central carbon metabolism, including elevated glycolysis and increased pentose phosphate activity, in order to promote the production of amino acids and nucleotides, as well as lipid synthesis [37,41]. However, in our work, 94.44% (34/36) metabolic pathways were enriched by downregulated DEGs at two time points, whereas only two metabolic pathways (including arginine and proline metabolism, as well as amino sugar and nucleotide sugar metabolism) were enriched by upregulated DEGs at 36 hpi rather than 24 hpi (Table 1). To our surprise, among the 34 metabolic pathways associated with downregulated DEGs, 25 pathways were found at 36 hpi ( Figure 4D), a time point with a mortality rate of 50% when the chicken embryo was infected with the highly virulent strain NA-1 ( Figure 1B). Together with the different pathological changes observed at 36 hpi in ovo, we speculated that the infection of highly virulent NDV induces a severe imbalance of reprogrammed cellular metabolic pathways between host cell supplies and virus requirements and aggravated the life-threatening symptoms in the infected cell, while avirulent LaSota results in reprogrammed cellular metabolic pathways with a suitable magnitude and duration.

Conclusions
In conclusion, our study clearly confirms that highly virulent NDV drives robust innate immune responses and severe metabolic disorders to promote the immunopathological damages of ND illness, and thus life-threatening symptoms in severely affected animals. However, more studies concerning the host's innate immune responses and metabolic disorders are needed to further explore the exact mechanisms involved in the interaction between the host and highly virulent NDV as compared to the avirulent NDV.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/v14050911/s1, Figure S1: Pathological changes of the chicken embryos infected with either NDV or mock. The typical hemorrhagic pathological changes were presented in NA-1-infected SPF chicken embryos after 24 hpi, while no pathological changes were observed in both LaSota-infected and mock-infected groups at all time points; Table S1: Detailed list regarding all DEGs in the groups between NA-1 and LaSota at 24 and 36 hpi; Table S2: GO enrichment analysis for all DEGs in the groups between NA-1 and LaSota at 24 and 36 hpi; Table S3: KEGG pathway enrichment analysis for all DEGs in the groups between NA-1 and LaSota at 24 and 36 hpi; Table S4: Primers used in this work.