Global gene expression analysis of peripheral blood mononuclear cells in rhesus monkey infants with CA16 infection-induced HFMD

Virus Research 214 (2016) 1–10

Contents lists available at ScienceDirect

Virus Research journal homepage: www.elsevier.com/locate/virusres

Global gene expression analysis of peripheral blood mononuclear cells in rhesus monkey infants with CA16 infection-induced HFMD Jie Song, Yajie Hu, Yunguang Hu, Jingjing Wang, Xiaolong Zhang, Lichun Wang, Lei Guo, Yancui Wang, Ruotong Ning, Yun Liao, Ying Zhang, Huiwen Zheng, Haijing Shi, Zhanlong He, Qihan Li ∗ , Longding Liu ∗ Institute of Medical Biology, Chinese Academy of Medical Science and Peking Union Medical College, Kunming, 650118, China

a r t i c l e

i n f o

Article history: Received 13 November 2015 Received in revised form 5 January 2016 Accepted 5 January 2016 Available online 8 January 2016 Keywords: Coxsackievirus A16 (CA16) Global gene expression analysis Peripheral blood mononuclear cells Rhesus monkey infants

a b s t r a c t Coxsackievirus A16 (CA16) is a dominant pathogen that results in hand, foot, and mouth disease and causes outbreaks worldwide, particularly in the Asia-Pacific region. However, the underlying molecular mechanisms remain unclear. Our previous study has demonstrated that the basic CA16 pathogenic process was successfully mimicked in rhesus monkey infant. The present study focused on the global gene expression changes in peripheral blood mononuclear cells of rhesus monkey infants with hand, foot, and mouth disease induced by CA16 infection at different time points. Genome-wide expression analysis was performed with Agilent whole-genome microarrays and established bioinformatics tools. Nine hundred and forty-eight significant differentially expressed genes that were associated with 5 gene ontology categories, including cell communication, cell cycle, immune system process, regulation of transcription and metabolic process were identified. Subsequently, the mapping of genes related to the immune system process by PANTHER pathway analysis revealed the predominance of inflammation mediated by chemokine and cytokine signaling pathways and the interleukin signaling pathway. Ultimately, co-expressed genes and their networks were analyzed. The results revealed the gene expression profile of the immune system in response to CA16 in rhesus monkey infants and suggested that such an immune response was generated as a result of the positive mobilization of the immune system. This initial microarray study will provide insights into the molecular mechanism of CA16 infection and will facilitate the identification of biomarkers for the evaluation of vaccines against this virus. © 2016 Elsevier B.V. All rights reserved.

1. Introduction Hand, foot, and mouth disease (HFMD) is a global and serious epidemic disease that frequently affects infants and children younger than 5 years of age, and it is a continuing threat to public health (Wong et al., 2010; Centers for Disease Control and Prevention, 1998). It is predominantly caused by viruses that belong to the enterovirus genus, including polioviruses (3 serotypes), coxsackievirus A (23 serotypes), coxsackievirus B (6 serotypes), enteroviruses (29 serotypes) and echoviruses (28 serotypes) (Solomon et al., 2010). Of these, CA16 and enterovirus 71 (EV71) are the principal etiological pathogens of HFMD (Ang et al., 2009). HFMD is generally a mild, self-limiting infection characterized by a typical maculopapular rash and blisters on the hands

∗ Corresponding authors. Fax: +86 871 68334483. E-mail addresses: [email protected] (Q. Li), [email protected] (L. Liu). http://dx.doi.org/10.1016/j.virusres.2016.01.002 0168-1702/© 2016 Elsevier B.V. All rights reserved.

and feet and in the mouth. However, numerous severe HFMD cases with central nervous system complications, such as encephalitis, meningitis and pulmonary edema, which can cause serious injury or even death in young children, have been reported (Solomon et al., 2010; Huang et al., 2012; Zhang et al., 2011a,b). In recent years, CA16 and EV71 have often been observed to circulate alternatively or together (Pan et al., 2009; Zhang et al., 2011a,b; Iwai et al., 2009). Large-scale HFMD outbreaks causing significant morbidity and mortality have occurred in mainland China. According to the Chinese Center for Disease Control, 4,607,238 cases were reported from 2013 to 2014, including 753 deaths that were primarily due to severe HFMD with complications. Many current studies have focused on EV71 because of its higher lethality rate, but relatively fewer studies have focused on CA16. Previous research performed by our laboratory has resulted in the successful development of an inactivated EV71 vaccine, for which a phase 3 clinical trial has been completed (Li et al., 2014a,b; Liu et al., 2013). With future availability of the EV71 vaccine, morbidity caused by CA16 and

2

J. Song et al. / Virus Research 214 (2016) 1–10

other enteroviruses may increase. Hence, efforts must be made to further understand CA16-induced pathogenesis and the associated immune response by profiling gene expression patterns in an infection model. In our previous studies, the basic CA16 pathogenic process was successfully mimicked in rhesus monkey infants. These monkeys exhibited the typical clinical manifestations, such as vesicles on the mucosa of the palate, tongue and limbs, fever, and the histopathological manifestations, such as a variety of viral loads, antigen expression and perivascular infiltration in the central nervous system (CNS) (Wang et al., 2014). In this study, we investigated the possible molecular mechanisms of CA16 infection by performing global gene expression analysis of peripheral blood mononuclear cells (PBMCs) in a CA16-infected neonatal rhesus monkey model.

samples collecting, all of the monkey infants in the experimental and control groups were euthanized via an overdose of anesthesia. 2.3. Infection of neonatal rhesus monkeys with CA16

2. Materials and methods

A total of 6 rhesus monkey infants were randomly divided into the following two groups: a control group (n = 3) and a CA16-infected group (n = 3). After the initial viral infection was transmitted via nasal atomization (4.5 lgCCID50/ml), all clinical symptoms were observed and recorded daily. Blood samples were collected on days 0 and 7 post-infection (7-1st dpi). The subsequent viral challenges were performed at 1-month intervals and were followed by the same sample collections and clinical observations on days 3, 7 and 14 post-infection (3-2nd dpi, 7-2nd dpi, 14-2nd dpi). The detailed process of CA16 infection is shown in Fig. S1. Supplementary material related to this article found, in the online version, at http://dx.doi.org/10.1016/j.virusres.2016.01.002.

