Thymic transcriptome analysis after Newcastle disease virus inoculation in chickens and the influence of host small RNAs on NDV replication

Newcastle disease (ND), caused by virulent Newcastle disease virus (NDV), is a contagious viral disease affecting various birds and poultry worldwide. In this project, differentially expressed (DE) circRNAs, miRNAs and mRNAs were identified by high-throughput RNA sequencing (RNA-Seq) in chicken thymus at 24, 48, 72 or 96 h post LaSota NDV vaccine injection versus pre-inoculation group. The vital terms or pathways enriched by vaccine-influenced genes were tested through KEGG and GO analysis. DE genes implicated in innate immunity were preliminarily screened out through GO, InnateDB and Reactome Pathway databases. The interaction networks of DE innate immune genes were established by STRING website. Considering the high expression of gga-miR-6631-5p across all the four time points, DE circRNAs or mRNAs with the possibility to bind to gga-miR-6631-5p were screened out. Among DE genes that had the probability to interact with gga-miR-6631-5p, 7 genes were found to be related to innate immunity. Furthermore, gga-miR-6631-5p promoted LaSota NDV replication by targeting insulin induced gene 1 (INSIG1) in DF-1 chicken fibroblast cells. Taken together, our data provided the comprehensive information about molecular responses to NDV LaSota vaccine in Chinese Partridge Shank Chickens and elucidated the vital roles of gga-miR-6631-5p/INSIG1 axis in LaSota NDV replication.


Results
In this text, differential expression profiles of circRNAs, miRNAs and mRNAs were examined by RNA-seq in chicken thymic tissues at 24 (B), 48 (C), 72 (D) or 96 (E) h after LaSota NDV vaccine inoculation compared to pre-inoculation (A) group, and corresponding results were presented in Supplementary Tables 1, 2 and 3, respectively.

