Newcastle disease is a highly infectious disease that has a devastating financial impact on the global chicken industry. It is caused by the Newcastle Disease Virus (NDV) or Avian Ortho avulavirus-1, which belongs to Kingdom Orthornavirae, Family Paramyxoviridae and Genus Avian Orthoavulavirus (El-Hamid et al., 2020). While vaccinations are widely available, Newcastle Disease is known to cause fatal effects in poultry. Although the immune response due to vaccination is vital, it is essential to study the differences in the resistance to Newcastle Disease due to the genetic variations among poultry breeds (Kapczynski et al., 2013).

In the previous study by M.S. Tarabany (Mahmoud S, 2019), the authors evaluated the immune response against NDV live attenuated and inactivated vaccines in purebred Leghorn, Fayoumi, male Fayoumi × female Leghorn (FL crossbred) and male Leghorn × female Fayoumi (LF reciprocal crossbred) Chicken. This study concludes that the purebred Fayoumi showed the highest antibody titer and lowest mortality rate, followed by the crossbreeds. In contrast, purebred Leghorn showed the lowest antibody titer and highest mortality rate. This differential immune response in the different breeds of Gallus (Chicken) laid the foundation for further research where the genetic differences were studied.

The harderian gland is the less studied tissue of the chicken, present behind the eye. Although its functions have not been studied well, it is known to play an essential role in immune response. The comparative analysis of the harderian gland to other crucial immune system organs showed that the harderian gland plays a key role in immune response (Vanamamalai et al., 2021). A study on the harderian gland of Leghorn and Fayoumi breed Chicken during Newcastle Disease showed the role of genes differentially expressing between the two breeds. This study shows that the Fayoumi breed, being comparatively more resistant to NDV than Leghorn, expressed more immune system-related genes than Leghorn during normal and Newcastle Disease Virus challenge conditions (Deist et al., 2018).

Long noncoding RNAs (lncRNAs) are transcripts with a length greater than 200 nucleotides and do not code for a functional protein. They are synthesised by RNA Polymerase - II and possess a 5′cap. Several lncRNAs, known as poly-adenylated lncRNAs, possess a 3′poly-A tail, while some non-polyadenylated lncRNAs are known to contain a triple helix (Wilusz et al., 2012). Although there is a debate regarding the functionality of lncRNAs, the growing number of lncRNA research shows the importance of lncRNAs. They can fold into complex structures, which enables them to interact with different molecules like DNA, other types of RNAs and even specific proteins (Anna et al., 2018). The primary function of lncRNAs is gene regulation at multiple levels, including DNA, RNA, and protein. The lncRNAs are involved in various biological activities, including gene activation and inactivation by modulating the chromatin, interacting with DNA and histone proteins, modulation of the transcriptional activity by interacting with mRNAs and transcriptional machinery including transcription factors, post-transcriptional regulation by sequestering the proteins to form lncRNA-protein complexes which modulate splicing and turnover of mRNAs and translational regulation by modulating translation assembly and binding directly with other types of RNAs. The lncRNAs are also involved in regulating miRNAs by blocking the miRNAs through miRNA sponging. The lncRNAs can also regulate proteins by binding through RNA-protein interactions. This gene regulation of lncRNAs, in turn, regulates various physiological activities of cells, including cell differentiation, cell growth and responses to various stimuli. The lncRNAs also play vital roles in muscular, nervous, cardiovascular and immune systems (Statello et al., 2021). In immune system-related processes, lncRNAs play significant roles in haematopoietic development, myeloid differentiation, CD4+ t-cell differentiation, inflammation activation and restriction through lncRNA-DNA, lncRNA-RNA and lncRNA-protein interactions (Chen et al., 2017).

In our recent study on trachea transcriptome (Vanamamalai et al., 2023), we have identified lncRNAs and performed differential expression analysis in three different types–challenge-based analysis, i.e., between non-challenged and challenged samples, breed-based analysis, i.e., between Leghorn and Fayoumi, and timepoint based analysis, i.e., between different timepoints–2, 6, 10 days post-challenge (DPC). However, in our previous study on the transcriptomic data of the harderian gland (Vanamamalai et al., 2021), we have identified the genes and long noncoding RNAs differentially expressing between normal samples and NDV-challenged samples of Leghorn and Fayoumi breeds separately, i.e., challenge-based analysis only, which made us understand the role of genes and lncRNAs during the infection stages in both breeds. In continuation to that, in this work, we analysed the same data to identify the genes and lncRNAs that were differentially expressed between Leghorn and Fayoumi samples at both normal and NDV-challenged conditions at three different timepoints, i.e., Breed based analysis and between different timepoints, when challenged with NDV, i.e., timepoint based analysis. We predict that this helps us to understand the role of genes and long noncoding RNAs, which could be involved in the Fayoumi breed becoming comparatively more resistant to NDV than the Leghorn breed chicken.

