Transcriptome Analysis Reveals the Immune Response of Chicken Erythrocytes to Marek’s Disease Virus Infection

Guan-bao Ning1, Sheng Niu1, Yue-jian Li1,2, Xiao-xiao Lu1,3, Shi-xiong Yang1, Ali-Raza Jahejo1, Ding Zhang1, Wei-fang Hao4, Wen-wei Gao1, Yu-jun Zhao1, Jian-hui Li1, Fang Yan1, Rong-kun Gao1, Yu-hai Bi1,5, Wen-xia Tian1* and Ling-xia Han6* 1College of Veterinary Medicine, Shanxi Agricultural University, Taigu 030801, China 2Weinan Animal Health Supervision Institute, Weinan 714000, China 3Wenxian Forestry Bureau, Wenxian 454850, China 4Taiyuan Center for Disease Control and Prevention, Taiyuan 030024, China 5CAS Key Laboratory of Pathogenic Microbiology and Immunology, Collaborative Innovation Center for Diagnosis and Treatment of Infectious Disease, Institute of Microbiology, Center for Influenza Research and Early-warning (CASCIRE), Chinese Academy of Sciences, Beijing 100101, China 6Unit of Laboratory Animal and Comparative Medicine, State Key Laboratory of Veterinary Biotechnology, Harbin Veterinary Research Institute, Chinese Academy of Agricultural Sciences, Harbin 150001, China

. The major histocompatibility complex (MHC) haplotypes of chickens is the locus with largest effect on MD genetic resistance and susceptibility, but other non-MHC genes are also thought to contribute to MD resistance through cellular interactions, variations in MDV-targeted cells, innate immunity, or cytokine regulation (McPherson and Delany, 2016). For example, the virulent strains of MDV are relatively sensitive to chickens with haplotypes B2, B5, B12, B13 and B15. In the initial phases of MDV pathogenesis, resistant and susceptible birds are similar in degree of viral replication. B-cell lysis and T-cell activation occur within 3-5 dpi and the viral latency begins around 7 to 10 dpi. (Boodhoo et al., 2016). However, subsequent phases of pathogenesis and host responses differ greatly. It has been described that the susceptibility of MD is involved in virus re-activation at about 14 dpi (Boodhoo et al., 2016). MDV re-activation is followed by the transformation and infiltration of CD4 + CD8 -T cells, which results in immunosuppression about 21 dpi (Osterrieder et al., 2006). Erythrocytes are the most abundant cells in circulation and the most significant cells in transportation and exchange of gases. However, the immunomodulatory roles of erythrocytes during viral infections have not yet been thoroughly investigated. A recent research found that the expression of cytokines in salmon erythrocytes was induced by salmon anemia virus (Workenhe et al., 2008). It has been reported that chicken erythrocyte expressed several toll-like receptors (TLRs) and upregulated the expression of some cytokines after TLR ligand treatments (Paolucci et al., 2013).
Data from other groups have explored the pathogenesis of Marek's disease and immune response to MDV by the transcriptional analysis in various organs, and many differentially expressed genes between MDV infected and uninfected chickens have been identified (Sarson et al., 2008;Kano et al., 2009;Heidari et al., 2010). For example, Heidari et al. (2010) have conducted the Gene Chip Chicken Genome Arrays to analyze the pattern of host gene expression in the spleen of MDV-infected chickens. Transcriptome studies have been used for exploring the host response to vaccination in MDV-infected chickens (Kano et al., 2009). Further, the host responses to MDV infection in genetically resistant and susceptible chickens have also been obtained by using transcriptional analysis (Sarson et al., 2008). However, there is little available information about the transcriptome expression of chicken erythrocyte in response to MDV infection. The present study aimed to elucidate the immune responses of chickens to MD by examining the transcriptional profiling in different phases of MDV infection.

