Transcriptome Analysis of the Midgut of the Chinese Oak Silkworm Antheraea pernyi Infected with Antheraea pernyi Nucleopolyhedrovirus

The Antheraea pernyi nucleopolyhedrovirus (ApNPV) is an exclusive pathogen of A. pernyi. The intense interactions between ApNPV and A. pernyi cause a series of physiological and pathological changes to A. pernyi. However, no detailed report exists regarding the molecular mechanisms underlying the interactions between ApNPV and A. pernyi. In this study, four cDNA libraries of the A. pernyi midgut, including two ApNPV-infected groups and two control groups, were constructed for transcriptomic analysis to provide new clues regarding the molecular mechanisms that underlie these interactions. The transcriptome of the A. pernyi midgut was de novo assembled using the Trinity platform because of the lack of a genome resource for A. pernyi. Compared with the controls, a total of 5,172 differentially expressed genes (DEGs) were identified, including 2,183 up-regulated and 2,989 down-regulated candidates, of which 2,965 and 911 DEGs were classified into different GO categories and KEGG pathways, respectively. The DEGs involved in A. pernyi innate immunity were classified into several categories, including heat-shock proteins, apoptosis-related proteins, serpins, serine proteases and cytochrome P450s. Our results suggested that these genes were related to the immune response of the A. pernyi midgut to ApNPV infection via their essential roles in regulating a variety of physiological processes. Our results may serve as a basis for future research not only on the molecular mechanisms of ApNPV invasion but also on the anti-ApNPV mechanism of A. pernyi.


Introduction
Antheraea pernyi is not only an important economic insect but also a rising model organism because of a variety of advantages, such as its ease of rearing and experimental manipulation compared with other Lepidoptera insects. However, A. pernyi nucleopolyhedrovirus (ApNPV), a major viral pathogen for A. pernyi, causes enormous damage to A. pernyi and also the sericulture industry. ApNPV belongs to the NPV subfamily, which displays two types of phenotypes with the same genetic resource during its infection cycle, including budded virus (BV) and occlusion-derived virus (ODV). ODV has high infectivity to the intestinal epithelial cells of the host. The ODV particles invade the peritrophic membrane of the midgut by oral infection, and drop off their capsules through an interaction between the microvilli of the intestinal epithelial cells and the capsule. The infection is initiated with the nucleocapsid of ODV releasing into the intestinal epithelial cells [1]. Then, the host initiates resistance strategies to block viral infection, such as through the immune system, apoptosis, and RNA interference [2]. The infected A. pernyi larva exhibits several symptoms in the late stage of infection, such as no feeding, various sizes of circular pus blotches on the epidermis, abnormal behavior, and finally death. However, no detailed report exists regarding the molecular mechanisms underlying these interactions between ApNPV and A. pernyi.
Changes in host gene expression in response to viral infection are of great interest. Research has been focused on gene expression changes in the midgut in insects because of its important role in resisting pathogen invasion and proliferation as the first line of innate immunity in the host [3]. In the midgut of Aedes aegypti, silencing of 5G1, which belongs to the midgut serine protease family, or the soybean trypsin inhibitor significantly increased midgut infection rates of Dengue virus 2 (DENV-2), suggesting that some midgut serine proteases may limit DENV-2 infectivity of A. aegypti [4]. Smartt et al. (2009) identified 26 differentially expressed cDNA clones in the midgut tissue of Culex pipiens quinquefasciatus after exposure to West Nile virus (WNV) [5]. Using Digital Gene Expression analysis, Gao et al. (2014) identified 752 differentially expressed genes, including 649 up-regulated and 103 down-regulated genes, in the B. mori cytoplasmic polyhedrosis virus (BmCPV) infected midguts of 4008 silkworm strains [6]. Kolliopoulou et al. (2015) obtained 308 differentially expressed genes in Bombyx mori larval midgut infected with BmCPV through transcriptome analysis, some of which are involved in physical barriers and immune responses [7].
RNA-sequencing (RNA-Seq) based on deep sequencing technologies is a powerful and highthroughput method for transcriptome analysis. RNA-Seq has been widely used to investigate the molecular mechanisms of the interactions between pathogens and their insect hosts [8][9][10][11]. In this study, RNA-Seq was applied to polyadenylate-enriched mRNAs from A. pernyi midguts to better understand the complexity of the molecular mechanisms underlying the interactions between ApNPV and the A. pernyi midgut at the transcriptional level. We first assembled the transcriptome sequences and then identified differentially expressed genes (DEGs), including both up-regulated and down-regulated genes in ApNPV-infected midguts compared with controls. We also screened immune-related genes from the DEGs via bioinformatic analysis. The results of this study provide useful information for further research on the molecular mechanisms underlying the ApNPV-induced nuclear polydedrosis virus disease of A. pernyi.