The transcriptomic data of the harderian gland of Leghorn and Fayoumi chicken was downloaded from the publicly available database EBI-ENA, with project ID PRJEB22672 (Deist et al., 2018). The dataset comprises 94 samples: 46 normal samples and 48 NDV-challenged samples. The details of the sample dataset are shown in Supplementary Table S1. The computational pipeline FHS pipe, as explained in our previous work (Vanamamalai et al., 2023).

The pipeline includes various steps – Quality filtering using Fastp (Shifu et al., 2018), mapping using HISAT2 v2.2.1 tool (Kim et al., 2019) and the latest reference genome of Chicken (GRCg7b), assembly using Stringtie v2.1.4 tool (Pertea et al., 2015). The assembled transcripts were then annotated with different class codes using the tool GFFCompare v0.11.2 (Pertea and Pertea, 2020). The sequences of class codes “I,” “U,” and “X” were selected and extracted using the Bedtools v2.28.0 (Quinlan and Hall, 2010). The extracted sequences were subjected to various filters–length filter (minimum 200 nucleotides), ORF filter (maximum 300 nucleotides/100 amino acids), Pfam filter (no hits) using RPSBlast v2.11.0+ (Christiam et al., 2009) against the Pfam database (Mistry et al., 2021) with an e-value 1e-3 and coding potential filter using Coding potential calculator 2 (CPC2) (Kang et al., 2017). The sequences with noncoding tags were filtered and subjected to BlastN v2.11.0+ (Christiam et al., 2009) with 100% query coverage and 100% identity against other noncoding RNA databases like eukaryotic transfer RNA database (Chan and Lowe, 2009), Silva ribosomal RNA database (Christian et al., 2013) and miRbase (Kozomara et al., 2019) as final potential long noncoding RNAs. These are searched against NONCODE database version 6 (Zhao et al., 2021) using BlastN (Christiam et al., 2009) with 100% query coverage and 100% identity to identify the known and novel lncRNAs. The identified lncRNAs were plotted across the Chromosomes according to their locus positions using Phenogram (Wolfe et al., 2013). Previously (Vanamamalai et al., 2021), the lncRNAs were extracted separately at each timepoint, which generated many duplicates between each timepoint. So, now lncRNAs are identified from the entire dataset at once, which gives a non-redundant list of lncRNAs.

FHSpipe generates read counts for the genes and the long noncoding RNAs, which were further used to perform differential expression analysis using edgeR v3.34.0 with Generalised linear model (GLM) (Robinson et al., 2010). Contrasts were written to identify the genes differentially expressed between Leghorn and Fayoumi at both conditions (Breed-based analysis) at three timepoints and between different timepoints of challenged samples (Timepoint-based analysis). The genes and lncRNAs with an FDR value below 0.05 were filtered as significant differentially expressed genes (DEGs) and significant differentially expressed lncRNAs (DElncRNAs). Using the tool Circos v0.69.8 (Krzywinski et al., 2009), the chromosomal localisation of differentially expressed genes and differentially expressed lncRNAs was plotted.

As described earlier (Vanamamalai et al., 2023), GO functional annotation of the identified DEGs was performed using a standalone tool, OmicsBox v2.1.14 (BioBam Bioinformatics, 2019), which includes several steps: BLAST (Altschul et al., 1990), Interpro scan (Philip et al., 2014), Gene ontology mapping and annotation (Gotz et al., 2008), eggNOG mapping (Huerta-Cepas et al., 2019), and combined pathway analysis against Reactome (Fabregat et al., 2018) and KEGG (Kanehisa and Goto, 2000) databases.

GO functional enrichment analysis was performed separately of all the conditions using the GSEA (Gene Set Enrichment Analysis) module in OmicsBox (Subramanian et al., 2005) to detect the significantly different Gene Ontologies and pathways in challenge-based, breed-based and timepoint-based analysis. Various options, including the number of gene set permutations, weighted enrichment statistic and cutoffs (p < 0.05, FDR<0.25, |NES|>1), were used as mentioned previously (Vanamamalai et al., 2023).

WGCNA (Weighted Gene Correlation Network Analysis) v1.70.3 (Langfelder and Horvath, 2008) was used to perform the co-expression analysis of DElncRNAs with the DEGs, which is further used for functional annotation of the DElncRNAs. Co-expression analysis was conducted individually for each condition for all the DEGs and DElncRNAs. The interactions between DEGs and DElncRNAs, involved in the Immune system processes and with a minimum weight cutoff of 0.01, were chosen and visualised using Cytoscape v3.9.1 (Shannon et al., 2003) in the form of networks.