Animal infection experiments
Twenty 1-day-old specific-pathogen-free (SPF) chickens (B15 haplotype) and the MD5 strain (very virulent MDV) were kindly provided by Dr. Lingxia Han of Harbin Veterinary Research Institute. The MD5 strain was diluted by sterile phosphate buffered saline (PBS) (Solarbio, Beijing, China). The chickens were randomly divided into two groups (10 chickens/group): Each chicken of the infected group was immunized by an intra-abdominal injection of MD5 with the dose of 500 plaque-forming units (PFU)/500 μL, and each chicken of uninfected control group was injected with the same volume of PBS. All chickens were reared in an animal biosafety level 2 laboratory, and the blood samples were collected from two groups (10 chickens each group) and pooled at 14 and 22 dpi respectively for total RNA extraction. The samples were designated as the control group (14C and 22C), and the treatment group (14T and 22T), where the sample names correspond to 14 dpi and 22 dpi respectively. The management of the experimental animals was in agreement with the welfare guidelines approved by the College of Animal Science and Veterinary Medicine of Shanxi Agricultural University, China (Number 88, 2010).

Erythrocyte isolation and RNA extraction
The erythrocytes were isolated from each blood sample collected at 14 and 22 dpi by a Histopaque-1119 solution (Sigma-Aldrich, Oakville, ON) and RNA extraction from these erythrocyte samples was performed as described previously .

Library construction and RNA sequencing
Library construction and RNA sequencing were performed by Beijing BioMarker Technologies (Beijing, China) according to the institute's protocols. The RNA purity was checked by using the NanoPhotometer spectrophotometer (IMPLEN, CA, USA), and then RNA degradation and contamination were monitored on 1% agarose gel electrophoresis. The RNA integrity was checked using an RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). A total of 1 μg purified RNA per sample was used as input material for the RNA library preparation. Sequencing libraries were generated using NEBNext®Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) according to the manufacturer's instructions, and index codes were added to attribute sequences to each sample. Subsequently, the libraries were sequenced on the Illumina HiSeq 2500. The detailed protocol has been described previously (Wei et al., 2017). All transcriptome sequences were submitted O n l i n e

F i r s t A r t i c l e
in the National Center for Biotechnology Information (NCBI) Sequence Read Archive under the SRA accession number SRP136403.

Quality control
Raw reads were firstly examined by using the inhouse perl scripts. Clean reads were obtained by deleting the low-quality reads and the reads containing adapter or ploy-N from raw data. Then Q30, GC-content and sequence duplication level of the clean data were calculated. All the downstream analyses were based on clean data of high quality.

Differential expression analysis
Prior to differential gene expression analysis, for each sequenced library, the read counts were adjusted by the edgeR program package through one scaling normalized factor. Differential expression analysis of two groups was performed using the DEGseq R package. P value was adjusted using q value for controlling the false discovery rate (Storey, 2003). Gene expressions with q-value < 0.005 and |log2 (foldchange)| ≥ 1 found by DESeq were assigned as significantly differential expression.

Transcriptome assembly, read mapping, and annotation of new genes
The sequencing data of the 4 groups was described in Table I. The data quality assessment has been shown that quality Q30 was higher than 87.31% for each sample and the GC content was 52.48%-54.73%. 77.88%-82.64% clean reads were matched to chicken genome after filtering ( Table I). The quality assessment showed that the transcriptome sequencing data was quite good and reliable. 1583 new genes were discovered that were not found or not annotated before in the reference genome. have not been annotated in chicken reference genome. Among these, 568 genes were not present in the databases so were not annotated, and 1015 new genes were annotated according to the databases including NR, Swiss-Prot, GO, COG, and KEGG, and 568 new genes have not been annotated (Table  II).

Functional annotation and classification of the DEGs between uninfected and infected groups
Total, up-and down-regulated DEGs were detected among these four groups (14T vs 14C, 22T vs 22C). The number of total DEGs (1880) was the highest between the treatment and the control group on 22 dpi, of which 1578 DEGs were up-regulated and 302 were downregulated ( Fig. 1). In the GO terms, three main functional classifications were determined from the obtained DEGs among the four groups: cellular component, biological process and molecular function (Fig. 2). KEGG pathway analysis indicated that the DEGs at each comparison could be grouped into 162 (14T vs 14C), 176 (22T vs 22C), pathways (Supplemental file 1).

