Transcriptome Analysis of Thymus Tissues from Chinese Partridge Shank Chickens Before/After NDV LaSota Vaccine Injection Via High-Throughput RNA Sequencing

Background: Newcastle disease (ND), caused by virulent Newcastle disease virus (NDV), is a huge threat for poultry and birds. LaSota strain of NDV is a common live attenuated vaccine to control ND. In this study, high-throughput RNA sequencing was performed to explore thymus tissue transcriptome change and better manage “vaccine failure” in Chinese Partridge Shank Chickens at 0 h and 48 h post LaSota vaccine injection. Results: 140 long non-coding RNAs (lncRNAs), 8 microRNAs (miRNAs) and 1514 mRNAs were differentially expressed post LaSota NDV vaccine inoculation. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis for vaccine-affected mRNAs and miRNA target genes was conducted. Moreover, 70 innate immune genes (31 up-regulated, 39 down-regulated) inGO, InnateDB and Reactome Pathway databases were found to be differentially expressed. Furthermore, the interaction network of proteins encoded by these dysregulated innate immune genes with known names was established by STRING website. Additionally, expression patterns of some dysregulated lncRNAs, miRNAs and mRNAs were further tested through RT-qPCR assay. Pearson correlation analysis for dysregulated transcripts validated by RT-qPCR assay was also performed. Based on the correlation data, 2 potential competing endogenous RNA (ceRNA) networks (ENSGALT00000060887/gga-miR-6575-5p/paraoxonase 1 (PON1) or MSTRG.188121.10/gga-miR-6631-5p/MMP9) were established. Conclusions: Our data presented many differentially expressed lncRNAs, miRNAs, mRNAs and innate immune genes, and established 2 potential ceRNA networks (lncRNAs/miRNAs/mRNAs) in Chinese Partridge Shank Chickens in response to LaSota vaccine inoculation, which could deepen our understanding on host responses to NDV LaSota vaccine in Chinese Partridge Shank Chickens.

Coding RNAs and non-coding RNAs including long non-coding RNAs (lncRNAs) and microRNAs (miRNAs) have been emerged as crucial players in various biological process and host responses to disease and pathogens in animals including poultry [9][10][11]. Recently, high-throughput transcriptome sequencing (RNA-Seq) has attracted much attention from researchers due to its outstanding advantages such as little genomic sequence limitation and low background signal [12]. RNA-seq is a powerful technique that can simultaneously capture almost all coding and noncoding transcripts and decipher the information of transcriptome including lncRNAs, miRNAs, and mRNAs [13,14]. RNA-seq not only can measure gene expression and discover novel RNAs especially non-coding RNAs, but also can identify differentially expressed genes under different conditions and investigate the host-pathogen interactions, which contributes to deepen our understanding on molecular pathology and help us to assess the changes associated with diseases [15,16].
In this text, we took advantage of RNA-seq technology to have an in-depth exploration on transcriptome changes in response to NDV LaSota vaccine injection in Chinese Partridge Shank Chickens. Our data revealed that 140 long noncoding RNAs (lncRNAs), 8 microRNAs (miRNAs) and 1514 mRNAs were differentially expressed (|log 2 FoldChange| > 1, P value < 0.05) in thymus tissue samples of Chinese Partridge Shank chickens post LaSota NDV vaccine inoculation. Moreover, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis for vaccineaffected mRNAs and miRNA target genes was conducted. Additionally, expression patterns of 712 innate immune genes in InnateDB database and 319 innate immune genes in Reactome Pathway Database were measured in our experiments. Moreover, 70 innate immune genes (31 up-regulated, 39 down-regulated) were demonstrated to be differentially expressed.
Also, expression patterns of some lncRNAs, miRNAs and mRNAs were further measured through RT-qPCR assay.