Based on chromosomal localisation, the DEGs and DElncRNAs co-expression interactions can be cis and trans interactions. The co-expressing DEG and DElncRNA were found on the same chromosome in cis, whereas they were found on separate chromosomes in trans interactions (Guttman and Rinn, 2012; Zhao et al., 2020). The type of interaction between the DEGs and the co-expressing DElncRNAs was determined through their chromosomal positions. The functions of DElncRNAs were predicted using the functions of the DEGs co-expressing with these DElncRNAs. The pathways for each gene obtained from Gene-ontology analysis were taken and assigned to the respective co-expressing lncRNA.

The gene-transcription factor interaction was determined using the standalone tool MEME-suite v5.3.2 (Bailey et al., 2009). As mentioned in our previous studies (Vanamamalai et al., 2021; Vanamamalai et al., 2023), the 5′UTR of 5 KB size of all the DEGs were extracted and analysed for the motifs using the MEME tool under MEME-suite. Using the Tomtom tool of MEME-suite, the obtained motifs were then compared to the JASPER2022 Vertebrates database to identify the potential transcription factors binding to these motifs (genes). These were scanned against the Animal transcription factor database v4.0 (Shen et al., 2023) to obtain the transcription factors of chicken. This data was used to plot a gene-TF network using Cytoscape v3.9.1 (Shannon et al., 2003).

The microRNAs interacting with DEGs were obtained using miRNet v2.0 (Chang et al., 2020), an online platform for miRNA target prediction. This tool uses the miRNA target gene data collected from well-annotated databases - miRTarBase v8.0, TarBase v8.0 and miRecords. This tool takes a list of either official gene symbols Entrez IDs or Ensembl IDs as input to predict gene-miRNA interaction. The output was in the CSV file format, which contains the details of miRNA ID and accession, target gene symbol and ID, Predicted miRanda Score, Literature ID and tissue type. This data is used to plot a gene-miRNA network using Cytoscape v3.9.1 (Shannon et al., 2003).

The coordinates of each QTL (Quantitative Trait Locus) were downloaded from the Animal QTL Database (Zhi-Liang et al., 2019), DEGs and DElncRNAs were collected, and QTLs were assigned to each of the DEGs using an in-home built Python script based on the positions of the QTL and DEGs.

The NCBI primer design tool (Ye et al., 2012) was used to design the primers to be used in the RT-PCR and the information was mentioned in Supplementary Table S2. The experimental design and methodology for the lentogenic NDV challenge were used as previously described in detail in our previous work on trachea transcriptome analysis (Vanamamalai et al., 2023). This work was performed adhering to the ethical guidelines approved by the Institute Animal Ethics Committee (IAEC/DPR/20/2). Briefly, chickens aged 21 days were grouped into two groups of 30 each–one group was inoculated with 200 μL of EID 50 ≥ 106 per dose of NDV of LaSota strain, with 50 μL inoculated into each eye and nostril. Similarly, the other group was inoculated with 200 μL of phosphate-buffered saline (PBS). Hemagglutination inhibition (HI) and indirect ELISA determined the ND antibody levels as described earlier (Vanamamalai et al., 2023). Harderian gland tissue was collected from both the groups on 2, 6 and 10 DPC. Total RNA was isolated from the tissue using the Viral RNA purification kit (Himedia Pvt. LTD.) and reverse transcribed to cDNA using the High-capacity cDNA Reverse transcription Kit (Applied Biosystems, United States). Insta Q96™ real-time PCR machine (Himedia, India) was used to perform real-time PCR using Maxima SYBR Green/ROX qPCR Master Mix (2X) (MBI Fermentas, United States). The 2−ΔΔCt method was used to calculate relative expression between respective timepoints with housekeeping gene β actin as the internal control. GraphPad Prism (GraphPad Software, 2024) was used to visualise the results as bar plots, including the standard error of the mean (SEM) and the significance values obtained using multiple unpaired t-tests.

In quality filtering, an average of about 95.67% of the reads had passed the quality filter in the 94 samples, ranging from a minimum of 73.05% to a maximum of 98.63% in different samples. The average Q30 base content was found to be 95.25% across all the samples, and the average GC content was found to be 47.49%. The mapping percentage of all the samples was in the range of 71.05%–96.28%. The average mapping percentage was 90.58%, with a median quality of 92.47%. The number of assembled transcripts in the first run of the assembly was in the range of 38,547–49,406 transcripts. After reassembly using merged GTF, the number of transcripts was found to be 28,788 in all the samples. The quality and assembly details of each sample are mentioned in Supplementary Table S3.