Differentially expressed immune-related genes caused by MDV infection
Aiming to identify the DEGs involved in immune pathways and explore the specific changes in the expression of genes caused by MDV infection, all of the DEGs were classified into different functional categories by using GO enrichment analysis (P < 0.05). In the 14T/14C group (Supplemental file 2), genes related to defense response (GO: 0006952), blood vessel development (GO: 0001568), cell morphogenesis involved in neuron differentiation (GO: 0048667), response to oxygen-containing compound (GO: 1901700), blood vessel morphogenesis (GO: 0048514), regulation of mitogen-activated protein kinases (MAPK) cascade (GO: 0043408), activation of immune response (GO: 0002253), and cytokine production (GO: 0001816) were significantly enriched. In the data obtained at 22 dpi (Supplemental file 3), many genes were significantly enriched in the GO terms: the blood vessel development (GO: 0001568), response to oxidative stress (GO: 0006979), regulation of cytokine production (GO: 0001817), cytokine production (GO: 0001816), immune response (GO: 0006955), regulation of T cell mediated immunity (GO: 0002709), positive regulation of immune response (GO: 0050778), inflammatory response (GO: 0006954).
We used common expression patterns to further analyze the DEGs in the groups of 14T/14C, 22T/22C. Based on this method, all of the DEGs were placed into 56 groups, and then three immune-related groups were selected to analyze the immune response of erythrocyte to MDV infection (Fig. 3, Table III  of the three groups classified as involved in these immune-related GO terms: regulation of B cell cytokine production, NK T cell differentiation, interleukin-4 production, positive regulation of I-kappaB kinase/NF-kappaB signaling, positive regulation of interleukin-13 and interleukin-10 biosynthetic process and some other process. The genes involved in these immune-related GO terms showed increased expression in 14T and 22T as compared to 14C and 22C.

The KEGG immune-related pathway mapping for DEGs
In the 14T/14C comparison, there are nine significantly enriched pathways, including endocytosis, ECM-receptor interaction, focal adhesion, herpes simplex infection, cytokine-cytokine receptor interaction, phagosome, cell adhesion molecules, influenza A, ribosome, most of the enriched pathways were related to immune function. In the 22T/22C comparison, pathways related to the cytokinecytokine receptor interaction, toll-like receptor signaling pathway, phagosome, calcium signaling pathway, MAPK signaling pathway, focal adhesion, JAK-STAT signaling pathway were significantly enriched.  More importantly, the multiple occurrences of Cytokine-cytokine receptor interaction, Toll-like receptor signaling pathway, MAPK signaling pathway and many other immune-related pathways (enriched in the 14T/14C and 22T/22C) were particularly striking. These pathways may be a good starting point to characterize the difference in immune reactions during MDV infection. We focused on the toll-like receptor signaling pathway, cytokine-cytokine receptor interaction, MAPK signaling pathway and Janus kinase/signal transducers and activators of transcription (JAK-STAT) signaling pathway (Figs. 4, 5, 6 and 7), and these pathways were enhanced for all groups. Of these, three main gene families (toll-like receptors, gallinacins, interleukins receptors) were upregulated in the 22T group O n l i n e

DISCUSSION
MDV, the etiological agent of MD, is a highly cellassociated oncogenic α-herpesvirus that causes the infiltration of lymphoid cells into various organs. To provide insights into the immune responses of chicken erythrocytes to MDV infection, we investigated the dynamic changes of immune molecules in chicken erythrocytes by transcriptome analysis.
Cytokines play an important role in the regulation of immune responses and inflammation, and cytokine receptors act as transmembrane proteins to transmit a signal into the cell upon ligand binding. Many genome-wide association researches have also demonstrated that polymorphisms and mutations of cytokine receptors lead to some autoimmune disorders (O'shea et al., 2008). It has been shown that high level of type I IFNs transcripts and moderate levels of interleukin (IL)-8 transcripts were expressed in activated chicken erythrocytes, and these cytokine transcripts were up-regulated by TLR3 and TLR21 ligands poly I:C and CpG ODN (Paolucci et al., 2013). In this study, the response of cytokines and cytokine receptors family members to MDV infection has been explored (Fig. 4), and we found that C-X-C motif chemokine (CXCL)-12, 13 and C-C motif chemokine (CCL)-19, 21 and interleukin 2 receptor beta (IL-2RB) were identified and found to be highly expressed in 14 and 22 d after MDV infection, but their expression levels in 22 dpi were lower compared with that of 14 dpi. Similarly, there were many components of cytokinecytokine receptor interaction that were upregulated in 14 and 22 dpi. These data suggested that the expressions of cytokines and their receptors in chicken erythrocytes were induced by MDV infection and suppressed in the phase of immunosuppression. Furthermore, vascular endothelial growth factor (VEGF)-C and VEGF-D act an important role in lymphangiogenesis through the binding of FMSlike tyrosine kinase 4 (FLT4), which promotes lymphatic vessel development and lymphatic endothelial cell (LEC) proliferation (Karkkainen et al., 2004;Mäkinen et al., 2001). In this present study, the transcripts of VEGF-C, VEGF-D and their receptor FLT4 were induced after the infection of MDV, indicating that chicken erythrocytes may participate in lymphangiogenesis by promoting these gene expressions during the development of the MDV infection. Overall, the cytokine-cytokine receptor interaction in chicken erythrocytes may play a significant role in the host immune response to MDV infection.

O n l i n e F i r s t A r t i c l e
Transcriptome Analysis Reveals the Immune Response of Chicken Erythrocytes 9 Fig. 6. Schematic diagrams of the JAK-STAT signaling pathway.
TLRs are a family of innate immune-recognition receptors that play important roles in recognizing molecular patterns associated with pathogen invasion, and they are expressed in various animals from protostomes to vertebrates (Jiménez-Dalmaroni et al., 2016). Pathogen associated molecular patterns (PAMPs) can be recognized by TLRs that signal to the host in the presence of infection via Toll/IL-1 receptor (TIR) domain-containing adaptors, such as MyD88 (Arleevskaya et al., 2019). TLR3, one important member of TLRs family, plays an important role in the antiviral innate immune response (Shah et al., 2016).
It has been reported that activation of TLR3 inhibited MDV infection in chicken embryo fibroblast cells (Hu et al., 2016). Similarly, our data showed that the expression of TLR3 was identified in chicken erythrocytes and upregulated in 14 and 22 dpi, suggesting TLR3 may recognize the doublestranded RNA (dsRNA) produced in the development of MDV infection, and then activate relate signaling pathways to limit the further spread. In addition, some other important TLR gene family members including TLR2, 4, 7, 8 were found to be highly expressed only in 22 dpi (Fig. 5), which indicated that these TLRs may focus on participating in

O n l i n e F i r s t A r t i c l e
the reactivation stage of MDV infection (Hu et al., 2016). The common MyD88-dependent pathway and MyD88-independent signaling pathway can be used by TLRs to activate IRF-3 and NF-κB, and then induce the production of type I interferons and pro-inflammatory cytokines (Ouyang et al., 2018). In this study, MyD88 was highly expressed after MDV infection in both 14 and 22 dpi, which was consistent with the results of TLR3 expression. All the above results supported the possible defensive mechanism of chicken erythrocytes that the dsRNA produced in the development of MDV infection could be recognized by TLRs, and then the common MyD88-dependent pathway and MyD88-independent signaling pathway were used by TLRs and modulate host immune response when MDV invade into the blood. As one of a handful of pleiotropic cascades, the JAK/STAT pathway was used to transduce a multitude of signals for a wide array of cytokines and growth factors in various animals (Morris et al., 2018). Some important cellular events including cell proliferation, differentiation and apoptosis can be stimulated by JAK activation, which is critical to haematopoiesis and immune development. STATs are latent transcription factors that reside in the cytoplasm until activated. JAK/STAT signaling pathway was significantly enriched in all of three DEGs groups, suggesting an important role of the JAK/STAT signaling pathways in the immune response of erythrocytes to MDV infection. In this present study, the results showed that the JAKs including JAK3 and TYK2 were highly expressed after the infection of MDV in both of 14 and 22 dpi, which is consistent with the expression pattern of STATs including STAT1, STAT4 and STAT6 (Fig. 6). These results suggested that JAK/STAT cascade was activated by inducing the JAKs and STATs after MDV infection and it provides a direct mechanism to transform an extracellular signal into a transcriptional response . In addition to JAK/STAT pathway effectors, suppressors of cytokine signaling (SOCS) family of proteins have been recognized as crucial negative regulators of JAK/STAT pathway and they also suppress cytokine signal transduction (Ghafouri-Fard et al., 2018), and a negative feedback loop of the JAK/STAT circuitry can be induced by SOCS. The transcription and protein of the SOCS can be stimulated by activated STATs, and then the JAK/STAT pathway is turned off by SOCS protein binding phosphorylated JAKs and the receptors. Our data showed that SOCS had a higher expression after MDV infection in 14 and 22 dpi, and the expression level of SOCS in 22 dpi is higher than that of 14 dpi. These results support the hypothesis that MDV infection could activated JAKs by various cytokines such as IL2/3/6 and erythropoietin. Activated JAKs subsequently phosphorylated and targets many sites, including STATs, tyrosine-protein phosphatase non-receptor type 11 (PTPN11) and phosphoinositide-3kinase (PI3K). The SOCS family of proteins can be activated by the high expression of STATs, and it may stimulate the negative feedback loop of the JAK/STAT circuitry. Mitogenactivated protein kinases (MAPKs) are evolutionarily conserved kinase modules that mediate intracellular signaling associated with fundamental cellular processes including growth, differentiation and apoptosis (Dhillon et al., 2007). MAPK signaling pathways are comprised of a three-tier kinase module: MAPK kinase kinase (MAP3K), MAPK kinase (MAP2K) and MAPK (Kim and Choi, 2010). In our experiment, many DEGs including tumor necrosis factor receptor (TNFR1), interleukin 1 receptor type II (IL1R2), MAP3K1, MAP3K3, MAP3K4, p38 MAP kinase (p38) and Ras are in the MAPK signaling pathway. Of these DEGs, members of the Ras family of proteins including K-Ras, H-Ras and M-Ras, were upregulated by the MDV infection in 22 dpi, which indicated that Ras plays a key role in transmission of extracellular signals into erythrocytes in the development of MDV infection (Malumbres and Barbacid, 2003). The components of p38 signaling pathway, p38, TNFR1 and IL1R2 were highly expressed after the MDV infection, which suggested that p38 could be activated through the proinflammatory cytokine receptors TNFR1 and IL1R2 (Kim and Choi, 2010). In addition to the DEG MAP3K1 of 14 dpi, the expression of MAP3K3 and MAP3K4 were increased in 22 dpi, indicating that the JNK and p38 signaling pathways were activated after the MDV infection.

CONCLUSION
The immune responses of chicken erythrocyte to MDV infection was well explored by transcriptome analysis. Many DEGs showed significant differential expressions in the Cytokine-cytokine receptor interaction, Toll-like receptor signaling pathway, JAK/STAT signaling pathway and MAPK signaling pathway, which are responsible for immune responses to MDV infection. These results improve our understanding of MDV pathogenesis and complement our knowledge about the role of erythrocytes immune response on MDV infection mechanism in chicken.

F i r s t A r t i c l e
Shanxi Province (20130311027-3), and "131" Leading Talent Project for Colleges and Universities of Shanxi Province.

Supplementary material
There is supplementary material associated with this article. Access the material online at: https://dx.doi. org/10.17582/journal.pjz/20191012061013

Statement of conflicts of interest
The authors have declared no conflict of interests.