Sample collection
Speci c pathogen-free (SPF) Chinese Partridge Shank Chickens were obtained from Beijing Merial Laboratory Animal Technology Co. Ltd (Beijing, China) and raised in the rooms at the biosafety level II facility. Chickens (n=20) were randomly divided into 2 groups (pre-inoculation and post-inoculation groups) with 10 chickens in each group. At the age of 30 days, the chickens in the post-inoculation group were inoculated with 0.2 ml of 10 5 50% egg-infectious dose (EID 50 ) Lasota suspension through eyes, respectively. Chickens were euthanized by the intravenous injection of pentobarbital sodium solution (90 mg/kg body weight) in each group. The thymus tissue samples were collected in RNase-free microtubes, immediately frozen in liquid nitrogen, and then stored at -80℃ until RNA extraction.
Our study was approved by Animal Care and Use Committee of the Henan University of Animal Husbandry and Economy and was performed in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of Henan University of Animal Husbandry and Economy. All efforts were made to minimize suffering of birds.
Total RNA extraction, cDNA Library Construction and Sequencing RNA was isolated from thymus tissue samples using Trizol Reagent (Thermo Fisher Scienti c, Waltham, MA, USA) and then treated with a TURBO DNA-free™ Kit (Thermo Fisher Scienti c) to remove DNA from RNA samples (Thermo Fisher Scienti c). Next, the concentration and purity of RNA was measured using NanoDrop 2000 spectrophotometer (Thermo Fisher Scienti c) and RNA quality was further tested using a RNA Nano 6000 Assay Kit (Agilent Technologies, Santa Clara, CA, USA) through Agilent Bioanalyzer 2100 (Agilent Technologies). After the removal of ribosomal RNA using the Ribo-Zero rRNA Removal Kit (Illumina), high-quality samples (RNA concentration≥400 ng/μl, OD260/280: 1.8-2.2, RNA Integrity number≥8) were used for the construction of following cDNA libraries of mRNAs and lncRNAs according to TruSeq RNA Sample Preparation Guide (Illumina, San Diego, CA, USA). Brie y, mRNAs and lncRNAs were enriched using magnetic beads and then fragmented using divalent cations. Next, fragmented mRNAs were reversely transcribed into cDNAs, followed by the conversion of residual overhangs into blunt ends. Thereafter, Illumina PE adapter oligonucleotides were ligated with adenylated DNAs (3'ends). Subsequently, cDNA fragments were enriched by PCR ampli cation and puri ed through Agencourt AMPure XP System (Beckman Coulter, Beverly, CA, USA). Next, cDNA libraries were quanti ed by the Agilent high sensitivity DNA assay on Agilent Bioanalyzer 2100 (Agilent Technologies). Finally, signal strand cDNA libraries were sequenced on Illumina HiSeq 2500 instrument (Illumina) by Shanghai Personal Biotechnology Co. Ltd (Shanghai, China). Small RNA libraries were conducted through NEBNext Multiplex Small RNA Library Prep Set for Illumina (New England Biolabs, Ipswich, MA, USA) following the protocols of manufacturer. Brie y, RNA samples were ligated to RNA 3' and 5' adapters and then reverse-transcribed and ampli ed into cDNA libraries. After quanti ed and quali ed by Agilent Bioanalyzer 2100 (Agilent Technologies), libraries were sequenced.

Raw data processing
The information of reference genome of Gallus_gallus.Gallus_gallus-5.0 (assembly GCA_000002315.3) established by Ensembl database (version 94.5) was shown in Table 1. Also, we provided the annotation information of Gallus_gallus.Gallus_gallus-5.0 (assembly GCA_000002315.3) genes in different databases including GO, KEGG and UniProtID databases ( Table 2). The information of raw data was presented in Table 3. The information of clean data without reads containing 3' adaptors or average base quality < Q20 was provided in Table 4. After ltering, high-quality clean data were aligned to reference genome of Gallus_gallus.Gallus_gallus-5.0 (assembly GCA_000002315.3) using Tophat2. The alignment results were presented in Table 5.

Differential Expression analysis
Expression levels of lncRNAs and mRNAs were measured using the fragments per kilobase of transcript per million mapped reads (FPKM) method. Differentially expressed genes were selected using the condition of |log 2 FoldChange| > 1 and Pvalue<0.05.

Target prediction
The potential interaction between miRNAs and lncRNAs or mRNAs was predicted by miRanda software.
Reverse transcription-quantitative PCR (RT-qPCR) RNA was extracted using Trizol Reagent (Thermo Fisher Scienti c) and quanti ed using NanoDrop 2000 spectrophotometer (Thermo Fisher Scienti c) as described above. Then, RNA was reversely transcribed into cDNA rst strands using an iScript cDNA Synthesis Kit (Bio-Rad, Hercules, CA, USA) and random primer (for mRNAs and lncRNAs) or stem-loop miRNA RT primers in Table 6. Next, cDNA was ampli ed using ITaq Universal SYBR Green Supermix (Bio-Rad) and speci c primers in Table 7. Expression levels of lncRNAs, miRNAs and mRNAs were calculated using the formula of 2 −ΔΔCT with GAPDH or 5S RNA as the internal control.

Statistics Analysis
Data was analyzed using GraphPad Prism 6 software (GraphPad Software, Inc., San Diego, CA, USA). Student's T test was used to compare the difference between groups. Difference was regarded as statistically signi cant when P value was less than 0.05.

Identi cation 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) (|log 2 FoldChange| > 1, P value<0.05) were signi cantly altered in thymus tissue samples of chickens after NDV vaccine treatment (Supplementary Table 1, Figure 1).
Based on the outcomes of GO enrichment analysis, vaccine-in uenced 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  Table 9, Sheet 2).
In view of conservation of innate immune responses among different organisms, we screened out more innate-immunityrelated 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   Establishment of interaction networks among proteins encoded by these ltered innate immune genes Next, the interactions among proteins encoded by ltered 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 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.
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.

