Changes in the Transcriptome Profile in Young Chickens after Infection with LaSota Newcastle Disease Virus

1. Introduction

Newcastle disease (ND) is a respiratory or systemic disease caused by avian orthoavulavirus serotype 1, commonly referred to as Newcastle Disease virus (NDV) [1]. ND has caused financial losses in the millions of dollars in different countries [2,3] due to increases in bird mortality in poultry flocks, expenses with the outbreak eradication, and trade restrictions during outbreaks [4]. Current ND vaccines do not prevent NDV infection and shedding and might thus mask outbreaks [5]. Improved vaccines would be helpful to avoid production and economic losses. Understanding the avian immune response after infection or vaccination with NDV live vaccines are necessary for the rational design and testing of vaccines.

Understanding early gene expression changes in chickens after vaccination against ND can identify immune-related genes linked to clinical disease prevention and virus shedding that could be used as immunity biomarkers. Transcriptome analysis allows in-depth analysis of gene expression [6] and permits identifying changes in the expression of immune-related genes after viral exposure with high sensitivity [7] without the need to pre-select certain genes and thus bias the results. Previous studies have identified differentially expressed genes in the Bursa, Harderian glands (Hg), lungs, spleen and thymus at two, three, six, and 10 days post infection (dpi) in three to six-week-old birds [8,9,10,11,12,13,14]. However, there is little information about early innate immune response at 12 and 24hours post-infection. To understand the early innate immune response, it is necessary to identify transcriptomic changes in tissues that are considered NDV entry sites. In that sense, lymphoid-associated tissues, and upper respiratory tract such as Hg and tracheas (Tc) have a key role.

Knowing the early innate immune response in day-old-chicks after ND vaccination is also important because, in the field, birds are frequently vaccinated right after hatching. Considering that the immune system of chickens is only fully mature after four to seven weeks of age [15], results obtained in older birds might not be applicable for hatchling chicks. Furthermore, due to the widespread use of vaccines, chicks are likely to have maternal antibodies against NDV, and there is no consensus if its presence can alter the immune response to ND vaccine [16,17].

In the current study, two experiments investigated the early innate immune response by transcriptome analysis in chicks upon vaccination with the NDV LaSota strain. In experiment 1, 13-days-old specific pathogen free (SPF) chickens were vaccinated. Changes in gene transcription 12, 24 and 48 hours post vaccination (hpv) were investigated by transcriptome analysis using RNA sequencing (RNA-seq) and by RT-qPCR. The innate immune genes identified by both methods were compared. Experiment 2 studied the transcriptome 24- and 48-hpv of 1 day-old and 14 days-old chickens with or without maternal antibodies against NDV. The main objective was to identify target innate immune genes expressed shortly after vaccination and recognize potential candidates to be further tested as immunity biomarkers for ND vaccines.

2. Materials and Methods

Birds were kept in Horsfall-type isolators in biosafety level 2 facilities. All birds in both experiments were fed a standard commercial chicken feed that met the National Research Council’s recommended minimum nutrient requirements for laying hens [18]. Feed and water were provided ad libitum. Housing temperature and humidity were set following the breed rearing manual. Birds were randomly selected from each isolator and euthanized with carbon dioxide for sample collection.

Experimental design of experiment 1. A total of 42 white leghorn specific-pathogen-free (SPF) fertile eggs (AVS Bio, Catskill, NY) were hatched and randomly distributed in two groups of 21 birds each. At 13 days- of age, each bird of one group received 10⁷ EID₅₀ of a commercial NDV LaSota vaccine in 100 µl via ocular route (50 µl in each eye). The other group served as an unvaccinated control and was mock-vaccinated with PBS. Twelve, 24- and 48- hpv, five birds per group were euthanized for sampling collection. Both Hg and the whole Tc from the same bird were collected, wrapped in individual aluminum foils, and flash frozen on dry ice immediately after collection. Samples were then stored at -80°C until further processing.

