Identification of differentially expressed lncRNAs
The effect of NDV vaccine on expression of 6491 lncRNAs was examined by RNA-seq in chicken thymus tissue samples. Results showed that expression levels of 140 lncRNAs (58 up-regulated and 82 down-regulated) (|log2FoldChange| > 1, P value<0.05) were significantly altered in thymus tissue samples of chickens after NDV vaccine treatment (Supplementary Table 1, Figure 1).
Identification of differentially expressed miRNAs and bioinformatics analysis of miRNA targets
Also, expression patterns of 837 miRNAs were measured in our project. Eight miRNAs (3 up-regulated, 5 down-regulated) (|log2FoldChange| > 1, P value<0.05) were found to be differentially expressed in thymus tissue samples isolated from chickens post vaccine treatment compared to control group (Supplementary Table 2, Figure 2A). In detail, expression of gga-miR-6631-5p, gga-miR-6575-5p and gga-miR-1727 was remarkably up-regulated and expression levels of gga-miR-1712-3p, gga-miR-12273-3p, gga-miR-6655-5p, gga-miR-1744-3p andgga-miR-6548-3p were notably down-regulated in thymus tissue samples of chickens post vaccine injection (Supplementary Table 2).It is well known that miRNAs can exert their function through regulating specific mRNAs. Hence, potential targets of differentially expressed miRNAswere examined by bioinformatics analysis using the miRanda software and the results were presented in Supplementary Table 3. To have a deep insight into the function of these miRNAs, GO and KEGG enrichment analysis for these miRNA target genes was carried out. The top 10 KEGG pathways enriched by these differentially expressed miRNA targets were Endocytosis (ko04144), Collecting duct acid secretion (ko04966), Epstein-Barr virus infection (ko05169), Viral carcinogenesis (ko05203), Peroxisome (ko04146), Plant-pathogen interaction (ko04626), Bacterial invasion of epithelial cells (ko05100), Longevity regulating pathway-multiple species (ko04213), Other types of O-glycan biosynthesis (ko00514), and Antigen processing and presentation (ko04612) (Supplementary Table 4 and Figure 2B). GO analysis for putative targets of differentially expressed miRNAs suggested that miRNA target genes were involved in the regulation of multiple biological processes, especially for cellular metabolic process (GO:0044237), metabolic process (GO:0008152), organic substance metabolic process (GO:0071704), primary metabolic process (GO:0044238), nitrogen compound metabolic process (GO:0006807), biosynthetic process (GO:0009058), organic substance biosynthetic process (GO:1901576), cellular biosynthetic process (GO:0044249), cellular nitrogen compound metabolic process (GO:0034641) and small molecule metabolic process (GO:0044281) (Top 10 GO terms in biological processes) (Supplementary Table 5, Sheet 1). Molecular function analysis in GO annotation showed that these miRNA target genes mainly exerted catalytic activity and binding activity with multiple molecules such as protein, enzyme, and nucleoside (Supplementary Table 5, Sheet 2).
Identification of differentially expressed mRNAs and related GO and KEGG enrichment analysis
Moreover, 1016 up-regulated mRNAs (Supplementary Table 3A, sheet 1) and 498 down-regulated mRNAs (Supplementary Table 3A, sheet 2) (|log2FoldChange| > 1, P value<0.05) were identified among 17131 mRNAs in thymus tissue samples of chickens at 48 h after NDV vaccine infection (Supplementary Table 6, Figure 3A). The top 10 up-regulated mRNAs with known symbols were cyclin dependent kinase 3(CDK3), thyrotropin releasing hormone receptor(TRHR), mesogenin 1 (MSGN1), chromosome 1 open reading frame 158 (C1ORF158), actin like 6B(ACTL6B), G protein subunit gamma transducin 2 (GNGT2), epidermal differentiation protein starting with MTF motif 4 (EDMTF4), myogenic factor 5(MYF5), avian beta-defensin 2 (AvBD2) and meiotic double-stranded break formation protein 4(MEI4)(Supplementary Table 6, sheet 2). The top 10 down-regulated mRNAs with known names were family with sequence similarity 181 member B(FAM181B), phosphatidylinositol glycan anchor biosynthesis class C (PIGC), ubiquitin like 4A (UBL4A), inka box actin regulator 1(INKA1), spermatogenesis associated 2 like (SPATA2L), NADH: ubiquinone oxidoreductase subunit B7 (NDUFB7), cyclin dependent kinase 2 associated protein 2 (CDK2AP2),adaptor related protein complex 5 subunit sigma 1(AP5S1), cytochrome c oxidase assembly factor COX14(COX14) and apolipoprotein L domain containing 1 (APOLD1) (Supplementary Table 3A, sheet 3).
Based on the outcomes of GO enrichment analysis, vaccine-influenced genes participated in the regulation of many biological processes, especially for transmembrane transport (GO:0055085), ion transport (GO:0006811), ion transmembrane transport (GO:0034220), inorganic ion transmembrane transport (GO:0098660), cation transport (GO:0006812), cation transmembrane transport (GO:0098655), calcium ion transmembrane transport (GO:0070588), calcium ion transport (GO:0006816), inorganic cation transmembrane transport (GO:0098662), metal ion transport (GO:0030001) (Top 10 terms in GO biological_process category) (Supplementary Table 3B, Sheet 2). Moreover, vaccine injection led to the abundant enrichment of genes belonging to cellular components such as nucleosome (GO:0000786), extracellular region (GO:0005576), DNA packaging complex (GO:0044815), integral component of plasma membrane (GO:0005887), intrinsic component of plasma membrane (GO:0031226), transporter complex (GO:1990351), transmembrane transporter complex (GO:1902495), protein-DNA complex (GO:0032993), plasma membrane part (GO:0044459) and ion channel complex (GO:0034702) (Top 10 terms in GO cellular_components category) (Supplementary Table 7, Sheet 3). Meanwhile, vaccine-affected genes were found to be involved in the regulation of molecular functions such as transmembrane transporter activity (GO:0022857), transporter activity (GO:0005215), inorganic molecular entity transmembrane transporter activity (GO:0015318), ion transmembrane transporter activity (GO:0015075), substrate-specific channel activity (GO:0022838), ion channel activity (GO:0005216), channel activity (GO:0015267), passive transmembrane transporter activity (GO:0022803), inorganic cation transmembrane transporter activity (GO:0022890), cation transmembrane transporter activity (GO:0008324) (Top 10 terms in GO molecular_function category) (Supplementary Table 3B, Sheet 4).KEGG analysis revealed that the injection of NDV vaccines led to differential expression of many genes, especially for those in Systemic lupus erythematosus (ko05322), Neuroactive ligand-receptor interaction (ko04080), Retrograde endocannabinoid signaling (ko04723), Staphylococcus aureus infection (ko05150), Complement and coagulation cascades (ko04610), Nicotine addiction (ko05033), Alcoholism (ko05034), Glutamatergic synapse (ko04724), GABAergic synapse (ko04727), and Asthma (ko05310) pathways (Top 10 KEGG pathway) (Supplementary Table 8, Figure 3B).
Screening and analysis of vaccine-influenced innate immune genes
Furthermore, differentially expressed innate immune genes were presented in Sheet 1 of Supplementary Table 4A. The information of these genes was integrated and shown in Sheet 3 of Supplementary Table 9. The innate immune genes with known names were transmembrane protein 173 (TMEM173), RNA binding motif protein 14 (RBM14), complement factor D (CFD), solute carrier family 11 member 1(SLC11A1), mannan binding lectin serine peptidase 2(MASP2), cytochrome b-245 alpha chain (CYBA), apolipoprotein A4 (APOA4), macrophage migration inhibitory factor (glycosylation-inhibiting factor) (MIF), retinoic acid receptor responder (tazarotene induced) 2 (RARRES2), collectin subfamily member 11 (COLEC11), tripartite motif containing 62 (TRIM62), otopetrin 1 (OTOP1), triokinase and FMN cyclase (TKFC) and serpin family G member 1 (SERPING1) (Supplementary Table 9, Sheet 2).
In view of conservation of innate immune responses among different organisms, we screened out more innate-immunity-related genes in InnateDB and Reactome Pathway databases. In this text, expression patterns of 712 innate immune genes in the InnateDB database (http://innatedb.sahmri.com/annotatedGenes.do?type=innatedb, Supplementary Table 10) were measured(Supplementary Table 4C sheet 1). Moreover, 50 innate immune genes(19 up-regulated, 31 down-regulated) were found to be differentially expressed in experimental group vs control group (Supplementary Table 11 sheet 2).Moreover, expression changes of 319innate immune genes in the Reactome Pathway Database (http://reactome.ncpsb.org/PathwayBrowser/#/R-GGA-168249&PATH=R-GGA-168256&DTAB=MT,Supplementary Table 12) were examined in our project (Supplementary Table 4E sheet 1). Among these genes, 15 innate immune genes (11 up-regulated, 4 down-regulated) were differentially expressed in experimental group vs control group (Supplementary Table 13 sheet 2). Finally, the information of 70 innate immune genes annotated by GO, Reactome pathway, and InnateDB databases (31 up-regulated, 39 down-regulated) was integrated into Supplementary Table 14.
Establishment of interaction networks among proteins encoded by these filtered innate immune genes
Next, the interactions among proteins encoded by filtered innate immune genes were analyzed through STRING: functional protein association networks (https://string-db.org/cgi/network.pl?taskId=s7sL6l9NyZlM). The image of interaction network was shown in Figure 4. The interaction edges with the combined score ≥ 0.4 were presented in sheet 1 of Supplementary Table 5A and related annotation information was displayed in sheet 2 of Supplementary Table 15. Moreover, KEGG function annotation analysis revealed that cytokine-cytokine receptor interaction pathway (matching proteins: Interleukin-10 (IL10), Interleukin-22 (IL22), Interleukin-5 (IL5), Interleukin-9 (IL9), prolactin (PRL), phagosome pathway (matching proteins: ATPase H+ transporting V0 subunit e2 (ATP6V0E2), collectin subfamily member 11 (COLEC11), cytochrome b-245 alpha chain (CYBA),nitric oxide synthase 1 (NOS1)) and Retinoic acid‑inducible gene I (RIG‑I)-like receptor signaling pathway (matching proteins: triokinase and FMN cyclase (TKFC/DAK), mitogen-activated protein kinase 10 (MAPK10),TMEM173 were significantly enriched by these screened innate immune genes (Supplementary Table 16 sheet 1). In addition, these innate immune genes led to the significant enrichment of 3 Cellular Component terms (Supplementary Table 16 sheet 2) and 41 Biological Process terms (Supplementary Table 16 sheet 3) in GO analysis. In addition, GO analysis for these dysregulated innate immune genes showed that CCAAT enhancer binding protein beta (CEBPB) and IL10 were involved in the negative regulation of neuron death.
RT-qPCR validation of some differentially expressed lncRNAs, miRNAs and mRNAs
Next, RT-qPCR validated that expression levels of MSTRG.22689.1, ENSGALT00000065826, ENSGALT00000059336 and ENSGALT00000060887 were notably up-regulated and MSTRG.188121.10expression was remarkably down-regulated in thymus tissue samples of 10 random chickens at 48 h post NDV challenge compared to pre-challenge group (Figure 5A). Also, we demonstrated that gga-miR-6575-5p, gga-miR-6631-5p and gga-miR-1727 were highly expressed, while gga-miR-6655-5p and gga-miR-6548-3p were low expressed in thymus tissue samples of 10 random chickens at 48 h after NDV vaccine inoculation relative to pre-inoculation group (Figure 5B). In addition, RT-qPCR results substantiated that paraoxonase 1 (PON1), mitogen-activated protein kinase 10 (MAPK10) and cystic fibrosis transmembrane conductance regulator (CFTR) were highly expressed and matrix metallopeptidase 9 (MMP9) were low expressed in thymus tissue samples of 10 random chickens post NDV vaccine injection relative to pre-injection group (Figure 5C). However, there was no conspicuous difference inthe expression of gga-miR-1744-3p (Figure 5B), gga-miR-1712-3p (Figure 5B), E2F transcription factor 1 (E2F1) (Figure 5C) and nuclear factor kappa B subunit 2(NF-κB2) (Figure 5C) between pre-inoculation and post-inoculation groups.
Co-expression analysis of dysregulated lncRNAs, miRNAs and mRNAs
Next, Pearson correlation analysis revealed that gga-miR-6631-5p expression was negatively associated with the expression of ENSGALT00000065826 or MSTRG.188121.10 in thymus tissue samples of chickens after NDV vaccine injection (Figure 6A). Moreover, there was a negative correlation between gga-miR-6575-5p expression and ENSGALT00000060887 or MSTRG.188121.10 expression in thymus tissue samples of 10 random chickens post NDV vaccine challenge (Figure 6A). In addition, MSTRG.188121.10 expression was positively related with MMP9 expression, but was inversely correlated with CFTR expression in thymus tissue samples of chickens following NDV vaccine treatment (Figure 6B). ENSGALT00000060887 expression was found to be positively associated with PON1 or MMP9 expression in thymus tissue samples of chickens after NDV vaccine administration (Figure 6B).Next, the correlation of mRNAs and miRNAs with shared related lncRNAs was further analyzed. Results showed that there was a negative correlation between gga-miR-6631-5p and MMP9 expression in thymus tissue samples of chickens post NDV vaccine administration (Figure 6C). Also, gga-miR-6575-5p expression was found to be inversely related withPON1 expression in thymus tissue samples of chickens post NDV vaccine challenge (Figure 6C). Based on the co-expression data of ENSGALT00000060887, gga-miR-6575-5p and PON1, we supposed that ENSGALT00000060887 might function as a molecular sponge of gga-miR-6575-5p to sequester gga-miR-6575-5p from PON1in chickens post NDV vaccine injection. MSTRG.188121.10 might regulate MMP9 expression through acting as a competing endogenous RNA (ceRNA) of gga-miR-6631-5p in chickens post NDV vaccine injection.