2.1. Virus and cells

2.4. Microarray screening of RNA samples

The CA16 virus G20 stain (sub-genotype B) was isolated from a patient during the 2010HFMD outbreak in Guangxi, China (GenBank: JN590244.1). The virus was grown in Vero cells (ATCC, Manassas, VA, USA) and harvested for freezing at −20 ◦ C. The Vero cells were maintained in Dulbecco’s modified Eagle’s medium (DMEM) (HyClone, Logan, UT, USA) supplemented with 10% fetal bovine serum (FBS) (Gibco, Grand Island, NY, USA) at 37 ◦ C in a CO2 incubator.

To exclude individual differences, three experimental and control samples were mixed for each time point. A Whole Rhesus Monkey Genome Microarray (G2519F-026806, GPL16026, Agilent, CA, USA) was chosen to screen for gene expression differences in monkey PBMCs. Gene Chip microarray experiments were conducted at the National Engineering Center for Biochip in Shanghai, China, according to the procedures described in the Agilent technical manual. Briefly, mRNA purified from total RNA after the removal of rRNA was amplified and transcribed into fluorescent cRNA by Low Input Quick Amp Labeling, according to the manufacturer’s protocol (Agilent). Labeled cRNA was purified using an RNeasy Mini Kit (QIAGEN, GmBH, Germany). Each slide was hybridized with 1.65 ml Cy3-labeled cRNA using a Gene Expression Hybridization Kit (Agilent) in a Hybridization Oven (Agilent). After 17 h of hybridization, the slides were washed in staining dishes (Thermo Scientific, MA, USA) with a Gene Expression Wash Buffer Kit (Agilent). The slides were then scanned using an Agilent Microarray Scanner (Agilent), and raw data were obtained using Feature Extraction Software 10.7 (Agilent) and normalized using the quantile algorithm with Gene Spring 11.0 (Agilent). Systemic bioinformatic analyses of the microarray data were performed with an SBC Analysis System (provided by Shanghai Biotechnology Corporation). Briefly, the normalization value was set to 1, and differentially expressed genes with fold changes of ≥2 or ≤0.5 were analyzed. Significantly enriched gene ontology (GO) terms and pathways were determined from these differentially expressed genes using the PANTHER Classification System, and then Gene MANIA was used to build a gene network according to the relationships among the genes, proteins and compounds in the database (Mi et al., 2013; Mi and Thomas, 2009). The raw microarray data were submitted to the Gene Expression Omnibus database and are available under the accession number GSE66757.

2.2. Ethics statement All animals were handled in accordance with the principles expressed in the “Guide for the Care and Use of Laboratory Animals” by the National Research Council of the National Academies (National Academy of Science, 2011) and “The Guidance to Experimental Animal Welfare and Ethical Treatment” by The Ministry of Science and Technology of the People’s Republic of China (The guidance to experimental animal welfare and ethical treatment, 2006) The Yunnan Provincial Experimental Animal Management Association (approval number: SCXK (Dian) 2011–0005) and the Experimental Animal Ethics Committee of the Institute approved the animal research (approval number: YISHENGLUNZI [2013] 2). All of infant rhesus macaques in this study wereseparated from their mothers. And they (2–3 months of age) were kept in single stainless steel cages (L: 560 mm, H: 767 mm, W: 700 mm) in a large room (BSL-2 conditions) with sufficient fresh air and natural light at approximately 25 ◦ C. Free access to food and water supplemented with the appropriate treats and vitamins were provided. The animals were given access to environmental enrichment (such as approved toys) to promote psychological well being. Clinical symptoms, including body temperature changes and vesicles in the oral mucosa or limbs, were monitored and recorded daily (Wang et al., 2014). Blood samples were collected under appropriate anesthesia to alleviate pain and minimize suffering, in compliance with the guidelines of the Institute of Medical Biology (IMB), Chinese Academy of Medicine Science (CAMS). Following recovery from anesthesia after blood samples collection, all of the animals were returned to the cages. They were all placed fully under the care of veterinarians at the IMB, CAMS. The housing conditions, experimental procedures and animal welfare were in accordance with the local laws and guidelines on the use of laboratory non-human primates and were in compliance with the recommendations of the Weatherall report. The absence of anti-CA16 antibodies was verified in all macaques prior to experimentation. No unexpected deaths occurred during the experiment. After the last time of blood

2.5. RNA extraction and qRT-PCR amplification To confirm the accuracy of the microarray data, we randomly selected various differentially expressed genes and determined their mRNA levels by qRT-PCR, using ␤-actin mRNA as an internal standard. Total RNA was extracted from PMBCs using a GeneJET RNA Purification Kit (Thermo Scientific, MA,USA) according to the manufacturer’s recommendations, and the total RNA concentrations of the samples were calculated by measuring optical density at 260 nm with a Nanodrop 2000 spectrophotometer (Thermo Scientific, USA). The extracted RNA was temporarily frozen in 75% ethanol until further testing. One microgram of total RNA from each sample was processed directly into cDNA with a First Strand

J. Song et al. / Virus Research 214 (2016) 1–10

3

Table 1 Statistics of differentially expressed genes. Combination

Up-regulated genes Down-regulated genes Total genes

7-1st dpi vs. 0 dpi 3-2nd dpi vs. 0 dpi 7-2nd dpi vs. 0 dpi 14-2nd dpi vs. 0 dpi

1,950 2,133 2,495 2,320

1,903 2,191 2,347 2,497

3,853 4,324 4,842 4,817

qRT-PCR validation of differentially expressed genes identified by microarray.

Fig. 1. Venn diagram showing a comparison of the differentially expressed genes at different time points. The overlapping differentially expressed genes are shown between 7-1st dpi vs. 0 dpi, 3-2nd dpi vs. 0 dpi, 7-2nd dpi vs. 0 dpi, and 14-2nd vs. 0 dpi. The list of 948 differentially expressed genes in the samples exhibits significant overlap.