RNA extraction and RNA-seq. Total RNA was extracted from Hg using the RNeasy kit (Qiagen, Hilden, Germany) in a QIAcube automated extraction system (Qiagen) accordingly to the manufacturer`s protocol. For the Tc, total RNA was extracted using the Quick RNA Miniprep kit (Zymo Research, Irvine, CA) following the manufacturer’s instructions.

The RNA concentration of each sample was measured by spectrophotometry (Nanodrop 2000, Thermo Scientific, USA). Three replicates with the highest RNA concentration and A₂₆₀:A₂₈₀ and A₂₆₀:A₂₃₀ ratios in the 2.0 to 2.2 range were selected. mRNA library preparation was done by poly A enrichment and sequencing on a NovaSeq PE150 system (Novogene Corporation Inc., Sacramento, CA). The sequencing data was uploaded to the Sequence Read Archive database (SRA accession number SAMN40643659).

Differential expression of genes, gene ontology, Kyoto Encyclopedia of Genes pathway enrichment, and protein-protein interaction analysis. The index for the chicken reference genome (NCBI GenBank accession number GCF_016699485.2) was generated using HISAT2 v2.0.5 [19] with default parameters. Raw reads were quality trimmed using Trimmomatic v0.39 [20], with standard settings prior to downstream analysis. Trimmed reads were aligned to the reference genome using HISAT2, and the read counts for gene expression were obtained using HTSeq v1.13 [21]. Counted mapped reads were further analyzed using the edgeR package v4.0 [22] in R v4.3.1 [23]. The reads were filtered to exclude genes with low counts and normalized by the trimmed mean of M-values method (TMM) [24]. For each organ and time point, the expression levels from vaccinated were compared to unvaccinated groups considering three replicates per group by Fishers’ exact test. The differential gene expression analysis results were considered significant with an adjusted p-value < 0.05 on the false discovery rate (FDR) and logarithm of fold change (logFC) > ± 0.5.

The top 50 differentially expressed genes (DEGs) from each comparison were selected based on their p-value. They were used as input for mapping Gene Ontology (GO) pathways related to each gene based on the genome wide annotation from chicken package (org.Gg.eg.db) by Bioconductor v3.18 [25]. For GO and Kyoto Encyclopedia of Genes (KEGG) pathway enrichment analysis limma v4.3 [26] was used to identify up- and down-regulated pathways. The GO terms and KEGG pathways with p-value < 0.05 were considered significant. In case the output result showed more than 50 significantly regulated pathways, only the top 50 based on their p-value were considered.

In addition, the top-ranked DEGs were used as input into the STRING database [27] for protein-protein analysis (PPI). The “multiple proteins” function was used with Gallus gallus as target organism. The program standard settings were used to determine the PPI network, with a minimum required interaction score of 0.400, excluding non-experimentally determined protein interactions. The PPI network for each significant contrast was imported to CytoScape, v3.10.1 [28]. The String Enrichment plugin together with the String: protein query database from STRINGapp v2.0.1 [29] were used to retrieve the functionality enrichment map. Genome was used as background. The connectivity cutoff (Jaccard similarity) was 0.8, and node Cutoff p-value < 0.05. CytoCluster app v2.1.0 [30] with the ClusterONE algorithm [31] was used to identify protein clusters. Minimum size = 5, minimum density = 0.05, edge weights = combined_score, EnrichmenMap: similarity_coefficient, similarity = Jaccard similarity, and overlap threshold = 0.8 were the settings used to identify clusters linked to immune-related pathways.

Differential expression of genes associated with immune response by RT-qPCR. cDNA was synthesized from three samples per group and total RNA was obtained using random hexamers and oligo-DT primers with the RT² First Strand Kit (Qiagen) following the manufacturer’s instructions. Quantity and quality of the resulting cDNA was measured by spectrophotometry. All samples were normalized to a concentration of 0.5 µg of total RNA followed by dilution to a total volume of 111 µl with RNAse-free water and 2x RT² SYBR green master mix and stored at -20°C until further analysis. For the RT-qPCR, the RT² Profiler™ PCR array chicken innate and adaptive immune responses kit (Qiagen) was used following the manufacturer`s instructions. Each well contained specific immune related primers to detect cDNA of mRNA from genes associated with immune responses against viruses. Tested genes and the five housekeeping genes used can be found in supplementary Table S1. Threshold cycles (Ct) calculated by the qPCRsoft software v 4.1 (Analytik-Jena AG, Jena, Germany) for each sample were loaded into a Microsoft^® Excel 2007 file provided by Qiagen [32]. Ct values were normalized against the housekeeping genes, and the relative expression of each gene in vaccinated and unvaccinated birds was calculated by t-test [32]. A genomic contamination control, reverse transcription controls, and positive PCR controls included in each plate for each sample were used to validate the results. Gene expression levels presenting at least a twofold difference with significance p-value < 0.05 were considered relevant.