A total of 16 different class codes were identified, of which 34,309 sequences of class codes I, U and X were extracted initially. Figure 1A shows the pie chart showing the proportions of the 16 class codes. After processing the data with length, ORF, PFAM and coding potential filters, 29,047 sequences were obtained. Out of the 29,047 final lncRNAs, about 13,261 sequences belong to class code I, 13,445 to class code U, and 2,341 to class code X, as shown in Figure 1B. The chromosomal localisation of the identified lncRNAs plotted using the tool Phenogram were shown in Figure 1C. With default blast parameters and no identity and coverage filter, about 11 sequences showed hits against the transfer RNA database of chicken and no hits against ribosomal RNA and miRbase databases. While using 100% identity and coverage filters, no sequences were similar to transfer RNA, ribosomal RNA and miRbase databases. The BlastN against NONCODE database version 6 showed that about 23,240 sequences showed no hits and 5,807 sequences showed hits, of which 160 sequences showed hits with 100% query coverage and 100% identity. These details are shown in Table 1 Detailed information on all the identified long noncoding RNAs is mentioned in Supplementary Table S4.

Figure 1. Figure showing the (A) Pie chart of the classification of all the transcripts into 16 class codes, (B) Pie chart of the classification of the identified lncRNAs into 3 different types—Intronic, Intergenic, and Anti-sense, (C) Chromosomal positions of all the identified long noncoding RNAs.

Table 1. Table showing the number of transcripts discarded in different steps in the lncRNAs identification pipeline.

Several DEGs and DElncRNAs were obtained at each timepoint, including upregulated and downregulated sequences, as shown in Table 2A for breed-based analysis and Table 2B for timepoint-based analysis. Supplementary Figure S1 shows the heatmap of the read counts of the DEGs, and Supplementary Figure S2 shows that of the DElncRNAs. Figure 2 shows the chromosomal localisation of the DEGs (A, C) and DElncRNAs (B, D) obtained in breed-based analysis (A, B) and timepoint-based analysis (C, D). More DEGs and DElncRNAs were identified in breed-based analysis than in timepoint-based analysis. The highest number of DEGs and DElncRNAs were identified in 2DPC in challenged, followed by 2DPC non-challenged. These results show that both breeds respond differently right from the initial stage of infection. In timepoint-based analysis, Leghorn showed the highest number of DEGs and DElncRNAs between 2 DPC and 6 DPC compared to Fayoumi. Although a similar number of DEGs and DElncRNAs were identified between 2 DPC and 10 DPC.

Table 2. Table showing the number of differentially expressed genes and lncRNAs at each condition. A: Breed-based analysis. B: Timepoint-based analysis.

Figure 2. Figure showing the synteny plot of the chromosomal localisation of Differentially expressed genes (A, C) and differentially expressed lncRNAs (B, D) obtained in breed-based analysis (A, B) and timepoint-based analysis (C, D).

The functional annotation of the DEGs showed they were involved in several pathways. Figure 3 shows the Level 2 Gene ontology annotation in (A) Biological process, (B) Molecular function and (C) Cellular component plotted using GraphPad Prism (GraphPad Software, 2024). The Biological process distribution shows that 19% of the sequences were annotated with Cellular process (GO:0009987), followed by 15% with metabolic process (GO:0008152) and 2% (78) of the DEGs with Immune system process (GO:0002376). The molecular function distribution shows that 60% of the sequences were annotated with Binding (GO:0005488) and 40% with Catalytic activity (GO:0003824). The cellular component distribution shows that 74% of the sequences were annotated with Cellular anatomical entity (GO:0110165) and 26% with protein-containing complex (GO:0032991). Figure 4 shows the number of pathways per category annotated to each of the DEGs using Reactome and KEGG databases plotted using GraphPad Prism (GraphPad Software, 2024). In KEGG categories, 91 out of 247 pathways were found to be under the “Metabolism” category, followed by the “Human diseases” category with 63 pathways. The least number of pathways were found in the “Genetic information processing” category, with only 8 pathways. In the Reactome categories, 170 out of 968 pathways were found to be under the “Signal transduction” category, followed by the “Metabolism” category with 138 pathways. The least number of pathways were found in “Digestion and absorption,” with only one pathway. In addition, about 117 pathways were found in the “Immune system” category. These 117 immune system pathways were annotated to 240 DEGs. In addition, 28 DEGs were found to be annotated with both Level 2 biological process GO Immune system process (GO:0002376) and immune system (Reactome) category pathways; out of these, five were novel genes with no known official gene symbol. About 290 DEGs were found to be annotated with either immune-related GO terms or immune system pathways. The details, including Gene Ontologies (BP, MF, CC), EGGNOG annotation, pathway name, pathway category, and expression values (log2 fold change) of the DEGs identified in breed-based analysis and timepoint-based analysis, were mentioned in Supplementary Table S5.

Figure 3. Figure showing the pie charts representing (A) the distribution of Level-2 Biological process GOs, (B) sequence distribution of Molecular function GOs and (C) Cellular component GOs annotated to DEGs using OmicsBox.

Figure 4. Figure showing the bar chart distribution of the categories of KEGG and Reactome pathways annotated to DEGs annotated using OmicsBox.