Identification of DE mRNAs.
In addition, 989 up-regulated mRNAs and 397 down-regulated mRNAs were found in chicken thymus tissue samples at 24 h post vaccine injection relative to control group (Supplementary Table 6, Fig. 3A). The expression levels of 1016 mRNAs were notably up-regulated and expression   Table 11). GO annotation analyses showed that 19 (Supplementary Table 12 sheet 2), 18 (Supplementary Table 13 sheet 2), 21 (Supplementary Table 14 sheet 2) or 7 (Supplementary Table 15 sheet 2) biological process terms were significantly enriched by the DE genes in B vs A, C vs A, D vs A, or E vs A group, respectively. These biological processes or pathways might function as crucial players in the responses to NDV vaccine.
Screening of vaccine-influenced innate immune genes. Additionally, 32 genes were found to be implicated in innate immune responses by seeking for the terms containing the key word of innate immune in GO annotation data (Supplementary Table 12 Table 16 sheet 2).   Table 18, sheet 5) genes were found to be differentially expressed in B vs A, C vs A, D vs A or E vs A group, respectively. Moreover, 4 innate immune genes (CCAAT enhancer binding protein beta (CEBPB), NFKB inhibitor beta (NFKBIB), granulin precursor (GRN) and frizzled class receptor 1 (FZD1)) were markedly down-regulated and 6 innate immune genes (amphiregulin (AREG), complement C8 alpha chain (C8A), interleukin 22 (IL22), GATA binding protein 4 (GATA4), interleukin 1 receptor accessory protein like 1 (IL1RAPL1), and coagulation factor XI (F11)) were notably up-regulated in all the groups of B vs A, C vs A, D vs A, and E vs A groups (Supplementary Table 18, sheet 6). vs A, C vs A, D vs A, and E vs A groups, respectively. Pathway in the Y-axis represents the KEGG pathway term enriched by differentially expressed genes. The "rich" in the X-axis represents the rich factor. Rich factor: the number of differentially expressed genes in the pathway term/the number of all annotated genes in the pathway term. FDR: P value after the correction of False Discovery Rate (FDR). DEG_number: number of differentially expressed genes in the KEGG pathway term. KEGG enrichment analysis was performed as previously described 48 Table 20, sheet 5) genes were found to be differentially expressed in B vs A, C vs A, D vs A, or E vs A group, respectively. Furthermore, 3 common innate immunity-related genes (complement factor I (CFI), complement C8 alpha chain (C8A), and fyn related Src family tyrosine kinase (FRK)) were found to be strikingly up-regulated in all the groups of B Table 20 sheet 6). The information of these innate immunity-associated DE genes annotated by GO, Reactome pathway, and InnateDB databases was integrated into Supplementary Table 21.
Prediction of gga-miR-6631-5p potential targets. It is well known that miRNAs can regulate the expression of protein-coding genes at post-transcriptional levels 16 . Prediction analysis showed that 1431 genes had the potential gga-miR-6631-5p binding sites (Supplementary Table 25 sheet 1). Among these possible targets, expression profiles of 1426 items were analyzed in our experiments (Supplementary Table 25 sheet 2). Our outcomes presented that mRNA levels of 14 potential targets were notably down-regulated and mRNA levels of 50 potential targets were markedly up-regulated in B vs A group (Supplementary Table 25  MiR-6631-5p overexpression promoted LaSota NDV replication in DF-1 cells. Considering the high expression of gga-miR-6631-5p in thymic tissues of Chinese Partridge Shank Chickens across all four time points (24,48,72, and 96 h) post NDV vaccine treatment, we supposed that gga-miR-6631-5p might play vital roles in the responses to LaSota NDV. RT-qPCR assay also validated that gga-miR-6631-5p was highly expressed in LoSota-infected DF-1 cells compared to uninfected cells (Fig. 6A), which was in line with the RNA-seq data.
To further explore the effect of gga-miR-6631-5p on LaSota NDV replication and related molecular basis, gga-miR-6631-5p mimic and its control NC mimic were synthesized. Transfection efficiency analysis presented that the transfection of gga-miR-6631-5p mimic led to the notable up-regulation of gga-miR-6631-5p level in DF-1 cells than that in cells transfected with NC mimic (Fig. 6B), suggesting that gga-miR-6631-5p mimic could www.nature.com/scientificreports/ be used for subsequent gain-of-function experiments. Next, the effect of gga-miR-6631-5p overexpression on LaSota NDV replication was further examined by TCID 50 assay. Results showed that the enforced expression of gga-miR-6631-5p facilitated LaSota NDV replication in DF-1 cells (Fig. 6C). RT-qPCR assay revealed that gga-miR-6631-5p overexpression led to the remarkable down-regulation of INSIG1, HSPBP1, SERPING1 and PSMD11 mRNA levels, but had no much effect on E2F1 mRNA expression (Fig. 7A). In view of the strongest inhibitory effect of gga-miR-6631-5p on INSIG1 mRNA expression, luciferase reporter assay was performed to further explore whether INSIG1 was a target of gga-miR-6631-5p in DF-1 cells. Results showed that gga-miR-6631-5p overexpression led to the notable down-regulation of luciferase activity of INSIG1-wt reporter, but had no much influence on luciferase activity of INSIG1-mut reporter with mutant gga-miR-6631-5p binding site (Fig. 7B,C). These data revealed that gga-miR-6631-5p could bind with INSIG1 3'UTR by predicted sites. Transfection efficiency analysis revealed that the transfection of in-miR-6631-5p could notably reduce gga-miR-6631-5p level in DF-1 cells compared with in-NC control group (Fig. 7D). The transfection of si-INSIG1 led to a dramatic decrease in INSIG1 mRNA expression in DF-1 cells relative to si-con control group (Fig. 7E). TCID 50 assay revealed that the depletion of gga-miR-6631-5p led to a noticeable reduction in LaSota NDV replicative potential in DF-1 cells (Fig. 7F). In addition, INSIG1 knockdown notably abrogated the detrimental effect of gga-miR-6631-5p loss on LaSota NDV replication in DF-1 cells (Fig. 7F). In a word, these outcomes suggested that INSIG1 was a target of gga-miR-6631-5p and gga-miR-6631-5p promoted LaSota NDV replication through down-regulating INSIG1 in DF-1 cells.  Additionally, 32 DE genes were found to be implicated in innate immune responses based on the GO annotation data. Considering the conservation of innate immune responses among different organisms 25 , more innate immunity-related genes were screened out through InnateDB and Reactome Pathway databases. Integration analysis showed that a total of 103 innate immune genes annotated by GO, Reactome pathway, and InnateDB databases were differentially expressed in response to NDV infection. Among these innate immunity genes, 4 genes (CEBPB, NFKBIB, GRN and FZD1) were low expressed and 8 genes (CFI, C8A, FRK, AREG, IL22, GATA4, IL1RAPL1 and F11) were highly expressed in all the groups of B vs A, C vs A, D vs A, and E vs A. AREG, a member of the epidermal growth factor (EGF) family, has been found to be expressed by multiple innate immune cell populations (e.g. basophils, mast cells, neutrophils) and adaptive immune cell populations (CD4 + T cells, CD8 + T cells) besides epithelial and mesenchymal cells 26 . Moreover, AREG was reported to be involved in mediating pathogen resistance, orchestrating tissue repair and homeostasis, inhibiting local inflammation and regulating immune responses 26 . Qiao et al. demonstrated that GATA4 loss inhibited inflammatory responses and NF-κB activation induced by lipopolysaccharide in human dental pulp cells 27 . Also, Kang et al. presented that GATA4 was involved in the regulation of inflammatory responses 28 .

The introduction of gga-miR
Furthermore, the interaction network of proteins encoded by DE innate immune genes was established through STRING website. KEGG enrichment analyses showed that these DE innate immune genes mainly participated in the regulation of 9 KEGG pathways including nucleotide-binding oligomerization domain-containing (NOD)-like receptor signaling pathway, retinoic acid-inducible gene I (RIG-I)-like receptor signaling pathway, cytokine-cytokine receptor interaction and phagosome pathway. Some researchers proposed that phagosome might function as the link organelle between innate and adaptive immunity 29,30 , suggesting the vital roles of phagosome pathway in immune responses. It has been well documented that RNA virus can be recognized by Toll-like receptors, RIG-I-like receptors and NOD-like receptors, whose activation in response to viral RNAs can initiate innate immune responses such as pro-inflammatory cytokine secretion and type I interferon (IFN) production [31][32][33] . Moreover, Fournier et al. pointed out that the intranasal application of NDV MTH68 strain could induce pro-inflammatory cytokine secretion in mouse bronchial lavage fluid and activate RIG-I and IFN pathways in mouse dendritic cells 34 . In addition, previous studies showed that RIG-I and IFN pathways played crucial roles in protecting host cells from NDV infection 35,36 .
Despite the discovery of multitudinous circRNAs over the past decades, the functions and regulatory mechanisms of most circRNAs remain unknown 37 . Recently, accumulating evidence shows that circRNAs can function as miRNA sponges (ceRNAs) to relieve the inhibitory activity of miRNAs on target mRNAs by sequestering miRNAs 38,39 . Considering the abnormal expression of gga-miR-6631-5p across all the time points post NDV vaccine infection, circRNAs or mRNA targets that had a possibility to interact with gga-miR-6631-5p were predicted using the miRanda software. Combined with expression data of circRNAs/mRNAs and prediction data of circRNAs/mRNAs and gga-miR-6631-5p, 6 DE circRNAs and 64 DE genes in B vs A group, 8 DE circRNAs and 73 DE genes in C vs A group, 1 DE circRNA and 64 DE genes in D vs A group, and 2 DE circRNAs and 33 DE genes in E vs A group were found to have the possibility to bind to gga-miR-6631-5p. Among these DE gga-miR-6631-5p targets, 5 genes (E2F1, INSIG1, HSPBP1, MAPK10, MASP2) in B vs A group, 5 genes (SERPING1, E2F1, HSPBP1, MAPK10, MASP2) in C vs A group, 4 genes (E2F1, HSPBP1, PSMD11, INSIG1) in D vs A group and 3 genes (SERPING1, MAPK10, MASP2) in E vs A group were found to be implicated in innate immune responses.
Moreover, our data revealed that gga-miR-6631-5p expression was markedly up-regulated in DF-1 cells infected with LaSota NDV virus compared with uninfected cells. Enforced expression of gga-miR-6631-5p facilitated LaSota NDV replication through targeting INSIG1 in DF-1 cells. INSIG1 has been found to be a bridge in the activation of TANK-binding kinase 1 (TBK1) by promoting K27-linked polyubiquitination of stimulator of interferon genes (STING) and INSIG1 loss curbed STING-mediated antiviral gene induction 40 . TBK1 was a vital kinase in the induction of antiviral innate immune responses 41 and STING was a negative regulator in the replication of multiple RNA viruses such as Vesicular Stomatitis virus and Sindbis virus 42 . In addition, Ran et al. pointed out that STING knockdown suppressed the expression of interferon regulatory factor 7 (IRF7), type I interferon (IFN)-α, and IFN-β in chicken embryo fibroblasts 43 .
The previous transcriptome analyses in embryos 44 , spleen 45 or trachea 46 from Leghorn and Fayoumi chickens and chicken embryo fibroblasts 47 mainly explored the influence of NDV vaccine treatment on immune cells, molecular pathways, and gene expression. In our project, we investigated the expression profile alterations of circRNAs, miRNAs and mRNAs in thymic tissues of Chinese Partridge Shank Chickens at 24, 48, 72, 96 h post NDV infection versus pre-inoculation group. In addition, we further analyzed the expression patterns of 712 innate immune genes annotated by InnateDB database and 319 innate immune genes annotated by Reactome Pathway database in thymic tissues of Chinese Partridge Shank Chickens. A total of 103 innate immune genes annotated by GO, Reactome pathway, and InnateDB databases were identified to be differentially expressed in response to NDV infection in our study. Moreover, we established the potential interaction network of proteins encoded by these dysregulated innate immune genes through String website. Additionally, some possible circR-NAs/gga-miR-6631-5p and gga-miR-6631-5p/genes pairs were constructed based on bioinformatics prediction analysis. Furthermore, our data revealed that gga-miR-6631-5p promoted LaSota NDV replication by targeting INSIG1 in DF-1 cells. These data might contribute to the investigation on interactions of NDV LaSota vaccine

Conclusion
In summary, our data provided the time-specific comprehensive information of circRNAs, miRNAs and mRNAs in response to NDV LaSota vaccine infection in thymic tissues of Chinese Partridge Shank chickens and elucidated the vital roles of gga-miR-6631-5p/INSIG1 axis in LaSota NDV replication. Also, multiple innate immune genes implicated in NDV LaSota vaccine infection were screened out.

Materials and methods
Animal experiment. All experimental procedures were approved by Institutional Animal Care and Use Committee of Henan University of Animal Husbandry and Economy. Specific pathogen-free (SPF) Chinese Partridge Shank Chickens (30-days-old) free of NDV antibody were reared in the rooms at the biosafety level II facility. These chickens were inoculated with or without 10 5 50% egg-infectious dose (EID 50  Data processing. Raw sequencing data were analyzed and filtered to obtain high-quality clean data, which were then aligned to the reference genome of Gallus_gallus.Gallus_gallus-5.0 (assembly GCA_000002315.3) using Tophat2.
Differential expression analysis. Expression analysis was conducted using the method of fragments per kilobase of transcript per million mapped reads (FPKM). Genes were considered as DE when the absolute value of fold change was greater than 2 and P value was less than 0.05. Reagents and cell transfection. The INSIG1-specific small interference RNAs (si-INSIG1) and its negative control si-con were designed and synthesized from Shanghai GenePharma Co., Ltd. (Shanghai, China). The gga-miR-6631-5p mimic and NC mimic control, gga-miR-6631-5p inhibitor (in-miR-6631-5p) and its negative control in-NC were purchased from Thermo Scientific Co., Ltd. Cell transfection was performed using Lipofectamine 3000 Reagent (Thermo Scientific) following the instructions of manufacturer.

Gene annotation analysis.
RT-qPCR assay. RNA was isolated and quantified as described above. cDNA first stands were synthesized from RNA template using the iScript cDNA Synthesis Kit (Bio-Rad, Hercules, CA, USA) and random primer (for mRNAs) or stem-loop gga-miR-6631-5p RT primer (5'-CTC AAC TGG TGT CGT GGA GTC GGC AAT TCA GTT GAG GCA GAA CC-3'). Next, real-time quantitative PCR reaction was performed using ITaq Universal SYBR Green Super mix (Bio-Rad) and quantitative primers. Relative expression levels of gga-miR-6631-5p and mRNAs were calculated using the method of 2 −ΔΔCT . GAPDH or 5S ribosomal RNA (5S rRNA) functioned as the endogenous inference to normalize the expression of mRNAs or gga-miR-6631-5p, respectively. Quantitative PCR primers were presented as follows: Statistics analysis. Statistics analysis was performed on GraphPad Prism 7 software (GraphPad Software, Inc., San Diego, CA, USA) with P < 0.05 as statistically significant. Difference between groups was compared using Student's T test. Results were presented as mean ± standard deviation (SD).
Ethical approval. All methods were carried out in accordance with relevant guidelines and regulations.

Data availability
The data displayed in this manuscript is available from the corresponding author upon reasonable request.