Sample preparation
The A. pernyi strain Jiaolan stored in our laboratory was used in this study. The larvae were maintained in a rearing chamber at 23 ± 2°C with 70 ± 5% relative humidity and fed with fresh leaves of Quercus mongolic collected from the research base of Shenyang Agricultural Uinversity (Shenyang city, 41.8 N°, 123.4°E). After molting twice, the third instar A. pernyi larvae were separated into two groups randomly. The experimental group marked as Ap_NPV was fed with the leaves added with 4.05×10 6 polyhedra/mL for three days and then fed with fresh leaves. For the non-infected controls, the same volume of 0.9% physiological saline was mixed in the feed for the A. pernyi larvae, and the rearing conditions were identical to those of the ApNPV-infected groups. Each group contained twenty five A. pernyi larvae. The midguts were dissected from the two groups when polyhedra could be observed in the hemocytes of ApNPV-infected A. pernyi larvae under a microscope. All experiments were performed with two independent biological replicates.

Total RNA extraction
Total RNA was extracted using TRIzol 1 Reagent (Invitrogen) according to the manufacturer's protocol. The RNA was quantified by measuring the absorbance at 260 nm using a NanoVue UV-Vis spectrophotometer (Bio-Science). RNA purity was checked using a NanoPhotometer 1 spectrophotometer (IMPLEN, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

cDNA library construction and Illumina RNA-Seq
Library construction and RNA-Seq were performed by the Novogene Experimental Department (Beijing, China). Briefly, a total of 3 μg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using the NEBNext 1 Ultra™ RNA Library Prep Kit for Illumina 1 (NEB, USA) following manufacturer's instructions. The mRNA was purified using poly-T oligo-attached magnetic beads. Fragmentation was performed using divalent cations under elevated temperature in 5×NEB Next First Strand Synthesis Reaction Buffer. First strand cDNA was synthesized using a random hexamer primer and M-MuLV Reverse Transcriptase (RNase H -). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. The remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of the 3'ends of the DNA fragments, NEBNext Adaptors with hairpin loop structures were ligated to prepare for hybridization. To select cDNA fragments of preferentially 150-200 bp in length, the library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, USA). Then, 3 μL USER Enzyme (NEB, USA) was added to size-selected, adaptor-ligated cDNA at 37°C for 15 min followed by 5 min at 95°C before PCR. PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and an Index (X) Primer. Finally, PCR products were purified (AMPure XP system) and library qualities assessed on an Agilent Bioanalyzer 2100 system. Clustering of the index-coded samples was performed on a cBot Cluster Generation System using the TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2500 platform and paired-end reads were generated.

De novo transcriptome assembly
Reads with only adaptors, > 5% unknown nucleotides, or low quality reads were first filtered. The quality reads were assembled into unigenes using the Trinity software [12]. The gene expression levels were estimated by RSEM [13] for each sample. Firstly, clean data were mapped back onto the assembled transcriptome. Then, read count for each gene was obtained from the mapping results. The FPKM measure (Fragments per Kilobase per Millions base pairs sequenced) [14,15] was used to calculate unigene expression. The transcripts spliced by Trinity were used as reference sequences (ref), and the clean reads of each sample were mapped onto the ref using the RSEM software. The read count of each unigene in each sample was converted to FPKM to obtain gene expression. A FPKM threshold of > 0.1 was applied to ensure that the genes were expressed in both of the two groups. Differential expression analysis was performed using the DESeq R package (1.10.1). The resulting P values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted Pvalue < 0.05 found by DESeq were assigned as differentially expressed.