Several enriched/over-represented Gene ontologies and pathways were identified using functional enrichment analysis after filtering the data with p < 0.05, FDR<0.25 and |NES|>1. Based on the normalised enrichment score (NES), the enriched GOs and pathways with positive NES were governed by upregulated genes, and those with negative NES were governed by downregulated genes. The highest number of enriched GOs were found to be biological process GOs and the highest number of enriched pathways were found to be Reactome pathways. In breed-based analysis, the highest number of enriched GOs were identified in challenge 2-day data, while there were no enriched GOs in challenge 6-day data. In the case of pathways, the highest number of enriched pathways were found in challenge 2-day data and the least in non-challenged 10-day data. In challenge 2-day data, a greater number of the enriched biological process GOs were obtained under Cellular process level 2 GO, followed by the metabolic process. Several developmental, regulatory GOs were also observed, while very few enriched GOs were identified under response to stimulus GO. Most of the enriched pathways were found to be under Metabolism (KEGG), Human diseases (KEGG), Metabolism (Reactome) and Disease (Reactome) categories. A few enriched pathways were found under Environmental information processing (KEGG), Organismal systems (KEGG), Signal transduction (Reactome), Organelle biogenesis and maintenance (Reactome), Transcription (Reactome), Cell cycle (Reactome), Haemostasis (Reactome) and Immune system (Reactome) categories. In challenge 6-day data, although no enriched GOs were mentioned above, several enriched pathways were under the Disease (Reactome) and Signal transduction (Reactome) categories. A few enriched pathways were found under Metabolism (Reactome), Metabolism (KEGG), Environmental Information Processing (KEGG) and Human Diseases (KEGG) categories. In the 10-day challenge data, the highest enriched GOs were identified in Localization Level 2 GO, followed by Metabolic process, and only 1 GO in Cellular process GO. Most of the enriched pathways were found under the Disease (Reactome) category, and only 1 enriched pathway was found under each of the Signal Transduction (Reactome) and Human diseases (KEGG) categories. In non-challenged 2-day data, all the enriched GOs were under Cellular process GO. Most of the enriched pathways were found to be under the Metabolism (KEGG) and Disease (Reactome) categories. One enriched pathway was found under the Signal Transduction (Reactome), Organelle biogenesis and maintenance (Reactome), and Immune System (Reactome) categories. In non-challenged 6-day data, 2 enriched GOs were found under Metabolic process GO and 1 each under Cellular process and Localisation GO. Most of the enriched pathways were found under the Disease (Reactome) category, and few enriched pathways were found under Metabolism (Reactome), Transcription (Reactome), Organelle biogenesis and maintenance (Reactome), Haemostasis (Reactome), Metabolism (KEGG), Human diseases (KEGG) and Organismal systems (KEGG) categories. In non-challenged 10-day data, 2 enriched GOs under Cellular process GO and 1 each under response to stimulus, Metabolic process and Localization GOs were found. Only 1 enriched pathway was found in this data under the Metabolism (KEGG) category. In timepoint-based analysis, in Fayoumi, there were no enriched GOs in 2 vs. 6-day data, and 1 each enriched pathway was found under Human diseases (KEGG), Vesicle-mediated transport (Reactome), Cellular responses to stimuli (Reactome) and transport of small molecules (Reactome). In 2 vs. 10-day data, more enriched GOs were identified under Cellular process GO, followed by developmental process and biological regulation GOs. Several enriched pathways were found under the Disease (Reactome) category, and a few enriched pathways were also found under Metabolism of RNA (KEGG), Signal transduction (Reactome) and Cellular responses to stimuli (Reactome). There were no enriched GOs and pathways in 6 vs. 10-day data. However, in Leghorn, in 2 vs. 6-day data, the highest number of enriched GOs were obtained under cellular process and developmental process GOs. Several metabolic, developmental, and regulatory GOs were also observed; only 1 enriched GO was identified under the response to stimulus GO. The most enriched pathways were found under the Disease (Reactome) category. A few enriched pathways were also found under Metabolism (KEGG), Human Diseases (KEGG), Metabolism (Reactome), Cellular responses to stimuli (Reactome), Transcription (Reactome), Muscle contraction (Reactome), Chromatin organisation (Reactome), Developmental biology (Reactome) and Haemostasis (Reactome). In 2 vs. 10-day data, only 2 enriched GOs were obtained, one each under response to stimulus and metabolic process GOs. The most enriched pathways were found under the Disease (Reactome) category. A few enriched pathways were also found under Metabolism (KEGG), Immune system (Reactome), Vesicle-mediated transport (Reactome), Metabolism of proteins (Reactome) and Haemostasis (Reactome). Finally, in 6 vs. 10-day data, only 2 enriched GOs were obtained, one each under biological regulation and metabolic process GOs. However, there were no enriched pathways. The number of the total and enriched Gene Ontologies (BP, MF, CC) and Pathways (KEGG, Reactome) are shown in Table 3. Detailed information on enriched Gene ontologies was mentioned in Supplementary Table S6, and enriched pathways were mentioned in Supplementary Table S7.

