The Entomopathogenic Fungi Isaria fumosorosea Plays a Vital Role in Suppressing the Immune System of Plutella xylostella: RNA-Seq and DGE Analysis of Immunity-Related Genes

Most, if not all, entomopathogenic fungi have been used as alternative control agents to decrease the insect resistance and harmful effects of the insecticides on the environment. Among them, Isaria fumosorosea has also shown great potential to control different insect pests. In the present study, we explored the immune response of P. xylostella to the infection of I. fumosorosea at different time points by using RNA-Sequencing and differential gene expression technology at the genomic level. To gain insight into the host-pathogen interaction at the genomic level, five libraries of P. xylostella larvae at 12, 18, 24, and 36 h post-infection and a control were constructed. In total, 161 immunity-related genes were identified and grouped into four categories; immune recognition families, toll and Imd pathway, melanization, and antimicrobial peptides (AMPs). The results of differentially expressed immunity-related genes depicted that 15, 13, 53, and 14 up-regulated and 38, 51, 56, and 49 were down-regulated in P. xylostella at 12, 18, 24, and 36 h post-treatment, respectively. RNA-Seq results of immunity-related genes revealed that the expression of AMPs was reduced after treatment with I. fumosorosea. To validate RNA-Seq results by RT-qPCR, 22 immunity-related genes were randomly selected. In conclusion, our results demonstrate that I. fumosorosea has the potential to suppress the immune response of P. xylostella and can become a potential biopesticide for controlling P. xylostella.


INTRODUCTION
Insects are surrounded by an environment rich with harmful microorganisms and recurring infections are common in the natural environment. In order to combat these potentially infectious pathogens, insects have evolved various defense systems, including the potent immune system. Unlike mammals, insects solely rely on innate immune responses for host defense. The innate immune responses are usually comprised of cellular and humoral defense responses. The former is best demonstrated by the action of hemocytes in the phagocytosis (Kanost et al., 2004) whereas the hallmark of latter is the synthesis of antimicrobial peptides (AMPs) (Hoffmann and Reichhart, 2002). Upon microbial infection, cellular, and humoral responses are activated by insects, to clear the infection, through different steps (Söderhäll and Cerenius, 1998). The invading pathogen is recognized via pattern recognition receptors (PRRs) (Hultmark, 2003) leading to the amplification of signal of infection by serine proteases following the activation of signaling pathways (Jiang and Kanost, 2000;Osta et al., 2004). Finally, the effector factors are induced in the specific tissues to combat the pathogens.
To counter the defense system of the host, insect pathogenic fungi have also developed their mechanisms. The pathogens use a set of enzymes to breach the cuticle (Butt, 2002) and also release secondary metabolites, to suppress the immune system of the host, during colonization (Vilcinskas et al., 1997;Vey et al., 2002). Among these entomopathogenic fungi, on one hand, Metarhizium anisopliae has developed a new technique to evade the immune system of host via masking the cell wall during hemocoel colonization (Wang and Leger, 2006), and on the other hand, Isaria fumosorosea releases chitinase, chitosanase, lipase, to physically penetrate the host and suppress its regulatory system, and a beauvericin compound to paralyze the host (Hajek and St. Leger, 1994;Ali et al., 2010).
The diamondback moth (DBM), Plutella xylostella (L.) (Lepidoptera: Plutellidae), is one of the devastating pests of brassicaceous crops and costs approximately US$4-5 billion per year worldwide (Zalucki et al., 2012). P. xylostella is commonly known to rapidly evolve resistance against almost all types of insecticides including products of Bacillus thuringiensis (Shakeel et al., 2017). Consequently, entomopathogenic fungi have received an increased attention as an environmentally friendly alternative control measure to insecticides for controlling P. xylostella. Several strains of fungi have been isolated and used to control various insect pests including P. xylostella (Altre et al., 1999;Leemon and Jonsson, 2008;Bukhari et al., 2011). Of these entomopathogenic fungi, I. fumosorosea is considered as one of the promising species of fungi to be used as biological control of insect pests and various mycopesticide based on I. fumosorosea have been developed worldwide (Zimmermann, 2008). Isaria fumosorosea, a well-known entomopathogenic fungi, is distributed worldwide. It was previously known as Paecilomyces fumosoroseus, however, now it has been transferred to Isaria genus (Zimmermann, 2008). Due to wide host range, it has become a promising biological control agent and its potential as a biological control agent, other than immunity, has been tested to control various insect pests, including Diaphorina citri (Avery et al., 2011), Bemisia tabaci , Trialeurodes vaporariorum (Gökçe and Er, 2005), and P. xylostella .
Previously, most of the reports on insect immunity preferred model insects, including Drosophila melanogaster (Wraight et al., 2010), Manduca sexta (Kanost et al., 2004), and Tenebrio molitor  against insect pathogenic fungi such as M. acridium and Beauveria bassiana (Xiong et al., 2015;Zhang et al., 2015). It is only recently that P. xylostella immunity has received the attention of few researchers, thanks to the availability of the genome sequence of P. xylostella (You et al., 2013). Although, a recent report on the immune response of P. xylostella to B. bassiana improved our information of insectpathogen interaction (Chu et al., 2016). However, the changes that occur in response to I. fumosorosea in P. xylostella are largely unclear, restricting the development of fungal species other than M. anisopliae and B. bassiana to be adopted as a biological control agent in P. xylostella and other lepidopteran pests' control.
To gain deep insight into the immunogenetics of P. xylostella, the present study conducted a genome-wide profiling analysis of I. fumosorosea challenged P. xylostella larvae at 12, 18, 24, and 36 h post-infection using RNA-Seq and digital gene expression (DGE). Additionally, a global survey of the activities of antifungal immune defense genes in P. xylostella may also contribute to the in-depth analysis of candidate genes in P. xylostella immunity.