Discussion
Vaccination has been well documented as an effective preventive strategy to protect animals from infectious diseases by stimulating the immune system [17,18]. To better utilize vaccines, it is requisite to gain a comprehensive knowledge and understanding on all possible vaccination effects (e.g. mode of action, risks, advantages, and disadvantages). In this study, transcriptome analysis by RNA-seq technology was performed to more deeply investigate the cellular and molecular responses to LaSota NDV vaccine in Chinese Partridge Shank Chickens.
In this text, expression patterns of 6491lncRNAs, 17131 mRNAs and837 miRNAs were examined in thymus tissue samples of Chinese Partridge Shank Chickens before/after LaSota vaccine injection by RNA-seq and miRNA-sEq. Among these transcripts, expression of 58 lncRNAs, 3 miRNAs and 1016 mRNAs were notably up-regulated and expression of 82 lncRNAs, 5 miRNAs and 498 mRNAs was markedly down-regulated in thymus tissue samples of Chinese Partridge Shank Chickens post LaSota NDV vaccine inoculation compared with pre-inoculation group. Moreover, GO and KEGG analysis for vaccine-affected mRNAs and putative target genes of dysregulated miRNAs was conducted in our project, which was presented in Supplementary Table 2B-2D, Supplementary Table 3B and Supplementary Table 3C.
Recently, ceRNA hypothesis has attracted much attention from researchers in elucidating the molecular basis of noncoding RNAs including lncRNAs [39]. ceRNA hypothesis proposes that lncRNAs can perform as the molecular sponge of miRNAs to attenuate the inhibitory effect of miRNAs on target genes [40,41]. Hence, co-expression patterns of these dysregulated lncRNAs, miRNAs and mRNAs validated by RT-qPCR assay were explored in our project. Pearson correlation analysis presented that gga-miR-6631-5p expression was negatively associated with the expression of ENSGALT00000065826 (lncRNA), MSTRG.188121.10 (lncRNA), or MMP9 in thymus tissue samples of chickens after NDV vaccine injection. Moreover, MSTRG.188121.10 expression was positively related with MMP9 expression in post-challenge group. These data suggested that MSTRG.188121.10 might act as a ceRNA of gga-miR-6631-5p to regulate MMP9 expression in chickens post NDV vaccine inoculation and there existed a possible link between ENSGALT00000065826 and gga-miR-6631-5p. Also, an inverse correlation between MSTRG.188121.10 and CFTR expression was found in postinoculation group, suggesting the close link of MSTRG.188121.10 and CFTR. In addition, gga-miR-6575-5p expression was inversely correlated with ENSGALT00000060887 (lncRNA), MSTRG.188121.10 (lncRNA) or PON1 expression in thymus tissue samples of chickens post NDV vaccine challenge. ENSGALT00000060887 expression was positively associated with PON1 expression in post-inoculation group. These outcomes suggested that ENSGALT00000060887 could regulate PON1 expression by functioning as a molecular sponge of gga-miR-6575-5p and there was a possible negative regulatory relationship between MSTRG.188121.10 and gga-miR-6575-5p. Furthermore, ENSGALT00000060887 expression was positively associated with MMP9 expression in post-challenge group, suggesting the relevance of ENSGALT00000060887 and MMP9 in chickens inoculated with NDV vaccine.
Although previous studies have performed the transcriptome analysis to explore host immune responses against NDV LaSota vaccine infection in embryos [42], spleen [43]