Experimental design of experiment 2. In the second experiment, fertile SPF (AVS Bio) and white leghorn eggs, obtained from a commercial hatchery, were hatched in the same conditions as in experiment 1. The commercial laying hens’ breeders were vaccinated repeatedly against ND and other diseases using live as well as inactivated vaccines. Maternal antibodies against NDV were detected by ELISA in commercial birds but not SPF birds [33]. After hatch, 56 birds of each type were split in two groups of 28 birds each, to be vaccinated when 1-day or 14 days-old. At least, each group of 28 birds was split into vaccinated or unvaccinated groups (n = 14). The experiment followed a completely randomized design in 2 x 2 x 2 factorial arrangement with bird type (SPF or commercial), age at vaccination (1 day- or 14 days-old) and vaccination status (vaccinated or control) as factors.

Vaccinated birds received the NDV Lasota strain following the same procedure as experiment 1. Tissue sampling followed experiment 1 protocol, besides samples were collected 24 and 48 hpv, and stored in DNA/RNA Shield (Zymo Research) for 24 hrs at 4°C. After that, the protective reagent was discarded, and samples stored at -80°C.

Sample processing and data analysis in experiment 2. Total RNA from seven Hg and Tc of each group were extracted as in experiment 1, and five samples were used for all analyses based on the same quality criteria from the previous experiment. The sequencing data was uploaded to the Sequence Read Archive database (SRA accession number SAMN40632518).

Analyses of DEGs, GO, KEGG pathways and PPI, as well as viral load by RT-qPCR and by transcriptome analysis were done as described for experiment 1 with the following modifications in the data analysis:

The DEGs analysis was performed for Harderian glands 24 hpv (hg24), Harderian glands 48 hpv (hg48), Tc 24 hpv (tc24) and Tc 48 hpv (tc48), comparing all vaccinated groups against all controls by glmQLFTest for multivariate analysis using edgeR package. The four lists of DEGs with adjusted p-value < 0.05 on the false discovery rate (FDR) and logarithm of fold change (logFC) > ± 1.5.were used to generate an upset plot with UpSetR v.1.4 [34]. That technique allowed to visualize the number of common genes between the analyzes performed. The DEGs identified in all four analyses were used to get the PPI and enrichment pathway map analysis, performed as mentioned for experiment 1. The pathways were obtained similarly to experiment 1, and the common pathways were obtained by the same method as the one used for common DEGs. In addition, genes and pathways that were differentially expressed when considering bird types and age at infection as well its interaction and the interactions between the three main effects: vaccination, bird type and age at infection were identified.

3. Results

3.1. Experiment 1

3.1.1. Differential Expression of Genes Associated Immune Response, GO, KEGG Pathways by RNA-Seq

Due to achieved RNA yields and integrity, not all samples of each organ, group and timepoint could be analyzed. The numbers of samples that were included are shown in Table 1. On average, each sample yielded 25 million raw read pairs with 95.5% clean reads/raw reads ratio (effective %). No DEGs were detected in the Hg collected 12 hpv comparing vaccinated group to the control. Because of the lack of DEGs, pathways were not analyzed. At 24 hpv, 22 up- and 112 down-regulated DEGs were identified in the Hg of vaccinated birds compared to the control group (Supplementary Table S2). Interferon-induced (IFI), and deltex (DXT) a gene that encodes for E3 ubiquitin-protein ligase present on ubiquitin-proteasome degradation pathway [35], and tumor necrose factor receptor associated factor (TRAF), were the innate-immune related genes up-regulated. No down-regulated genes were identified as immune related. The pathway enrichment analysis identified more than 50 up- and 17 down-regulated pathways. The top fifteen significant pathways, based on their p-values, from vaccinated birds are shown in Supplementary Table S3. Several JAK-STAT related pathways were up-regulated (Supplementary table S3). The PPI did not identify protein clusters linked to immune response. For Tc at 12 hpv, there were no up- or down-regulated DEGs on the vaccinated birds compared to the unvaccinated group.

3.1.2. Differential Expression of Genes Associated with Immune Response by RT-qPCR

Samples from Hg at 12-, and 24-hpv, and from Tc at 24 hpv were compromised due to low quality samples or not enough tissue available to perform the analysis. None of the tested genes were positive in Hg and Tc 24 hpv by RT-qPCR. None of the genes tested in the RT-qPCR analysis were identified as up- or down-regulated by the transcriptome analysis.

3.2. Experiment 2

Differential Expression of Genes, GO, KEGG Pathways and PPI in Response to NDV Vaccination in Harderian Glands and Trachea