Insect Stock
The population of susceptible P. xylostella was kindly provided by Institute of Plant Protection, Guangdong Academy of Agricultural Sciences, China and was maintained in the Engineering Research Centre of Biological Control Ministry of Education, South China Agricultural University, Guangzhou, Guangdong province, P. R. China for five years without exposure to pesticides. Larvae were maintained at 25 ± 1 • C with a light: dark cycle of 16:8 h and 60-70% relative humidity.

Fungus Culture, Conidia Suspension Preparation, and Samples Collection
The I.fumosorosea IfB01 strain (China Center for Type Culture Collection access number: CCTCC M 2012400) was cultured on a potato dextrose agar (PDA) plate at 26 • C. The conidia were collected from 10 days old culture and suspended with 0.05% Tween-80 into standardized 1 × 10 8 spores/mL . Healthy P. xylostella larvae (third-instar) were selected and separated into two groups. One group (treatment) was treated with the 1 × 10 7 spores/ mL suspension, whereas the other group (control) was treated with sterile deionized water containing 0.05% Tween-80. The samples of 50 surviving larvae were collected from the treatment group and the control group at 12, 18, 24, and 36 h, respectively, forming three pairs of hour post-treatment infection and hours post treatment control. Different time-points of sampling were selected to observe infection dynamics (Abkallo et al., 2015) and dynamical changes (Bar-Joseph et al., 2012) of differentially expressed genes (DEGs) in response to Isaria fumosorosea in Plutella xylostella, as the gene expression profiling of different time points can provide DEGs dynamical behavior information.
from each treatment and control for the isolation of poly (A) + mRNA using oligo (dT) magnetic beads. cDNAs (First-and second-strand) were prepared using random hexamers, RNase H, and DNA polymerase I. Magnetic beads were used to purify the double strand cDNA and finally, ligation of fragments was carried out with sequencing adaptors. To quantify and qualify the libraries of samples, Agilent 2100 Bioanalyzer and ABI Step One Plus Real-Time PCR System were employed and then sequencing was done on the Illumina HiSeq TM 2000 system (Illumina, USA). Illumina sequencing was carried out at the Beijing Genomics Institute (BGI-Shenzhen, China).

Mapping and Functional Analysis of Differentially Expressed Genes
The process of filtration was performed in such a way that raw reads with adopters and unknown bases (more than 10%) were removed. After filtering, the remaining clean reads were mapped to reference gene using Bowtie (Langmead et al., 2009) and HISAT (Kim et al., 2015) was employed to map the reference genome. Finally, normalization of all data was done as fragments per kilobase of transcript per million fragments mapped (FPKM). The analysis of differential expression was employed by a rigorous algorithm. The false discovery rate (FDR) methodology was adopted in multiple tests (Kim and van de Wiel, 2008) for determination of threshold of P-value. Genes with significant differential expression were searched out according to a standard threshold having an FDR value of <0.001 and the absolute value of log2 ratio ≥ 1. The genome database of P. xylostella was used as the background to determine significantly enriched GO terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enriched within the DEG dataset using hypergeometric test and a corrected P-value (≤0.05) as a threshold.