Validation of data reliability by quantitative RT-PCR
Quantitative RT-PCR (qRT-PCR) was used to confirm the expression profiles of genes that were identified from the Illumina sequencing analysis. Eight DEGs, including 5 down-regulated and 3 up-regulated genes, were selected for qRT-PCR. The gene-specific primers were designed using the predicted CDSs as reference sequences. The details of the primers are provided in S1 Table. qRT-PCR was performed on a LightCycler 480 Real-time Detection System (Roche Diagnostics). SYBR 1 Premix Ex TaqTMII (Tli RHaseH Plus) was mixed with 1 ng of template cDNA. The reaction was performed in a total volume of 10 μL, containing 5 μL of 2×SYBR Premix Ex Taq TM (TaKaRa), 1 μL of diluted cDNA mix, 0.4 μL of each primer (10 mM) and 3.2 μL of Milli-Q water. Cycling conditions were as follows: 30 s at 95°C, 5 s at 95°C, followed by 40 cycles with denaturation for 20 s at 60°C and annealing/elongation for 15 s at 65°C. Melting curves were generated after each run to confirm a single PCR product. All reactions were performed in triplicate. The housekeeping gene actin3 was used as the endogenous control. After the PCR program, data were analyzed with Light Cycler 480 software (Applied Biosystems). The comparative CT method (2 -ΔΔCT method) [16] was used to analyze the expression levels of different genes.

Ethics Statement
The leaves of Quercus mongolic were collected from the research base of Shenyang Agricultural University. We obtained the permission from Shenyang Agricultural University.

Transcriptome analysis of larval midgut samples
To obtain a global view of the transcriptome related to the response of the A.pernyi larvae midgut to ApNPV infection, we constructed four cDNA libraries, numbered correspondingly as Ap_CK1/ Ap_NPV1 and Ap_CK2/Ap_NPV2, in which Ap_NPV refers to the tested samples and Ap_CK the control ones with two replicates for each group. The four cDNA libraries (Ap_CK1, Ap_CK2, Ap_NPV1 and Ap_NPV2) produced 71,989,178, 57,985,902, 70,871,812 and 62,299,960 raw 150 bp paired-end reads (We submitted these raw reads to NCBI Sequence Read Archive, and obtained their NCBI accession numbers: SRR2919240, SRR2919241, SRR2919242 and SRR2919243, respectively) with GC percentages of 45.59%, 45.73%, 44.32% and 43.87%, respectively. After removal of adapter sequences, low quality reads and N-containing reads, 61,595,738 (85.56%), 49,324,114 (85.06%), 60,792,564 (85.78%) and 52,512,994 (84.29%) clean reads were obtained in the four cDNA libraries, respectively (Table 1). We identified 5,172 differentially expressed genes (DEGs) from the ApNPV-infected groups compared with the controls, including 2,183 up-regulated (S2 Table) and 2,989 down-regulated genes (S2 Table). The distribution of the DEGs is shown in Fig 1.

Analysis of the pearson correlation between different samples
The pearson correlation between different samples is an important parameter to test the reliability of the experiment and whether the sample selection is reasonable. In this study, we carried out two independent technology replications in each group. The pearson correlation between the two control samples was 0.902 with that was 0.84 between the two ApNPV infected samples, showing the high repeatability. The control and ApNPV infected samples showed low pearson correlation with each other, indicating the significant difference between them (Fig 2).