cDNA Synthesis Kit (Thermo Scientific, USA) in accordance with the manufacturer’s instructions. Subsequently, the cDNA was used to perform SYBR Green qPCR with a7500 Fast Real-time RT-PCR system (Applied Biosystems, Foster City, CA, USA) according to the manufacturer’s protocol. The primer sequences, amplification lengths and Tm values used in this study are listed in Table S1. The primers were designed by Oligo7 software. The amplification reactions were all carried out for 40 cycles as follows: denaturation at 95 ◦ C for 5 s and annealing at Tm for 34 s. Supplementary material related to this article found, in the online version, at http://dx.doi.org/10.1016/j.virusres.2016.01.002. 2.6. Statistical analysis The data are showed as the mean and standard deviation. Individual analyses were performed in triplicate. GraphPad Prism software (San Diego, CA, USA) was used for statistical analyses. 3. Results 3.1. Analysis of differentially expressed genes in CA16-infected rhesus monkey infant PBMCs Microarray technology allows for the simultaneous examination of thousands of genes and for a better understanding of gene expression patterns, which can reveal a range of biological characteristics and the potential molecular basis of clinical diseases (Zhang et al., 2014). Using Agilent 4 × 44 K Rhesus Monkey WholeGenome Microarrays, we identified 948 differentially expressed genes with significant changes at all 4 time points tested during CA16 infection. However, pair-wise analyses between the CA16infected samples, including the 7-1st dpi, 3-2nd dpi, 7-2nd dpi, and 14-2nd dpi samples, and the CA16 un-infected sample revealed the altered expression of 3,853, 4,324, 4,842 and 4,817 genes, respectively (Fig. 1). As shown in Table 1, the overall gene expression patterns showed remarkable up- or down-regulation at different time points following CA16 infection to varying degrees compared with 0 dpi. Additionally, principal component analysis (PCA) and hierarchical clustering analysis of the sample-level data were both performed to determine the similarity among the different samples under each condition. PCA demonstrated that the 7-1st dpi, 3-2nd dpi, 7-2nd dpi, and 14-2nd dpi samples were highly correlated; further, they were segregated from the 0 dpi sample (Fig. 2A).

Subsequently, 948 differentially expressed genes were analyzed by hierarchical clustering among the 5 samples (Fig. 2B), revealing the marked up- or down-regulation of these distinctive genes. To confirm the results obtained from microarray analysis, 16 of the differentially expressed genes, including C3, CD4, IL1A, MAPK12, MMP1, NRGN, PID1, CXCL11, CXCR6, IFNG, IL10, IRF7, JAK2, NTS, NEUROD4, and STAT1, were chosen for qRT-PCR with a randomly selected subset of 5 PMBC samples, including 0 dpi, 71st dpi, 3-2nd dpi, 7-2nd dpi, and 14-2nd dpi samples. We found that the CA16 infection groups had notably increased C3, CD4, IL1A, MAPK12, MMP1, NRGN and PID1 expression levels and markedly decreased CXCL11, CXCR6, IFNG, IL10, IRF7, JAK2, NTS, NEUROD4, and STAT1 expression levels (Fig. 3), consistent with the microarray data. 3.2. Functional categories and significant ontology analysis of differentially expressed genes Functional category and significant ontology analyses were performed to help to clearly and comprehensively understand the involvement of the differentially expressed genes in a wide variety of biological processes related to the occurrence and development of the disease. A total of 948 genes were differentially expressed, as determined by hierarchical clustering analysis. Next, we investigated whether these distinct genes have particular functions via functional classifications based on GO terms from GO Slim using the PANTHER gene list analysis tool. These 948 genes were primarily divided into 5 groups according to their GO biological process descriptions, including cell communication, cell cycle, immune system process, regulation of transcription and metabolic process, which comprised 135, 53, 115, 59 and 312 genes, respectively (Fig. 4). The log 2 of fold change related to these genes and numbers of up- regulated or down- regulated genes in different time point of the five groups are shown in Fig. 5. There are 60, 55, 53, 58 of upregulated genes and 54, 59, 61, 56 of down- regulated genes in 71st dpi, 3-2nd dpi, 7-2nd dpi, and 14-2nd dpi samples, respectively, in immune system process. And the details on the gene expression changes in the enriched functional categories in the five groups are presented in Fig. S2. All of the five functional categories were closely associated with the progressive pathogenesis of CA16 infection. Supplementary material related to this article found, in the online version, at http://dx.doi.org/10.1016/j.virusres.2016.01.002. 3.3. Co-expression network analysis of the top 10 differentially expressed genes in 7-1st dpi rhesus monkey infants Cellular processes are complex; thus, to fully understand the functions of these processes, biological networks depicting interactions among genes are particularly important. Many genes with different expression levels that are closely correlated may display common regulatory mechanisms or participate in similar biological processes, which ultimately induce different patterns of co-expression. We created gene co-expression networks describing the inter-relationships among the genes using GeneMANIA to assess the data obtained for the top 10 differentially expressed genes, which were chosen from the five ontology functional cat-

4

J. Song et al. / Virus Research 214 (2016) 1–10

Fig. 2. PCA, heat map and supervised hierarchical clustering analysis. (A) PCA demonstrated that the 7-1st dpi, 3-2nd dpi, 7-2nd dpi, and 14-2nd dpi samples were highly correlated and that these samples were segregated from the 0 dpi sample. (B) A heat map and supervised hierarchical clustering analysis revealed that 948 genes were differentially expressed at the 5 time points after CA16 infection. Each row represents a gene, and the samples are depicted in the columns. Red indicates genes that were expressed at higher levels, and blue denotes genes that were expressed at lower levels compared with the expression levels of 0 dpi. The color bars represent the log 2 of fold change.