Validation of DEGs Libraries by RT-qPCR
In order to validate mRNA expression levels exhibited by RNA-Seq results, Real-time quantitative PCR (RT-qPCR) was performed with 22 immunity-related DEGs chosen from the comparison of control vs. treatments. Total RNA isolation method was same as described earlier. In total, 1 µg of total RNA was treated with DNaseI (Fermentas, Glen Burnie, MD, USA) according to the instructions of the manufacturer and then cDNA was prepared using M-MLV reverse transcriptase (Promega, USA). The RT-qPCR was performed on a Bio-Rad iQ2 optical system (Bi-Rad) using SsoFast EvaGreen Supermix (Bio-Rad, Hercules, CA, USA) according to guidelines provided by the manufacturer. To confirm the PCR products purity, the amplification cycling parameters were set as; 95 • C for 30 s, 40 cycles of 95 • C for 5 s, and 55 • C for 10 s with a dissociation curve generated from 65 to 95 • C (Shakeel et al., 2015). For normalization, ribosomal protein S13 (RPS13) was used as an internal control (Fu et al., 2013) and the relative expression of genes was calculated using the 2 − CT method (Livak and Schmittgen, 2001). The primers used for RT-qPCR are listed in Table 1.

Features of the Sequenced cDNA Libraries
To identify genes involved in P. xylostella immunity in response to I. fumosorosea, five cDNA libraries were constructed from 3rd larval instar of P. xylostella at 12, 18, 24, 36 h after fungal treatment and control. A total of 11,652,857, 11,819,310, 12,051,947, 11,744,46, and 11,683,647 reads were generated from these five libraries (12, 18, 24, 36 h, and control respectively), from which 70. 01, 73.55, 73.23, 70.11, and 71.94% reads could be successfully mapped to the reference genome ( Table 2).

Dynamics of Differentially Expressed Immunity-Related Genes in Response to I. fumosorosea
To study the gene expression of P. xylostella larvae infected with I. fumosorosea, the pairwise comparison was carried out between libraries to determine the DEGs. The analysis of five libraries was carried out by determining the number of fragments per kb per million (FPKM) of clean reads. Relative to control, genes with (FDR) ≤ 0.001 and |log2Ratio| ≥ 1 were recognized as differentially expressed. Our results exhibited that, compared to the control, there were 53 (15 up-regulated and 38 downregulated), 64 (13 up-and 51 down-regulated), 109 (53 upregulated and 56 down-regulated), and 63 (14 up-and 49 down-regulated) immune-related genes that were significantly changed in P. xylostella after 12, 18, 24, and 36 h, respectively (Figure 1). A Venn diagram analysis showed that only 11 immunity-related DEGs were commonly expressed among all the treatments, whereas 7, 13, 45, and 12 immunity-related DEGs were specifically expressed in 12, 18, 24, and 36 h, respectively (Figure 2).

Go and KEGG Classification and Enrichment Analysis of Immunity-Related Genes in Response to I. fumosorosea
Following GO annotation, the immunity-related genes were classified into 26 different groups belonging to biological process, cellular component, and molecular function categories (Figure 3). In the biological process category, the two most enriched groups were the response to stimulus and biological regulation, whereas membrane and regulation of biological process were the top two enriched groups in the cellular component. The number of genes involved in catalytic activity and binding were the dominant groups in the category of   molecular function (Figure 3). The KEGG classification system categorized immunity-related genes into 21 different groups (Figure 4). The top five enriched groups among KEGG categories included infectious diseases (viral), signaling molecules and interaction, digestive system, infectious diseases (parasitic), and signal transduction (Figure 4).

Verification of DEG Results by RT-qPCR
To validate DEG results, 13 randomly selected immunity-related DEGs were analyzed by RT-qPCR. The results exhibited that the trend of expression level for all the selected genes was in consistence to that of RNA-Seq ( Figure 5).