Table 3. Table showing the number of Total and Enriched Gene Ontologies, i.e., Biological Process (BP), Molecular Function (MF) and Cellular Component (CC), Pathways, i.e., KEGG and Reactome identified in all the conditions in (A) Breed based and (B) Timepoint based analysis.

In both breed-based and timepoint-based analysis, several co-expression modules were identified between DEGs and DElncRNAs. Supplementary Table S8 shows the statistics of the co-expression analysis of all the datasets. There were no outliers in all the datasets. Scale-free topology model, fit index cutoff, was set to 0.8, and soft power was chosen according to the index. In a few cases where the index was below the cutoff, a soft power of 9 was chosen, as mentioned in the user manual. The minimum module size was 15 for Breed-based analysis and 9 for timepoint-based analysis. In breed-based analysis, the highest number of modules were identified in challenged 2-day and 6-day samples. In timepoint-based analysis, most modules were identified in Leghorn 2 vs. 6-day samples.

Most of the interactions were found to be trans-regulatory interactions, as shown in Table 4. About 93% were trans interactions, and 7% were found to be cis interactions. The lncRNAs were found to be downstream of the co-expressing gene in about 49.11% of cis interactions, upstream of the co-expressing gene in about 48.16% interactions, within the co-expressing gene in about 2.28% interactions, lncRNAs were spanning the 5′end in about 0.06% interactions and 3′end in about 0.38% interactions. The lncRNAs involved in cis interactions were functionally annotated from the functional annotation of the DEGs co-expressing with this lncRNA at all the different conditions. The results are mentioned in Supplementary Table S9. In breed-based analysis, challenge 2-day data showed the highest number of cis-interactions, DElncRNAs and DEGs, while non-challenged 10-day data showed the least number of cis-interactions, DElncRNAs and DEGs. In timepoint-based analysis, Leghorn 2 vs. 6-day data showed the highest number of cis-interactions, DElncRNAs and DEGs. In contrast, the Fayoumi 6 vs. 10 days data showed the least cis-interactions, DElncRNAs, and DEGs. There were no cis-interactions in Leghorn 6 vs. 10-day data. About 852 different lncRNAs were found to co-expressing with DEGs annotated with either an immune-related GO term or an immune system pathway, of which 435 lncRNAs were involved in more positive coregulations and 417 lncRNAs were involved in more negative coregulations.

Table 4. Table showing the number of cis and trans pairs of co-expressing genes and lncRNAs identified in (A) Breed-based and (B) Timepoint-based analysis.

A total of 10 motifs were identified in the 5’ UTR of the DEGs expressed across all the conditions. Of these, 4 motifs (motifs 1, 2, 4 and 5) showed hits with transcription factors with a q-value of <0.05. About 47 transcription factors were identified for all the DEGs across all the conditions. Motif 1 showed hits with 21 transcription factors, motif 2 with 17 transcription factors, motif 4 with 30 transcription factors and motif 5 with 13 transcription factors. These transcription factors belong to 4 families: bHLH, E2F, ZBTB and zf-C2H2. The details of all the transcription factors identified using TomTom were mentioned in Supplementary Table S10A and details of transcription factors associated with immune-related genes were mentioned in Supplementary Table S10B.

A total of 661 miRNAs were found to interact with 431 DEGs across all six conditions. A total of 688 DEGs showed no miRNAs associated with them. About 79 miRNAs were found to target immune-related genes. The genes KIF13B and PRELP were targeted by the highest number of miRNAs, 50 miRNAs, and the miRNA gga-mir-1776 was found to target the highest number of DEGs, 47 DEGs. The details of microRNAs associated with DEGs were mentioned in Supplementary Table S10C.