Fig. 3. Confirmation of the gene expression changes by qRT-PCR. A total of 16 individual genes were randomly selected and analyzed by qRT-PCR. The y-axis indicates the relative quantities of the specific mRNAs in the 7-1st dpi vs. 0 dpi (A), 3-2nd dpi vs. 0 dpi (B), 7-2nd dpi vs. 0 dpi (C) and 14-2nd vs. 0 dpi samples (D). The results are normalized by the endogenous ␤-actin expression level. Individual analyses were performed in triplicate. The error bars indicate the mean ± SD.

J. Song et al. / Virus Research 214 (2016) 1–10

5

pathways in all and the proportion of these pathways was shown in Fig. 7. Of these pathways, inflammation mediated by chemokine and cytokine signaling pathways and the interleukin signaling pathway comprised a relatively large proportion, accounting for 21.1% and 10.5%, respectively. The former pathways mainly contain cyclooxygenase (12.5%), chemokines (25.0%) and chemokine receptor (37.5%), whereas the latter primarily contains the interleukin receptor beta subunit (22.2%), signal transducers and activators of transcription (22.2%) and interleukins (22.2%), which may imply that they exert strong impacts on the potential modulation of the immune response. 4. Discussion Fig. 4. CA16 infection induced gene expression changes in infected rhesus infant PBMCs. Gene ontology analysis of the 948 overlapping differentially expressed genes in PBMCs from the 0 dpi, 7-1st dpi, 3-2nd dpi, 7-2nd dpi, 14-2nd dpi samples were divided into 5 groups, including cell communication, cell cycle, immune system process, regulation of transcription and metabolic process. The values are presented as the gene numbers.

egories described above in the gene expression profiles of the 7-1st dpi rhesus monkey infants (Fig. 6). Moreover, the top 10 differentially expressed genes are listed in Table 2. As shown in Fig. 6, the top 5 key genes were modulated by the top 10 differentially expressed genes in these gene co-expression networks. 3.4. Enriched PANTHER pathways of differentially expressed genes involved in the immune response To further explore the underlying regulatory mechanisms associated with the immune response after CA16 infection, a PANTHER pathway analysis pie chart was created for the 115 differentially expressed immune system process genes that were identified in ontology analysis (Fig. 7). These pathways included 27 different

A high mortality rate of HFMD caused by EV71 has been reported in different countries in the West Pacific region, while CA16 infection has only received limited attention because of its mild symptoms (Wang and Liu, 2014). Remarkably, many reports have shown that CA16 is responsible for nearly 50% of all HFMD cases in mainland China and that it can cause severe complications and even death (Yen et al., 2009; Wright et al., 1963). With increasing public concern about CA16 infections, the profiling of host gene expression patterns related to viral infection could be valuable for the understanding of the pathogenesis of CA16. In our previous work, we have revealed a close correlation between the mechanism of EV71 pathogenicity and an abnormal immune response (Zhang et al., 2014). Further, the findings of this study obtained using a neonatal rhesus monkey model indicated that CA16 infection resulted in the differential expression of 948 host genes related to cell communication, the cell cycle, immune system processes, regulation of transcription and metabolic processes. Recently, accumulating evidence declared that several viruses utilized multifaceted strategies to evade multiple arms of the innate

Fig. 5. The fold change of different genes in different time point of the five groups. A, B, C, D, and E represent log 2 of fold change and numbers of up- regulated or down- regulated genes of different genes in different time point compared to 0 dpi related to cell communication, cell cycle, immune system process, regulation of transcription and metabolic process, respectively, as shown.

6

J. Song et al. / Virus Research 214 (2016) 1–10

J. Song et al. / Virus Research 214 (2016) 1–10

7

Table 2 Top 10 differentially expressed genes. Gene ontologies

Gene symbols of top 10 up-regulated genes