On average, each sample yielded 23 million raw read pairs with average Effective% of 97.2. The number of DEGs for each main factor and interactions, categorized by organ and timepoint, is present in Table 2. The number of DEGs from vaccination was between 294 in the Tc at 24 hpv, 168 up-regulated genes and 126 down-regulated genes, and 1175 in the Hg at 48 hpv, 114 up-regulated genes and 1,061 down-regulated genes. Among the four vaccinated groups, there were 2169 unique DEGs. In general, vaccination resulted in an increased number of DEGs in the Hg compared to Tc, and more DEGs 48 hpv than 24 hpv. A similar trend was observed for age at vaccination, with more DEGs due to age in the Hg than in the Tc, and more at 48 hpv compared to 24 hpv. The number of DEGs due to bird type was also higher in the Hg. However, an inverse effect was observed regarding hpv, with a higher number of DEGs due to bird type in the Hg at 24 hpv than at 48 hpv. In the Tc, the previously observed trend of having more DEGs at 48 hpv was seen. Furthermore, there were several interactions between two main factors (vaccination and bird type/vaccination and age) and all three main effects together (vaccination with bird type, and age). For instance, some genes were up-regulated in vaccinated SPF birds, but their expression was down-regulated in vaccinated commercial birds or even genes were up-regulated in 14 days-old vaccinated SPF birds, but down-regulated in 1 day-old vaccinated commercial birds.

Considering all four analyses, Hg at 24 and 48 hpv and Tc at 24 and 48 hpv, twenty-one DEGs were present in both organs and both sampling time points (Figure 1). Twenty of these were consistently either up- or down-regulated in all four analyses, while one, myosin, light chain 4 (MYL4), was up-regulated at Hg 24 hpv but down-regulated in the other three analyses and because of that, was used for PPI analysis but not considered a core gene. The group formed by Interferon alpha inducible protein 6 (IFI6), 2′ -5′ oligoadenylate synthetase like (OASL), poly(ADP-ribose) polymerase family member 14 (PARP14), DExD/H-box helicase 60 (DDX60), DExH-box helicase 58 (DHX58), and interferon regulatory factor 7 (IRF7), among fourteen others common DEGs will be further on referred as the core genes (Supplementary Table S4). The PPI interaction and enrichment map from core genes formed three clusters. They were related to host innate immune response processes either by interferon, cytokine, and caspase regulation or regulation of viral genome replication (Figure 2, Supplementary Table S5).

The number of pathways for each main factor and their interactions, categorized by organ and timepoint, can be found in Table 3. There were 130 unique pathways regulated only by vaccination accounting all organs and hpv. These pathways regulation did not follow the pattern observed for DEGs. Interestingly, an equal number of pathways were regulated due to age at vaccination in the Hg compared to Tc, and 24 hpv compared to 48 hpv. However, more pathways were regulated due to age at vaccination in the Hg compared to Tc, and at 48 hpv than 24 hpv. For bird type, the number of regulated pathways was higher in the Hg than Tc, and at 24 hpv than 48 hpv. Additionally, several pathways showed interactions between two (vaccination and bird type / vaccination and age at vaccination) or all three main effects (vaccination with bird type and age at vaccination). For example, pathways that were up-regulated in 14 days-old vaccinated SPF birds are down-regulated in 1 day-old vaccinated commercial birds. No pathways regulated due to vaccination were found to be common across all organs and hpv’s within all vaccinated groups (Figure 3).

4. Discussion

RNA-seq is renowned for its in-depth transcriptome analysis, whether from tissue samples or single cells. It offers a more comprehensive view of gene expression changes than PCRbased assays [36]. However, this technique is still being refined, and factors like sample size, mRNA extraction and sequencing quality, and complex statistical analysis can impact results. Conversely, RT-qPCR is often regarded as the gold standard for tracking gene transcription. However, RT-qPCR relies on well-defined target genes, so without a prior knowledge, RT-qPCR may fail to detect differentially expressed genes, and the selection of target genes might bias the results.

Experiment 1 aimed to compare RNA-seq and RT-qPCR outputs in chicks to assess the correlation between the two methods. Due to the above-mentioned technical issues, the focus will be on the hg24 hpv RNA-seq results, as it was the group with the higher number of samples.