Identification of Immunity-Related Genes
To identify immunity-related genes in response to I. fumosorosea, we searched out the genome of P. xylostella and combined BLAST search and GO annotation results. A number of genes having fold change less than one and those annotated as hypothetical or unknown proteins were not selected. Finally, a good number (161) of immunity-related genes were identified and classified as immune recognition families, toll and Imd signaling pathways, melanization, AMPs, and others ( Table 3).
The entomopathogenic fungi are recognized as an environmentally friendly tactic for controlling the insect pests. Previously, the entomopathogenic fungi like M. anisopliae and B. bassiana have received an increasing attention due to wide host range and capability of mass production (Butt et al., 2001). Recently, it has been shown that I. fumosorosea also has the potential to control various insect pests (Gökçe and Er, 2005;Huang et al., 2010;Avery et al., 2011). Therefore, considering the importance of I. fumosorosea, a genomic analysis of immune response of P. xylostella following infection of I. fumosorosea at different time points using high-throughput sequencing Illumina was performed in the present study.

Immune Recognition Families
Recognition of pathogen is the initial step in the defense against invading microbes, eliciting cellular and humoral responses. Pathogens produce conserved pathogen-associated molecular patterns (PAMPs) and the host produces pattern-recognition receptors (PRRs) in response (Mogensen, 2009). PRRs like peptidoglycan recognition proteins (PGRPs), β -Glucan binding proteins (GNBPs), lectins, scavenger receptors, and hemolin bind to the PAMPs (Hultmark, 2003). Insect PGRPs can trigger signal transduction through the Toll pathway, leading to the activation of AMP production (Zaidman-Rémy et al., 2011). In the present report, 14 PGRPs were identified and most of them were down-regulated after treatment with I. fumosorosea (Figure 6 and Table 3). Among the down-regulated PGRPs, two PGRP transcripts (px_105387866 and px_105386207) were down-regulated up to 2-fold (−2.60 and −2.21), respectively at 12 h post-treatment. Previously, it has also been shown that PGRPs were down-regulated after the injection of secondary metabolite (destruxin) of M. anisopliae in D. melanogaster (Pal et al., 2007), whereas in contrast, PGRPs were up-regulated in response to M. acridium and Beauveria bassiana in Helicoverpa armigera and Ostrinia furnacallis (Liu et al., 2014;Xiong et al., 2015). The expression of other PRRs like lectins, hemolin, and GNBPs was also down-regulated after treatment with I. fumosorosea (Figure 6 and Table 3). Among PRRs, only the scavenger receptors were up-regulated at all-time points postinfection. Our results are in accordance with a previous report showing that among PRRs, only scavenger receptors were upregulated in response to destruxin A in D. melanogaster (Pal et al., 2007). Our results suggest that PRRs like PGRPs, GNBPs, lectins, and hemolin may be the target of I. fumosorosea and scavenger receptors are responsible for the activation of the immune response of P. xylostella to I. fumosorosea.

Toll and Imd Signaling Pathways
The Toll pathway is primarily activated by fungi and Grampositive bacteria while the Gram-negative bacteria triggers the activation of Imd pathway leading to the production of AMPs (Aggarwal and Silverman, 2008;Hetru and Hoffmann, 2009). Here, in our study, we found that only spatzle and MyD88 showed differential expression while the other immune genes of toll pathway were not induced after treatment with I. fumosorosea (Figure 6 and Table 3). Of note, Imd pathway was also not induced after the treatment with I. fumosorosea at different time points. The expression of MyD88 was upregulated whereas, spatzle showed down-regulated expression after treatment (Figure 6 and Table 3). Previously, a similar phenomenon was observed in D. melanogaster where only pelle and toll showed differential expression in the Toll pathway, and Imd pathway was not induced after injection of destruxin A (Pal et al., 2007). Thus, our results show that I. fumosorosea has the ability to suppress the expression of toll pathway genes and in the meantime P. xylostella could resist the infection of I. fumosorosea.