The DEGs annotated with Immune-related Gene Ontology terms and Immune system pathways, along with co-expressing DElncRNAs, the associated transcription factors and miRNAs, were selected and plotted as network figures using Cytoscape (Shannon et al., 2003). In the network plot, Immune-annotated DEGs were represented by blue diamond nodes, co-expressing DElncRNAs with black spherical nodes, associated transcription factors with red triangular nodes, microRNAs with pink arrow nodes and biological process GOs with green rectangular nodes. The lncRNAs, TFs and microRNAs surrounding a specific gene were interacting with that particular gene and those in a separate cluster were interacting commonly with multiple genes. In breed-based analysis, six different network plots were generated and represented in Figure 5. In this figure, A represents network data obtained in non-challenged 2-day data. Of the 10 DEGs, the gene LGMN showed highest number of interactions with lncRNAs (557) and gene CD47 showed the least (96). The genes CD47, GPRC5B, LGMN and TF showed several lncRNAs which were co-expressing only with the respective gene. The lncRNAs and TFs in the central region denote the common one between all the DEGs. In Figure 5B represents network obtained from non-challenged 6-day data. Of the 11 DEGs, the gene DDOST showed highest interactions (179) and gene PECAM1 showed the least (11). The genes FLNB, GPRC5B, LGMN, MSTRG.22587, PECAM1 and TF showed several lncRNAs which were co-expressing only with the respective gene. In Figure 5C represents network obtained in non-challenged 10-day data. Only 1 gene GPRC5B was identified. In Figure 5D represents network obtained in challenged 2-day data. Of the 19 DEGs, the gene DDOST showed highest interactions (1,035) and gene FGL2 showed the least (51). The genes GPRC5B, LGMN and TF showed several lncRNAs which were co-expressing only with the respective gene. In Figure 5E represents network obtained in challenged 6-day data. Of the 4 DEGs, the gene GPRC5B showed highest interactions (104) and gene MSTRG. 19088 showed the least (1). All the 4 genes - ANXA2, GPRC5B, MSTRG.19088 and PADI2 showed several lncRNAs which were co-expressing only with the respective gene. In Figure 5F represents network obtained in challenged 10-day data. Of the 2 DEGs, MSTRG.9319 showed highest interactions (100) and gene GPRC5B showed the least (93). Both the genes showed several lncRNAs which were co-expressing only with the respective gene.

Figure 5. Figure showing the Network plot between immune annotated differentially expressed genes (Blue colour Diamond shaped), co-expressing differentially expressed lncRNAs (Black colour spherical shaped), associated transcription factors (Red colour triangular shaped), microRNAs (Pink colour arrow shaped) and biological process GOs (Green colour rectangular shaped) identified in breed based analysis. (A) non-challenged 2-day, (B) non-challenged 6-day, (C) non-challenged 10-day, (D) challenged 2-day, (E) challenged 6-day, (F) challenged 10-day.

In timepoint-based analysis, only four different network plots were generated and represented in Figure 6, as there were no interactions in Leghorn 6 v/s 10 DPC and Fayoumi 6 v/s 10 DPC data. In Figure 6A represents network obtained in Leghorn 2 v/s 6-day data. Of the 18 DEGs, the gene DDOST showed highest interactions (558) and genes CD47 and FGL2 showed the least (33). The genes ANXA2, CD74, DDOST showed several lncRNAs which were co-expressing only with the respective gene. In Figure 6B represents network obtained in Leghorn 2 v/s 10-day data. Of the 10 DEGs, the genes DDOST and TF showed highest interactions (214) and genes CD74 showed the least (11). The genes MSTRG. 22587 and TF showed several lncRNAs which were co-expressing only with the respective gene. In Figure 6C represents network obtained in Fayoumi 2 v/s 6-day. Only 1 gene MSTRG. 20675. In Figure 6D represents network obtained in Fayoumi 2 v/s 10-day. Of the 5 DEGs, the gene TF showed highest interactions (264) and the gene MSTRG. 20675 showed the least (113). The genes DDOST, SH3BP5 and TF showed several lncRNAs which were co-expressing only with the respective gene.

Figure 6. Figure showing the Network plot between immune annotated differentially expressed genes (Blue colour Diamond shaped), co-expressing differentially expressed lncRNAs (Black colour spherical shaped), associated transcription factors (Red colour triangular shaped), microRNAs (Pink colour arrow shaped) and biological process GOs (Green colour rectangular shaped) identified in timepoint based analysis. (A) Leghorn 2 v/s 6-day, (B) Leghorn 2 v/s 10-day, (C) Fayoumi 2 v/s 6-day, (D) Fayoumi 2 v/s 10-day.

About 3,972 QTLs were obtained, including 729 Exterior, 1,006 Health, 654 Physiology and 1,583 Production QTLs. Of all six conditions, challenged 2-day samples showed the highest number of all types of QTLs (2032), and normal 10-day showed the lowest number of all QTLs (81). The details of the number of different types of QTLs obtained in each of the conditions are shown in Table 5.

Table 5. Table showing the number of differentially expressed genes associated with different types of QTLs obtained in (A) Breed-based analysis and (B) Timepoint based analysis.

The relative expression values between respective timepoints were calculated using the 2-ΔΔCt method individually for each condition, i.e., 2 vs. 10-day, 2 vs. 6-day and 6 vs. 10-day Figure 7 shows the relative expression values of DElncRNAs and co-expressing DEGs identified in A: 2 vs. 10-day, B: 2 vs. 6-day and C: 6 vs. 10-day timepoints. The relative expression values determined through RT-PCR were consistent with the in silico analysis. In addition, the unpaired t-test results show that all the 7 lncRNAs and 5 genes showed significant differences.

Figure 7. Figure showing the relative expression (2-DDCT) values of DElncRNAs and co-expressing genes identified in (A) 2 vs. 10-day, (B) 2 vs. 6-day and (C) 6 vs. 10-day conditions, with 2 days data represented as green bars, 6 days data as orange bars and 10 days data as blue bars. SEM values were plotted as error bars, and significance values (p) were represented by–*−p≤0.05, **−p≤0.01 and ***−p≤0.001.