Gene ontology (GO) analysis of the differentially expressed genes
The DEGs were assigned to various GO categories to determine their functional classifications (Fig 3). A total of 2,965 DEGs were classified into different GO categories, including 1,821  In biological process, most DEGs were involved in transport, establishment of localization, localization and single-organism transport. In cellular component, the DEGs were enriched in cytoplasm, organelle membrane and cytoplasmic part. In molecular function, most GO terms were involved in enzyme activity, such as oxidoreductase activity, peptidase activity, serine hydrolase activity and catalytic activity, indicating that the enzyme system of A. pernyi midgut may be degenerated by ApNPV infection.   Table) and 154 pathways for the up-regulated genes (S7 Table). The pathways were classified into the following four categories: metabolism, genetic information processing, environmental information processing and cellular processes (Fig 4). In the mapped pathways of the down-regulated genes, the abundant genes mapped onto protein processing in endoplasmic, oxidative  phosphorylation, carbon metabolism and biosynthesis of amino acids. The up-regulated genes were mainly enriched in spliceosome, endocytosis, ubiquitin-mediated proteolysis and several signaling pathways (Fig 5).

Validation of Illumina sequencing results by qRT-PCR
To verify the RNA-seq results, we selected 8 DEGs, including 2 genes encoding HSPs, 1 gene encoding apoptosis-related protein, 2 genes encoding serine proteases, 2 genes encoding serpins and 1 gene encoding CYP, by designing specific primers for quantitative RT-PCR. The qRT-PCR results revealed different expression trends for the analyzed genes (Fig 6). Among the 8 DEGs, comp45066_c0, comp42863_c0, comp18400_c0, comp44655_c0 and comp37472_c0 displayed downward trends, whereas comp63819_c0, comp32864_c0 and comp39389_c0 exhibited rising trends. Besides, we also investigated the 25 top up-regulated and top down-regulated genes by qRT-PCR (S10 Table,

Discussion
Antheraea pernyi nucleopolyhedrovirus is an exclusive pathogen of A. pernyi. Like other NPVs, ApNPV also undergoes four temporally regulated transcription phases in the infection process: immediate early, early, late, and very late. The initial steps of viral infection mainly help the virus establish infection and selectively modulate certain categories of host genes. The movement of the nucleocapsid from the cytoplasm to the nucleus is essential for controlling host gene expression and producing viral progeny. The late transcription phase is focused mainly on DNA replication and structural gene expression [17,18]. In this study, to ensure that the A. pernyi larvae acquired the NPV disease during the 3 rd instar, a high dose of polyhedra (4.05×10 6 polyhedra/mL) was added to the meal to feed the A. pernyi larvae for three days.
The time point at which polyhedra can be observed in the hemocytes in ApNPV-infected larva under the microscope was selected to compare the genes in the ApNPV-infected midguts that were differentially expressed relative to the controls. A total of 5,172 DEGs, including 2,183 up-regulated and 2,989 down-regulated genes, were identified through transcriptome sequencing. These DEGs identified during infection provide potential insights into the complex molecular mechanisms of the response to viral infection at the transcriptional level. Up-regulated genes may represent host cell responses to pathogens, whereas down-regulated genes may be partially attributed to the shut off of host macro molecular synthesis in favor of viral replication [19]. Quantitative RT-PCR analysis was used to validate the differentially expressed genes according to RNA-seq. X-axis represents the eight genes selected for qRT-PCR validation. comp45066_c0, heat shock protein 60 (Chilo suppressalis); comp42863_c0, heat shock protein 19.9 (Bombyx mori); comp63819_c0, apoptosis-inducing factor 3 (Homo sapiens); comp18400_c0, serine protease 3 (Lonomia obliqua); comp44655_c0, serine protease 5 (Mamestra configurata); comp32864_c0, serine protease inhibitor 12 (Bombyx mori); comp39389_c0, serine protease inhibitor 5 (Bombyx mori); comp37472_c0, cytochrome CYP324A1 (Spodoptera littoralis). Y-axis represents the relative expression level. The Ct values of each reaction were normalized to the endogenous control actin3.
doi:10.1371/journal.pone.0165959.g006 The midgut is the first line of innate immunity in insects and plays an important role in defending against pathogen invasion. In this study, several homologous genes possibly implicated in A. pernyi's immune response against ApNPV infection were found by deep sequencing to be highly differentially expressed between the ApNPV-infected midguts and the controls. These genes are involved in several categories, including heat-shock protein, apoptosis, serpin, serine protease and cytochrome P450, as outlined in detail in Table 2. Together, these genes constitute a complex response to ApNPV infection in A. pernyi larvae.
Heat-shock proteins (HSPs) are expressed constitutively in all cells and essential for several important cellular processes, such as protein folding and protecting proteins from denaturation or aggregation [20]. HSPs play an important role in antimicrobial and autoimmune responses and have potent effects in inducing antigen-specific immunities to bound materials upon immunization [21]. When insects are invaded by pathogens such as virus, HSPs function as molecular chaperones in preventing the accumulation of damaged proteins to maintain cellular homeostasis by refolding, stabilization, intracellular translocation and degradation of proteins. HSPs can also act as a 'danger'-signaling molecule and send 'danger signals' to the immune system to generate a response to external stimuli, such as pathogens or environmental changes [21][22][23]. In this study, we identified 12 DEGs encoding different HSPs (Table 2), and the different expression trends of these genes indicate that some may be activated and others inhibited by ApNPV infection.
Apoptosis is the process of programmed cell death that is a normal component of the development of multicellular organisms. Apoptosis is required to destroy cells that may be infected with viruses or contain DNA damage that represent a threat to the integrity of the organism [24]. Apoptosis is one of the strategies by which antiviral defense mechanisms function in insects [25]. Some Lepidopteron insects resist baculovirus infection by sloughing off or selectively apoptosis the infected cells in the midgut epithelium [26]. We screened 7 DEGs involved in apoptosis from the obtained DEGs (Table 2). The up-regulation of the apoptosis-related genes (comp28678_c0, comp19334_c1 and comp45307_c0) and the down-regulation of the apoptosis inhibitors (comp18389_c0, comp42011_c0 and comp63250_c0) indicate that the host cells may initiate apoptosis to defend against the infection of ApNPV.
Serine proteases (SPs) are important enzymes that are not only involved in digestion, embryonic development and cell differentiation but also play important roles in resisting pathogen invasion [27]. SPs activate prophenoloxidase (PPO) and induce melanotic encapsulation of organisms by catalyzing PPO to active phenoloxidase (PO) and then initiate the innate immune response to pathogen invasion [28]. Serpins are serine protease inhibitors that exert tight regulation of proteolytic cascades important for many biological processes, such as the complement cascade, inflammation, and innate immunity in different organisms [29]. Serpins control the signals of innate immunity by regulating the activities of SPs to protect the host from pathogen infection [30]. By screening the DEGs, 18 SPs and 4 serpins were obtained ( Table 2). We observed that the transcript level of comp39389_c0, encoding serine protease inhibitor 5 (Serpin5), in the ApNPV-infected midgut was higher than that in the control, consistent with the results of Bao et al. (2009) and Wu et al. (2011) [31,32]. Serpin5 is a negative regulator of the Toll pathway and functions extracellularly, likely by blocking the proteolytic activation of Spaetzle, the Toll receptor ligand [33]. The differentially expressed SPs and serpins identified in this study indicate that they may play important roles in the immune response to ApNPV infection in the midgut of A. pernyi.
Cytochrome P450s (CYPs) are a superfamily of heme proteins that are involved in the metabolism a wide range of both endogenous and exogenous compounds [34]. In insects, CYPs are important for the detoxification of plant allelochemicals and insecticides [35,36]. Twelve CYP DEGs were identified in the current study (Table 2), and the different expression levels of these genes imply that they may participate in the immune response induced by ApNPV infection.
The results of KEGG analysis showed that the down-regulated genes were mainly involved in several metabolic pathways such as oxidative phsophorylation, carbon metabolism, the biosynthesis of amino acids, galactose metabolism and fatty acid degradation (Fig 5A, S6 Table), suggesting that the functions of digestion and absorption of the A. pernyi midgut were degenerated due to ApNPV infection. Among the mainly pathways related to the up-regulated genes, several important immune pathways involved in defending against the virus infection were found (Fig 5B, S6 Table), including PI3K-Akt signaling pathway, Hippo signaling pathway, and JAK-stat signal pathway, indicating that the A. pernyi midgut initiated the immune system after the virus invasion. In conclusion, our results provide new clues for exploring the molecular mechanisms of ApNPV infection and the anti-ApNPV mechanisms of A. pernyi.