Melanization
Melanization is considered as a vital component of the immune system of insects. It regulates the melanin cascade mediated by prophenoloxidases (PPO) (Taft et al., 2001). When a pathogen invades, PPO gets activated and transformed into PO following transformation of phenolic substances into quinone intermediates and ultimately killing pathogens. Here, only three PPO were found after the treatment of P. xylostella with I. fumosorosea and two of them were up-regulated up to 2-fold at 12 h post-infection.
Serine proteases represent a very large group of enzymes in almost all organisms and are involved in various biological processes (Ross et al., 2003). The structure of serine proteases consists of His, Asp, and Ser amino acid residues forming a catalytic triad (Perona and Craik, 1995). Generally, serine proteases are inactive pro-enzymes and need proteolytic cleavage for activation (Ross et al., 2003). Notably, many serine proteases identified in our study showed up-and down-regulated expression with a serine protease (px_105393891) showing highly up-regulated expression (10.77) and a serine protease (px_105381636) showing downregulated expression (−9.26) after treatment with I. fumosorosea at 18 h post-infection (Figure 6 and Table 3). It has been reported that the serine proteases showed same up-and down-regulated expression in P. xylostella and D. melanogaster after treatment with destruxin A (Pal et al., 2007;Han et al., 2013).
Serpins, a super-family of proteins, are found in nearly all organisms (Gettins, 2002). In general, they contain 350-400 amino acid residues. The similarity of amino acid sequence ranges from 17 to 95% among all serpins. They contain three β-sheets and seven to nine α-helices folding into a conserved tertiary structure with a reactive center loop (RCL) (Gettins, 2002). The RCL of these serpins binds to the specific active site of the target proteinase. When the cleavage of the serpin takes place at scissile bond, it goes through an important conformational change, trapping the target proteinase covalently (Dissanayake et al., 2006;Ulvila et al., 2011). Interestingly, almost all the serpins were down-regulated at an early stage of treatment at 12 h post-infection. In contrast, the expression level of serpins was up-regulated in P. xylostella after treatment with Diadegma semiclausum parasitization (Etebari et al., 2011;Han et al., 2013). The activation of serpins by D. semiclausum in previous reports may be a strategy to suppress the activity of PPO in the host defense system.

Antimicrobial Peptides
AMPs are evolutionarily conserved low molecular weight proteins and play a vital part in the insect defense system against microorganisms (Bulet and Stocklin, 2005). Here, in the present study, lysozyme, moricin, gloverin, and cecropin were identified after the treatment of P. xylostella with I. fumosorosea at different time periods. Interestingly, all the AMPs were down-regulated after treatment with I. fumosorosea (Figure 6 and Table 3) The expression of lysozyme (px_105381977) was decreased up to 10-fold (−10. 87) at 12 h post-infection, moricin (px_105392532) expression was decreased up to 9-fold (−9.57)    at 12 h post-infection, gloverin (px_105389810) expression was reduced up to 4-fold (−4.80) at 18 h post-infection, and the expression of cecropin (px_105394859) was down-regulated up to 6-fold (−6.03) at 18 h post-infection of I. fumosorosea. Previously, most of the reports on immune response of insects to entomopathogenic fungi identified that the expression of AMPs was up-regulated after the treatment leading to a conclusion that the entomopathogenic fungi were unable to suppress the immune system (Liu et al., 2014;Xiong et al., 2015;Zhang et al., 2015). However, varroa mites and destruxin A were reported to suppress the expression of AMPs in Apis mellifera and D. melanogaster (Gregory et al., 2005;Pal et al., 2007). The immune response suppression in host by an entomopathogenic fungi such as, I. fumosorosea would have obvious benefits for success of pathogenic fungi. Previously, it was observed that when mutations were introduced in Toll and IMD pathways, the D. melanogaster was unable to produce AMPs resulting in extreme vulnerability to fungal challenge (Lemaitre et al., 1996;Tzou et al., 2002). Thus, the ability to reduce AMP production is likely to aid fungal survival in a variety of insect hosts. A similar suppression of AMPs in our study by I. fumosorosea adds a new dimension to the dynamics of host-pathogen interactions.

CONCLUSION
Concluding our findings, the present study adopted genomic analysis with RNA-Seq and DGE technology to find out DEGs especially focusing on immunity-related DEGs after treatment with I. fumosorosea. It is speculated that the entomopathogenic fungi I. fumosorosea not only down-regulated the expression of PRRs and other immune genes but also the activity of AMPs was inhibited leading to the ultimate suppression of the immune system of P. xylostella. Thus, it shows that I. fumosorosea has the potential to suppress the immune system of P. xylostella and can be adopted as a bio-pesticide for P. xylostella control. Our study explores a new avenue in research to develop bio-pesticides for controlling P. xylostella by focusing on the insect immune system.

ETHICS STATEMENT
Our work confirms to the legal requirements of the country in which it was carried out.