The findings suggest that chicks vaccinated when 13-day-old do not show significant response to NDV up to 12 hpv, but the response begins to rise between 12 and 24 hpv. One of the identified up-regulated DEGs there was IFI, part of host innate immune system, responsible for reducing viral replication by blocking single-stranded RNA virus fusion protein receptors in chicken cells [37]. This gene’s ability to inhibit viral propagation makes it a potential biomarker for immunity, particularly in relation to virus shedding. Additionally, TRAF gene was overexpressed, indicating enhanced type I interferon production [38]. Previous research identified TRAF downregulation caused by NDV M protein plasmids, as a strategy to mocking the host inflammatory response in vitro [39]. Differences between the current and cited research could have happened due to using different experimental models. TRAF expression profile identified in here indicate a well-functioning host early innate immune response when following NDV vaccination.

Another replication strategy of viruses involves inhibiting cell apoptosis by promoting the degradation of proteins associated with cell death upregulation [40]. NDV specifically manipulates the ubiquitin-proteasome pathway to achieve this [41]. In the current study, there was an DXT upregulation, suggesting an increase in the degradation of proteins related to cell death by the ubiquitin-proteasome route. Even though PPI analysis did not find an increase of protein functional pathways related to these genes, gene ontology results showed upregulation of cell-death pathways. With that said, more studies are needed to elucidate DXT role in host cell death during NDV infection. All significant DEGs can be found at supplementary table S2.

Moving on to experiment 2, the main goal was to examine the early regulation profile of DEGs and pathways following NDV vaccination. This was done to identify core genes and evaluate if they have potential as NDV vaccine biomarkers candidates. Strong biomarkers are the ones that can be used to track protective parameters, such as innate immune response, under various conditions. In that sense, experiment 2 investigated four combinations of age at vaccination and presence of maternal antibodies as factors. The results confirmed changes in the gene expression profile between the groups, including interactions between the three main effects. Many factors like age of the birds at vaccination and maternal immunity were seen interacting. However, 20 genes that were consistently up- or down-regulated, not only in all groups, but also in the Hg and Tc at 24- and 48 hpv were identified.

The output indicated the core genes to be involved with innate immune response regulation by changes in the interferon-I (INF-I) cascade in two ways. Firstly, by modulating genes from the pattern recognition-receptors (PRR) system, that were seen up-regulated. Among them, e Retinoic acid-inducible gene-I (RIG-I)-like receptor (RLR) and Ubiquitin ligases protein families are key components of the PRR system, serving as primary tools for eukaryotic cells to recognize cytoplasmatic viruses [42]. Once a viral particle is recognized, PRR triggers several intracellular routes to activate INF-I production by regulating gene expression of interferon inducible genes in attempt to reduce viral propagation [35]. In experiment 2, PRR-related genes PARP14, OASL, DHX58, and DDX60, were upregulated together with IRF7, a non PRR-related interferon regulator gene.

As previously mentioned, increased INF-I translation is crucial for activating the innate immune responses regulated by INF-I, including the gene expression of the IFI gene family [43]. Regarding this family, IFI6 was seen overexpressed in the present study. Previous research reported IFI6, also known as ISG12, to be upregulated during NDV infection. This gene enhances mitochondrial membrane permeability to induce cell-death in an attempt to shut down several single-stranded RNA viruses infection [44]. Given the importance of apoptosis in decreasing viral shedding, IFI6/ISG12 shows potential as a immunity marker for NDV vaccines but that warrants further investigation. Other identified interferon-inducible genes, IFITM1 and IFIT5, are linked to NDV propagation and have potential as vaccine biomarkers. The former hinders the permeability of the host cell membrane to NDV [45], while the later directly binds to viral nucleotides sequences. This intracellular antiviral capability of IFIT5 leads to the translation of non-functional viral proteins impacting virus replication [46].

It is worth noting that the functional profile of the core genes aligns with three clusters identified by PPI analysis. Furthermore, the enrichment map (supplementary table S5) has identified proteins that, to date, have not been specifically linked to NDV vaccinal response. Among the defense response to virus cluster genes, which include DDX60, DHX58, IFIT5, OASL, CMPK2, and RSAD2, the latter two, have little characterized antiviral action upon NDV infection [11,47,48]. Moreover, USP41 which has been identified as part of interferon regulatory factor, and negative regulation of viral genome replication clusters has not yet been linked to the immune response against NDV in chickens, according to current knowledge.

In conclusion, using RNA-seq, this study identified potential marker candidates for NDV vaccine immunity. The PRR involving RIG-I-like receptor, their interaction with ubiquitination-linked processes, and the subsequent regulation of the immune cascade modulated by INF-I appear to be key biological processes in NDV vaccination, but they require further investigation. The core genes highlighted in the present research could be useful in developing markers with focus on viral replication. However, the relationship between core genes expression profile and viral shedding pattern after vaccination needs to be addressed in further experiments.