Apart from this, several transcription factors and microRNAs interacting with the DEGs were also identified. Transcription factors interacting with the 28 immune annotated genes were identified, of which 5 DEGs showed no motifs with transcription factors, 2 DEGs showed all 4 motifs and all 47 transcription factors were found to interact with them. The remaining 21 DEGs showed 1-3 motifs and 12–46 transcription factors. Of the identified 47 transcription factors, 17 belong to zf-C2H2 (zinc finger Cys2–His2) family, which are known to play a vital role in disease and developmental process (Han et al., 2016), 14 belong to bHLH (basic helix-loop-helix) family, which play vital role in cell determination, differentiation and proliferation (Ledent and Vervoort, 2001), 1 belongs to ZBTB (Zinc finger and BTB domain containing) family, which play a vital role in differentiation, development and functioning of T-lymphocytes and related molecular mechanisms (Cheng et al., 2021), 1 belongs to E2F family, which regulate cell cycle and apoptotic processes (Xanthoulis and Tiniakos, 2013) and 14 were of unknown families. Along with the positive co-expressing lncRNAs, further studies can also use this transcription factor information to enhance the expression of the immune-related DEGs in the susceptible Leghorn breed to improve the overall resistance. Along with this, several microRNAs were identified to be associated with DEGs. In the case of the 28 immune-related genes, 13 genes showed no miRNAs. The remaining 15 genes were identified to be targeted by 79 different miRNAs. Of them, the miRNA gga-mir-2127 was identified to be targeting 3 genes–CD74, DUSP7 and PTPRC. Similarly, the miRNAs gga-mir-1647 (DUSP7 & GPRC5B), gga-mir-6556-5p (DDOST & LGMN), gga-mir-6561-5p (DUSP7 & IRF8), gga-mir-6587-3p (DUSP7 & LOC419851), gga-mir-6637-3p (ACTB & DUSP7), gga-mir-6649-3p (DDOST & TMEM263), gga-mir-6656-5p (DDOST & PDIA3), gga-mir-7456-5p (GPRC5B & TMEM263) and gga-mir-7460-5p (TF & TMEM263) were identified to be targeting 2 different genes each. Future research, including studies like small RNA-seq and chip-seq, can give further insights into miRNA and transcription factor expression levels during NDV challenge in chickens. This data, along with the information on transcription factors, microRNAs and co-expressing DElncRNAs identified in this study, can be used further to work on any specific genes that play a vital role in the immune system processes during NDV challenge in chickens to enhance their expression thereby improving the resistance in the susceptible breeds.

In our previous study on harderian gland data, it was observed that Leghorn showed a higher number of DEGs and DElncRNAs, especially those annotated with immune-related GO terms and pathways mainly in the comparison–Leghorn 6 DPC challenged vs. non-challenged. In comparison, resistant Fayoumi showed a lesser number of DEGs. The expression pattern of DElncRNAs was also found to be similar. In addition to this information, in the current study, Leghorn showed a higher number of upregulated immune-related genes than Fayoumi in timepoint-based analysis (between timepoints during the challenge), especially during the initial stages. This finding shows that Leghorn, being comparatively more susceptible to NDV than Fayoumi, showed several immune-related genes and positive co-expressing cis lncRNAs upregulated in challenged Leghorn when compared to non-challenged Leghorn and also in different timepoints during challenge. Surprisingly, Fayoumi, being comparatively resistant, showed little difference between challenged and non-challenged conditions and different time points of the challenge. This finding could also mean that Fayoumi might have similar immune gene and cis-lncRNA expression levels even during non-challenged conditions. The breed-based analysis, which compared Leghorn with Fayoumi in both challenged and non-challenged conditions separately, identified several immune-related genes and positive co-expressing cis lncRNAs to be upregulated in Fayoumi when compared to Leghorn in both challenged and non-challenged conditions. This result clearly shows that the differences in the expression of genes annotated with immune-related GO terms and pathways, i.e., immune-related genes and the co-expressing cis-lncRNAs between Leghorn and Fayoumi, and their role in the presence of differences in the resistance of Leghorn and Fayoumi chicken against NDV. This study is limited to in silico analysis of harderian gland tissue transcriptome data and experimental validation of expression of a few selected lncRNAs and co-expressing genes. The analysis of different tissues helps to understand the overall response of the host during disease. Further studies need to be done to unravel the mechanisms of co-expression and interaction of DEGs and DElncRNAs, which might be used for breed improvement.

Publicly available datasets were analysed in this study. This data can be found here: EBI-ENA database with project ID PRJEB22672. The pipeline used in the study – FHSpipe can be found here: https://github.com/Venky2804/FHSpipe.

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2024.1365888/full#supplementary-material