Fold changes (7-1st dpi/0 dpi

Gene symbols of top 10 down-regulated genes

Fold changes (7-1st dpi/0 dpi

Cell communication

FPR1 MYOM2 CXCL3 IL1A FN1 CSF3 GDF15 PKIB CXCL1 CXCL2

85.05858795 41.25082878 36.69877495 29.81475342 27.06999282 24.84493737 23.9397304 23.11380707 22.72726224 18.21926684

MOS MYLK4 GABRA2 NTS ELF5 CXCL9 CXCL11 ARHGAP19 CXCL10 IFNG

0.048253995 0.043603804 0.035197286 0.034572705 0.034501039 0.026708023 0.019628665 0.015975662 0.015437019 0.013910663

Cell cycle

IL1A BAX EGR1 FOSB RBM47 MYO7B FOSL2 OTUD1 RALGDS ACTA2

29.81475342 12.60590246 12.03572631 11.23358728 10.72646503 6.030482763 5.363203516 5.1050196 4.202572768 4.073692829

CDC42BPA DDX60 WISP3 TMF1 ANKS1B EXO1 MS4A5 MOS MYLK4 ELF5

0.200792808 0.094629598 0.087713578 0.083967969 0.070745696 0.062429828 0.056704412 0.048253995 0.043603804 0.034501039

Immune system process

PTGS2 FPR1 CXCL3 IL1A CSF3 CXCL1 OLR1 CCL2 CXCL2 C3

177.3102957 85.05858795 36.69877495 29.81475342 24.84493737 22.72726224 21.54288706 20.94806382 18.21926684 17.07011651

COLQ GBP2 MOS FASLG ELF5 CXCL9 CXCL11 ARHGAP19 CXCL10 IFNG

0.064498374 0.06176278 0.048253995 0.04303262 0.034501039 0.026708023 0.019628665 0.015975662 0.015437019 0.013910663

Regulation of transcription

EGR1 FOSB KLF4 IRF7 NLRP3 GTF2B FOSL2 SDPR SOX14 SREBF1

12.03573 11.23359 10.32259 7.004189 6.498986 5.372919 5.363204 5.249095 3.645066 3.215241

TCF7L1 NLRP4 WDR72 PHF16 PDX1 CDX1 RAX NEUROD1 ZSCAN2 ELF5

0.109449 0.108157 0.081591 0.070242 0.064208 0.060477 0.053982 0.051253 0.05096 0.034501

Metabolic process

RPS4Y2 MMP1 PTGS2 EIF1AY APOC2 PYGL MYOM2 ELA2B SERPING1 SERPINB2

1889.7292 306.02846 177.3103 111.99197 53.442404 43.512756 41.250829 33.535226 25.205321 24.677032

MYLK4 PSAT1 ELF5 PFKFB3 ADCY2 GUCY1A3 TRIM9 GBP1 HAO2 ARHGAP19

0.0436038 0.0401777 0.034501 0.0337928 0.0300327 0.0250808 0.0242034 0.0211374 0.0199511 0.0159757

immune response and promoted their replication and spreading, including hiding viral from cytosolic pattern recognition receptors (PRRs), blocking the anti-viral pathway, producing microRNA and protease to destroy the function of host cell and so on (Wang et al., 2013; van Gent et al., 2015; Xu et al., 2014). In our study, IFN-␥, which is involved in the first line of defense against viral infections, possesses potent antiviral, anti-proliferative, and immunomodulatory activities (Saha et al., 2010; Schoenborn and Wilson, 2007) was the most highly down-regulated among the top 10 differentially expressed genes and was correlated with both cell communication and immune system processes. Moreover, it was reported that the EV71 3C protein, which disrupts the formation of a func-

tional retinoid acid-inducible gene I (RIG-I) complex and inhibits the function of IRF7, suppressed type I IFN responses and finally leaded to exacerbate the virus-induced disease. Therefore, the dramatic decrease in IFN-␥ expression was concluded to be due to the CA16 infection, which disrupted the first line of defense in the rhesus monkey infants, thereby leaving host cells susceptible to CA16 invasion. Chemokines are a type of cytokine that can induce the movement of some immune cells or inflammatory cells to target sites of inflammation (Chen and Ichinohe, 2015). We found that the expression of some chemokines were altered, such as CXCL1, CXCL2, CXCL3, CXCL9, CXCL10, and CXCL11, which were all hypothesized to participate in the initiation of the immune

Fig. 6. Co-expression network analysis of top 10 differentially expressed genes in 7-1st dpi rhesus monkey infants. The top 5 key regulatory genes were regulated by the top 10 differentially expressed genes in the 7-1st dpi rhesus monkey infants.

8

J. Song et al. / Virus Research 214 (2016) 1–10

Fig. 7. Enriched PANTHER pathways of differentially expressed genes involved in the immune response. A PANTHER pathway analysis pie chart is shown for the 115 overlapping differentially expressed immune system process genes subjected to ontology analysis after CA16 infection at the0 dpi, 7-1st dpi, 3-2nd dpi, 7-2nd dpi, and 14-2nd dpi time points.

response to CA16 infection. Formylpeptide receptor 1 (FPR1) was the second highest up-regulated gene among the top 10 differentially expressed genes and was classified as an immune system process gene. However, FPR1, a G protein-coupled receptor, is not only an essential player in signal transduction for phagocyte accumulation, but it is also necessary for promoting innate immunity as an immune system process gene (Liu et al., 2012). Many viruses have been reported to modulate the cell cycle to increase their replication efficiency and to cause host cell apoptosis by cell cycle arrest (Rodrigues et al., 2014). Furthermore, it has been discovered that CA16 and EV71 both prevented the cell cycle transition from S phase into the G2/M phase in diverse cell types, which facilitated their own growth (Yu et al., 2015). However, in our study, Bax, which exerts predominantly pro-apoptotic effects during the early phase of apoptosis (Li et al., 2014a,b), was notably increased, further indicating that CA16 infection may mediate the host cell cycle and lead to cell apoptosis. Autophagy, as well as apoptosis, is a formation of cell death which is benefit for maintaining the homeostasis. Also, there is evidence that CA16 could elicit an incomplete autophagy by inhibiting Akt/mTOR signaling and triggering the extrcellular signal-regulated kinase (ERK) signaling to enhance its own replication (Shi et al., 2015). In our results, VPS54, which is a subunits of Golgi-associated retrograde protein (GARP) complex, was remarkable down-regulated after CA16 infection and it was presumed that this change may impair protein retrieval to the trans-Golgi network (TGN), break lysosomal enzyme sorting and endosomal cholesterol traffic, further affect the occurrence of autophagy (Reggiori et al., 2003). Nevertheless, viruses can also take advantage of host cell metabolite production to guarantee their replication and dissemination; on the other hand, they can destroy host cell metabolic machinery. For example, viruses use lipid scaffolds as a physical support to accomplish their replication and assembly in host cells; meanwhile, they use the nucleic acids of host cells for their own replication, resulting in the blockade of RNA (or DNA) synthesis (Rodrigues et al., 2014). RPS4Y2, which encodes a ribosomal pro-

tein, thereby accelerating protein synthesis (Andrés et al., 2008), was the most highly up-regulated metabolic process gene; therefore, we postulate that CA16 infection may cause aberrant protein metabolism. Transcriptional regulation is necessary for gene expression in biological processes involved in virus-host interactions, including all of the above mentioned biological processes. Therefore, transcriptional regulatory factors are also especially important in the modulation of gene expression. For example, IRF7, one of the top 10 up-regulated differentially expressed genes, is a key regulator of the type I interferon signaling pathway in protecting against pathogenic infections (Ning et al., 2011). Therefore, genes associated with transcriptional regulation may play important roles in the modulation of other biological processes. The differentially expressed genes that are involved in the different biological processes described above were considered to have pivotal impacts on the pathogenesis of CA16 infection in this study. Additionally, of the 948 differentially expressed genes, the expression of NRGN and NTS, which have been implicated in learning and memory functions, was markedly elevated and reduced, respectively, after CA16 infection. NRGN was especially enriched in CA1 pyramidal neurons in the hippocampus, as it regulates the calcium-calmodulin (Ca-CaM) signaling pathway by binding to the Ca2+ -free form of CaM and subsequently activates NMDA receptor, which can induce long-term potentiation (LTP), which is related to the formation of spatial learning and hippocampal plasticity (Ohi et al., 2012). While NTS is abundantly expressed in the entorhinal cortex, which is important for learning and memory, NTS may have a therapeutic effect by improving cognitive function (Xiao et al., 2014). Nevertheless, most clinical cases with severe neurologic complications caused by EV71 infection exhibit delayed neurodevelopment or reduced cognitive function (Wang et al., 2012; Weng et al., 2010; Chang et al., 2007). Accordingly, these data indicate that cases of serious HFMD caused by CA16 are accompanied by neurologic complications that may be closely linked to expression changes in NRGN and NTS to some extent.

J. Song et al. / Virus Research 214 (2016) 1–10

What’s more, it has been reported that CA16 infection could be vulnerable to a repeated infection in clinic and the underlying mechanism of this phenomenon was still unknown (Zou et al., 2012). However, in our research, the numbers of up- regulated genes and down- regulated genes related to immune system process after the initial infection are similar to the second infection, which maybe further suggest that CA16 has successfully established repeated infection in neonatal rhesus monkeys. To further explore the role of the immune response in the pathogenesis of CA16 infection, we focused our study on changes in the immune response by employing the PANTHER Classification System. We observed that some inflammatory markers belonged to the chemokine and cytokine signaling pathways, including IL-6, IL8, IL-18, CXCL5, and CCL20, all of which were significantly altered after CA16 infection. IL-6, a multifunctional cytokine, has several diverse functions, and it not only promotes monocyte differentiation and increases neutrophils via various signaling pathways in innate immunity, but it also stimulates antibody synthesis and secretion by B lymphocytes (Okino et al., 2014). IL-8 is a chemokine that is involved in the local inflammatory response, helping to recruit neutrophils to sites of infection (Morales-Garcia et al., 2012). IL-18 is a cytokine that primarily enhances IFN-␥ production (Serti et al., 2014). CXCL5 and CCL20 are the ligands of CCR2 and CCR6, respectively, and these chemokines are crucial for inflammatory cell recruitment (Nouailles et al., 2014). Furthermore, the CCR6/CXCL20 axis is involved in immune surveillance and autoimmunity pathways in the central nervous system (Reboldi et al., 2009). Overall, the cytokines and chemokines described above are highly associated with inflammation mediated bychemokine and cytokine signaling pathways, as well as the interleukin signaling pathway. Furthermore, it was reported that EV71-initated intracellular signaling pathway, including p38 MAPK, JNK, p42/p44 MAPK, NF-␬B and AP-1 pathway tightly linked to inflammation, caused increased cyclooxygenase-2 (COX-2) expression and PGE2 generation in a time- and virus titer-dependent manner and consequently evoked a severe inflammation in brain (Tung et al., 2010). However, the data in our microarrays is consistent with this report that PTGS-2, also known as COX-2, is the top one of up-regulated gene in immune response process. Thereby, the genes above descriptions referred to inflammation were whether intimately related to the serious degree of this disease induced by CA16-infection and this underlying mechanism needed to further investigate. Additionally, some genes in the co-expression network of the immune system, such as STAT1 and TAP1, which act as transcription factors, were down-regulated post-CA-16 infection. Whereas STAT1 acts downstream of IFN-␥ and performs its antiviral function via an intercellular transcriptional pathway (Saha et al., 2010), TAP1 is a transporter that is associated with antigen processing and has been implicated in immune response regulation (Schmitt and Tampe, 2000). Thus, the changes in the expression of these genes that are associated with the immune response may partly explain the abnormal immune response elicited by CA16 infection, which ultimately leads to the pathogenesis of CA16 infection. Furthermore, the different categories in the co-expression network have a disparity in the degree of integration. Nevertheless, the degree of integration in the co-expression network largely depends on that these genes participated in co-expression network referred to numerous different pathways and then involved in the regulation of these pathways and finally generated a complex network. Therefore, IL1A, CXCL1, CXCL2 and CXCL3, as cytokines, which have a close interaction relationship and similar regulation pathways, including cell migration and inflammation pathway and so on, are all involved in the co-expression networks of cell communication and immune responses which had a higher integration and generated a more complex network. However, the others co-expression network generated by the genes, which have a less interaction

9

relationship, are scattered. Thereby, it was presumed that these relatively loose genes have a remote regulation which leaded to a less linking in the co-expression networks of cell cycle, immune system process and regulation of transcription. The enormous information from the analysis of the global gene expression of PBMCs in a CA16-infected rhesus monkey infant model give us multiple cues to explore the deeply mechanisms of pathogenesis resulted from CA16-infection. For example, the differences expression profile of CA16-infection and the published expression profile of EV71-infection should be compare to explore the mechanisms that caused different symptom by EV71 and CA16 infections on host. Moreover, how CA16 infection caused CNS impairment and whether NRGN or NTS was involved in this impairment needed further investigation. In conclusion, we first examined the global gene expression of PBMCs in a CA16-infected rhesus monkey infant model, which not only helped to further the understanding of CA16 infection but also provided new insights that warrant further detailed investigations regarding the patho-mechanisms of CA16 infection. Conflict of interest The authors have declared that no competing interests exist. Acknowledgements This work was supported by National Natural Sciences Foundation of China (81373142), Important National Science & Technology Specific Projects (2014ZX09102042) and Science and Technology Projects of Yunnan Province (2012ZA009). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. References Andrés, O., Kellermann, T., López-Giráldez, F., Rozas, J., Domingo-Roura, X., Bosch, M., 2008. RPS4Y gene family evolution in primates. BMC Evol. Biol. 8, 142, http://dx.doi.org/10.1186/1471-2148-8-142. Ang, L.W., Koh, B.K., Chan, K.P., Chua, L.T., James, L., Goh, K.T., 2009. Epidemiology and control of hand, foot and mouth disease in Singapore, 2001–2007. Ann. Acad. Med. (Singapore) 38, 106–112. Centers for Disease Control and Prevention, 1998. Deaths among children during an outbreak of hand, foot, and mouth disease—Taiwan, Republic of China, April–July 1998. Morb. Mortal Wkly. Rep. 47, 629–632. Chang, L., Huang, L., Gau, S.S., Wu, Y., Hsia, S., Fan, T., et al., 2007. Neurodevelopment and cognition in children after enterovirus 71 infection. N. Engl. J. Med. 356, 1226–1234, http://dx.doi.org/10.1056/NEJMoa065954. Chen, I., Ichinohe, T., 2015. Response of host inflammasomes to viral infection. Trends Microbiol. 23, 55–63, http://dx.doi.org/10.1016/j.tim.2014.09.007. China. The guidance to experimental animal welfare and ethical treatment, 2006. Available http://www.most.gov.cn/fggw/zfwj/zfwj2006/200609/t20060930 54389.htm (accessed 30.08.06.). Huang, H., Weng, K., Shih, S., 2012. Viral and host factors that contribute to pathogenicity of enterovirus 71. Future Microbiol. 7, 467–479, http://dx.doi. org/10.2217/fmb.12.22. Iwai, M., Masaki, A., Hasegawa, S., Obara, M., Horimoto, E., Nakamura, K., et al., 2009. Genetic changes of coxsackievirus A16 and enterovirus 71 isolated from hand, foot, and mouth disease patients in Toyama, Japan between 1981 and 2007. Jpn. J. Infect. Dis. 62, 254–259. Li, R., Liu, L., Mo, Z., Wang, X., Xia, J., Liang, Z., et al., 2014a. An inactivated enterovirus 71 vaccine in healthy children. N. Engl. J. Med. 370, 829–837, http://dx.doi.org/10.1056/NEJMoa1303224. Li, B., Yadav, R.K., Jeong, G.S., Kim, H., Chae, H., 2014b. The characteristics of Bax inhibitor-1 and its related diseases. Curr. Mol. Med. 14, 603–615, http://dx.doi. org/10.2174/1566524014666140603101113. Liu, M., Zhao, J., Chen, K., Bian, X., Wang, C., Shi, Y., et al., 2012. G protein-coupled receptor FPR1 as a pharmacologic target in inflammation and human glioblastoma. Int. Immunopharmacol. 14, 283–288, http://dx.doi.org/10.1016/ j.intimp.2012.07.015. Liu, L., Zhang, Y., Wang, J., Zhao, H., Jiang, L., Che, Y., et al., 2013. Study of the integrated immune response induced by an inactivated EV71 vaccine. PLoS One 8, e54451, http://dx.doi.org/10.1371/journal.pone.0054451. Mi, H., Thomas, P., 2009. PANTHER pathway: an ontology-based pathway database coupled with data analysis tools. Methods Mol. Biol. 563, 123–140, http://dx. doi.org/10.1007/978-1-60761-175-2 7.

10

J. Song et al. / Virus Research 214 (2016) 1–10

Mi, H., Muruganujan, A., Casagrande, J.T., Thomas, P.D., 2013. Large-scale gene function analysis with the PANTHER classification system. Nat. Protoc. 8, 1551–1566, http://dx.doi.org/10.1038/nprot.2013.092. Morales-Garcia, G., Falfan-Valencia, R., Garcia-Ramirez, R.A., et al., 2012. Pandemic influenza A/H1N1 virus infection and TNF, LTA, IL1B, IL6, IL8, and CCL polymorphisms in Mexican population: a case-control study. BMC Infect. Dis. 12, 299. National Academy of Science, 2011. Guide for the Care and Use of Laboratory Animals. National Academy Press, Washington, DC. Ning, S., Pagano, J.S., Barber, G.N., 2011. IRF7: activation, regulation, modification and function. Genes Immun. 12, 399–414, http://dx.doi.org/10.1038/gene. 2011.21. Nouailles, G., Dorhoi, A., Koch, M., Zerrahn, J., Weiner 3rd., J., Faé, K.C., et al., 2014. CXCL5-secreting pulmonary epithelial cells drive destructive neutrophilic inflammation in tuberculosis. J. Clin. Invest. 124, 1268–1282, http://dx.doi.org/ 10.1172/JCI72030. Ohi, K., Hashimoto, R., Yasuda, Y., Nemoto, K., Ohnishi, T., Fukumoto, M., et al., 2012. Impact of the genome wide supported NRGN gene on anterior cingulate morphology in schizophrenia. PLoS One 7, e29780, http://dx.doi.org/10.1371/ journal.pone.0029780. Okino, C.H., dos Santos, I.L., Fernando, F.S., Alessi, A.C., Wang, X., Montassier, H.J., 2014. Inflammatory and cell-mediated immune responses in the respiratory tract of chickens to infection with avian infectious bronchitis virus. Viral Immunol. 27, 383–391, http://dx.doi.org/10.1089/vim.2014.0054. Pan, H., Zhu, Y.F., Qi, X., Zhang, Y.J., Li, L., Deng, B., et al., 2009. Analysis on the epidemiological and genetic characteristics of enterovirus type 71 and coxsackie A16 virus infection in Jiangsu, China. Zhonghua Liu Xing Bing Xue Za Zhi. 30, 339–343. Reboldi, A., Coisne, C., Baumjohann, D., Benvenuto, F., Bottinelli, D., Lira, S., et al., 2009. C-C chemokine receptor 6—regulated entry of TH-17 cells into the CNS through the choroid plexus is required for the initiation of EAE. Nat. Immunol. 10, 514–523, http://dx.doi.org/10.1038/ni.1716. Reggiori, F., Wang, C.W., Stromhaug, P.E., et al., 2003. Vps51 is part of the yeast Vps fifty-three tethering complex essential for retrograde traffic from the early endosome and Cvt vesicle completion. J. Biol. Chem. 278 (7), 5009–5020. Rodrigues, A.F., Carrondo, M.J.T., Alves, P.M., Coroadinha, A.S., 2014. Cellular targets for improved manufacturing of virus-based biopharmaceuticals in animal cells. Trends Biotechnol. 32, 602–607, http://dx.doi.org/10.1016/j.tibtech.2014.09. 010. Saha, B., Jyothi, S., Chandrasekar, B., Nandi, D., 2010. Gene modulation and immunoregulatory roles of interferon gamma. Cytokine 50, 1–14, http://dx. doi.org/10.1016/j.cyto.2009.11.021. Schmitt, L., Tampe, R., 2000. Affinity, specificity, diversity: a challenge for the ABC transporter TAP in cellular immunity. Chembiochem 1, 16–35. Schoenborn, J.R., Wilson, C.B., 2007. Regulation of interferon-gamma during innate and adaptive immune responses. Adv. Immunol. 96, 41–101, http://dx.doi.org/ 10.1016/S0065-2776(07)96002-2. Serti, E., Werner, J.M., Chattergoon, M., Cox, A.L., Lohmann, V., Rehermann, B., 2014. Monocytes activate natural killer cells via inflammasome-induced interleukin 18 in response to hepatitis C virus replication. Gastroenterology 147, 209–220, http://dx.doi.org/10.1053/j.gastro.2014.03.046. Shi, Y., He, X., Zhu, G., et al., 2015. Coxsackievirus A16 elicits incomplete autophagy involving the mTOR and ERK pathways. PLoS One 10 (4), e0122109. Solomon, T., Lewthwaite, P., Perera, D., Cardosa, M.J., McMinn, P., Ooi, M.H., 2010. Virology, epidemiology, pathogenesis, and control of enterovirus 71. Lancet Infect. Dis. 10, 778–790, http://dx.doi.org/10.1016/S1473-3099(10)70194-8.

Tung, W.H., Hsieh, H.L., Yang, C.M., 2010. Enterovirus 71 induces COX-2 expression via MAPKs, NF-kappaB, and AP-1 in SK-N-SH cells: Role of PGE(2) in viral replication. Cell Signal. 22 (2), 234–246. Wang, S., Liu, C., 2014. Update of enterovirus 71 infection: epidemiology, pathogenesis and vaccine. Expert Rev. Anti Infect. Ther. 12, 447–456, http://dx. doi.org/10.1586/14787210.2014.895666. Wang, S., Lei, H., Liu, C., 2012. Cytokine immunopathogenesis of enterovirus 71 brain stem encephalitis. Clin. Dev. Immunol. 2012, 1–8, http://dx.doi.org/10. 1155/2012/876241. Wang, B., Xi, X., Lei, X., et al., 2013. Enterovirus 71 protease 2Apro targets MAVS to inhibit anti-viral type I interferon responses. PLoS Pathog. 9 (3), e1003231. Wang, J., Qi, S., Zhang, X., Zhang, Y., Liu, L., Che, Y., et al., 2014. Coxsackievirus A 16 infection does not interfere with the specific immune response induced by an enterovirus 71 inactivated vaccine in rhesus monkeys. Vaccine 32, 4436–4442, http://dx.doi.org/10.1016/j.vaccine.2014.06.062. Weng, K., Chen, L., Huang, P., Shih, S., 2010. Neural pathogenesis of enterovirus 71 infection. Microb. Infect. 12, 505–510, http://dx.doi.org/10.1016/j.micinf.2010. 03.006. Wong, S.S.Y., Yip, C.C.Y., Lau, S.K.P., Yuen, K.Y., 2010. Human enterovirus 71 and hand, foot and mouth disease. Epidemiol. Infect. 138, 1071–1089, http://dx. doi.org/10.1017/S0950268809991555. Wright Jr., H.T., Landing, B.H., Lennette, E.H., McAllister, R.M., 1963. Fatal infection in an infant associated with coxsackie virus group A, type 16. N. Engl. J. Med. 268, 1041–1044, http://dx.doi.org/10.1056/NEJM196305092681904. Xiao, Z., Cilz, N.I., Kurada, L., Hu, B., Yang, C., Wada, E., et al., 2014. Activation of neurotensin receptor 1 facilitates neuronal excitability and spatial learning and memory in the entorhinal cortex: beneficial actions in an Alzheimer’s disease model. J. Neurosci. 34, 7027–7042, http://dx.doi.org/10.1523/ JNEUROSCI.0408-14.2014. Xu, C., He, X., Zheng, Z., et al., 2014. Downregulation of microRNA miR-526a by enterovirus inhibits RIG-I-dependent innate immune response. J. Virol. 88 (19), 11356–11368. Yen, F.B., Chang, L.Y., Kao, C.L., Lee, P.I., Chen, C.M., Lee, C.Y., et al., 2009. Coxsackieviruses infection in northern Taiwan—epidemiology and clinical characteristics. J. Microbiol. Immunol. Infect. 42, 38–46. Yu, J., Zhang, L., Ren, P., et al., 2015. Enterovirus 71 mediates cell cycle arrest in S phase through non-structural protein 3D. Cell Cycle 14 (3), 425–436. Zhang, Y., Cui, W., Liu, L., Wang, J., Zhao, H., Liao, Y., et al., 2011a. Pathogenesis study of enterovirus 71 infection in rhesus monkeys. Lab. Invest. 91, 1337–1350, http://dx.doi.org/10.1038/labinvest.2011. Zhang, W., Wang, W.Y., Yang, Z.H., Pang, B.D., Wu, H., Yang, Q.Z., et al., 2011b. Severe hand-foot-and-mouth disease caused by mixed infection of enterovirus 71 and coxsackieA16: reports of 6 cases. Chin. Gen. Pract. 2011, 3341–3343. Zhang, Y., Yang, E., Pu, J., Liu, L., Che, Y., Wang, J., et al., 2014. The gene expression profile of peripheral blood mononuclear cells from EV71-infected rhesus infants and the significance in viral pathogenesis. PLoS One 9, e83766, http:// dx.doi.org/10.1371/journal.pone.0083766. Zou, X.N., Zhang, X.Z., Wang, B., et al., 2012. Etiologic and epidemiologic analysis of hand, foot, and mouth disease in Guangzhou city: a review of 4,753 cases. Braz. J. Infect. Dis. 16 (5), 457–465. van Gent, M., Gram, A.M., Boer, I.G., et al., 2015. Silencing the shutoff protein of Epstein–Barr virus in productively infected B cells points to (innate) targets for immune evasion. J. Gen. Virol. 96, 858–865.