Identification of Immunity-Related Genes in Dialeurodes citri against Entomopathogenic Fungus Lecanicillium attenuatum by RNA-Seq Analysis

Abstract

Dialeurodes citri is a major pest in citrus producing areas, and large-scale outbreaks have occurred increasingly often in recent years. Lecanicillium attenuatum is an important entomopathogenic fungus that can parasitize and kill D. citri. We separated the fungus from corpses of D. citri larvae. However, the sound immune defense system of pests makes infection by an entomopathogenic fungus difficult. Here we used RNA sequencing technology (RNA-Seq) to build a transcriptome database for D. citri and performed digital gene expression profiling to screen genes that act in the immune defense of D. citri larvae infected with a pathogenic fungus. De novo assembly generated 84,733 unigenes with mean length of 772 nt. All unigenes were searched against GO, Nr, Swiss-Prot, COG, and KEGG databases and a total of 28,190 (33.3%) unigenes were annotated. We identified 129 immunity-related unigenes in transcriptome database that were related to pattern recognition receptors, information transduction factors and response factors. From the digital gene expression profile, we identified 441 unigenes that were differentially expressed in D. citri infected with L. attenuatum. Through calculated Log₂Ratio values, we identified genes for which fold changes in expression were obvious, including cuticle protein, vitellogenin, cathepsin, prophenoloxidase, clip-domain serine protease, lysozyme, and others. Subsequent quantitative real-time polymerase chain reaction analysis verified the results. The identified genes may serve as target genes for microbial control of D. citri.

About the Authors:

Shijiang Yu

Affiliation: Citrus Research Institute, Southwest University/Chinese Academy of Agricultural Sciences, National Engineering Research Center for Citrus, Chongqing 400712, China

Lili Ding

Affiliation: Citrus Research Institute, Southwest University/Chinese Academy of Agricultural Sciences, National Engineering Research Center for Citrus, Chongqing 400712, China

Ren Luo

Affiliation: Citrus Research Institute, Southwest University/Chinese Academy of Agricultural Sciences, National Engineering Research Center for Citrus, Chongqing 400712, China

Xiaojiao Li

Affiliation: Citrus Research Institute, Southwest University/Chinese Academy of Agricultural Sciences, National Engineering Research Center for Citrus, Chongqing 400712, China

Juan Yang

Affiliation: Citrus Research Institute, Southwest University/Chinese Academy of Agricultural Sciences, National Engineering Research Center for Citrus, Chongqing 400712, China

Haoqiang Liu

Affiliation: Citrus Research Institute, Southwest University/Chinese Academy of Agricultural Sciences, National Engineering Research Center for Citrus, Chongqing 400712, China

Lin Cong

Affiliation: Citrus Research Institute, Southwest University/Chinese Academy of Agricultural Sciences, National Engineering Research Center for Citrus, Chongqing 400712, China

Chun Ran

* E-mail: [email protected]

Affiliation: Citrus Research Institute, Southwest University/Chinese Academy of Agricultural Sciences, National Engineering Research Center for Citrus, Chongqing 400712, China

Introduction

Dialeurodes citri Asahmaed (D. citri) belongs to the family Aleyrodidae, of the order Hemiptera. This important, widely spreading pest is found in the world’s citrus producing areas. It originates from Southeast Asia and has been found in Asia, South and North America, and Europe [1, 2]. D. citri has a piercing/sucking mouthpart, enabling them to suck the branches and leaves of citrus, which will cause leaf and fruit abscission. More seriously, D. citri can secrete honeydew, which can adsorb dust in the air and provide nutrition for some fungi. This leads to the occurrence of sooty mold and seriously affects photosynthesis and the quality of citrus fruits. In the state of Florida in the USA, D. citri was once the dominant pest for citrus trees [3]. In Oceania, D. citri was detected for the first time in New Zealand in 2000 [4], and it spread swiftly to other citrus-producing regions, leading to a reduction in output by about 90% [5]. The citrus-producing regions in Chongqing, China, such as Changshou, Tongnan, Tongliang, and Nanchuan, are at risk for large-scale outbreaks of D. citri, which bring much trouble to growers. At present, the main method for preventing D. citri infection is chemical control. However, it is difficult for a chemical pesticide to effectively control an outbreak of D. citri infection in a field. In addition, excessively frequent chemical spraying may lead to environmental pollution, pesticide residue, and drug resistance [6–8]. Therefore, an environmentally, friendly and effective control strategy for D. citri is needed.

Lecanicillium attenuatum (L. attenuatum), a member of the order Hypocreales, is an important entomopathogenic fungus. L. attenuatum parasitizes various pests, such as whiteflies, scale insects, aphids, and nematodes [9–12]. Recent research suggested that the fungus can also parasitize fungi, causing cucumber powdery mildew [13]. At present, L. attenuatum has been separated from crops in many countries [14–17]. Our laboratory successfully separated the fungus from D. citri corpses collected in a citrus nursery of Chongqing. The obtained fungus showed a strong infection ability in D. citri.

However, it is difficult for pathogenic fungi to infect D. citri. Pests, similar to those higher animals, have a complete immune defense system, which decreases greatly the lethality of pathogenic fungi. It is well known that insects possess three major lines of immune defense: the body wall, cell-mediated immunity, and humoral immunity. The body wall of insects is their first barrier against invaders. The main functional parts of the body wall are chitin, protein, wax, and some recognition factors. The cellular immunity of the insect relies mainly on the ability of the hemocytes in the body to enwrap and devour the invaders and antigens. In addition, hemocytes are also involved in wound healing and blood coagulation, acting to prevent pathogenic microorganisms from entering the body through a wound [18]. Gillespie et al. found that, after treatment by the Metarhizium acridum, the desert locust shows a sharp rise in the number of hemocytes and begins, playing a rold in the immune response [11]. This means the hemocytes respond rapidly, playing an important role in cellular immunity. The insect humoral immune response includes melanization, lysozyme, antimicrobial peptides (AMPs), lectins, antiviral factors (AvFs), and proteinase inhibitors (Pis). In case of microorganism infection, the insect can activate the cellular and humoral immune responses through four steps [19]: 1) identify the invaders. The insect has special pattern recognition receptors (PRRs), such as peptidoglycan recognition proteins (PGRPs), thioester-containing proteins (TEPs), gram-negative binding proteins (GNBPs), scavenger receptors (SCRs), C-type lectins (CTLs), and galectin (GALE), that can identify the external pathogenic microorganisms, leading to a downstream immune response [20]; 2) an extracellular cascade reaction to activate serine proteases and remove serine protease inhibitor; 3) amplification of the signal of infection or removal of the false alarms; 4) stimulation of the transcription of the effecting factor through a signal transduction pathway to produce the immune response. The AMPs are mainly subjected to the Toll pathway, Imd pathway, and Jas/Stat signal transduction pathway.

At present, genes and proteins related to the immune response is mainly studied in a model insect. However, there have been few studies on the identification of genes in D. citri that are responsible for the immune response to L. attenuatum infection. Our analyses of the transcriptome of D. citri (no reference genome) and differential expression of genes in the digital gene expression (DGE) profile provide information on immunity-related genes of D. citri and a theoretical basis for research into the molecular mechanism of the immunity of D. citri against pathogenic fungi. With gene function annotation, it is possible to find the major genes responsible for the defense of D.citri against L. attenuatum, and thus, the limitations and blindness of research on a single gene are reduced. This research will benefit the biological control of insects using pathogenic fungi.

Materials and Methods

Fungus culture and conidia suspension preparation

L. attenuatum strain TL001 was cultured on potato dextrose agar (PDA) plates at 25°C and 80% humidity. Conidia (spores) used for infection were harvested from 3–4 weeks old cultures by scraping the surface of the mycelia with sterile cell scrapers into sterile deionized water containing 0.1% Tween-80. Conidia were separated from other mycelial structures over a sterile funnel packed with autoclaved glass wool, washed twice with ddH₂O by centrifugation at 4,000 rpm, counted, and diluted to 1×10⁸ spores/ml. Freshly prepared conidia were used for all experiments.

Collection of D. citri and RNA extraction

D. citri eggs, nymphs of all ages, and adults were collected in the net house of D. citri and transferred to a 1.5-ml EP tube in equal proportions. The EP tube was numbered #1, frozen in liquid nitrogen, and stored at -80°C.

The citrus leaves to which the D. citri larvae were attached were collected and washed with clear water. These leaves were separated into two groups. One group was treated with the spore suspension for 3 days, whereas the other group, as the control group, was sprayed with sterile water. The successful infection of D. citri larvae with L. attenuatum was confirmed under a dissecting microscope by observation that the attached spores on the cuticle of larvae had germinated and enwrapped pests’ bodies. D. citri larvae in the control and treatment groups were picked and transferred to two 1.5-ml EP tubes, numbered #2 and #3, respectively. The contents of the EP tubes were frozen with liquid nitrogen and stored at -80°C.

Using the RNA Isolater Total RNA Extraction Reagent (Vazyme, China) method, total RNA was extracted from samples #1, #2 and #3. The genomic DNA in the total RNA was removed using DNase I (Invitrogen, USA). The quality of the extracted RNA was evaluated through 1% agarose gel electrophoresis, and the concentration was determined with an Nanodrop 2000N spectrophotometer (Thermo Fisher Scientific, USA), recorded as optical density ratios OD260/OD280 and OD260/OD230.

Construction and sequencing of cDNA library

The RNA (#1) with a poly(A) tail was purified with magnetic beads containing polynucleotide T and ruptured with buffer solution. Enriched poly(A) RNA of each sample was fragmented into 200–700 nt pieces with RNA Fragmentation Reagents. Then, the first-strand cDNA was synthesized with reverse transcriptase and arbitrary primers. Next, the second-strand cDNA was synthesized with DNA polymerase I and RNA enzyme. The double-stranded DNA was modified with Klenow fragment and T4 polynucleotide kinase successively. A base A was added on the 3'→5' exonuclease of the Klenow fragment. Then the base was connected to the corresponding adaptor using T4 ligase. Finally, fragments with around 200bp length were purified with QiaQuick GelPurify Kit (Qiagen, Hilden, Germany), and used as templates for PCR amplification to create the cDNA library. The library was paired-end sequenced using PE100 strategy on Illumina HiSeqTM 2500 (Illumina, San Diego, CA, USA) in Biomarker Technologies (Beijing, China).

Assembly and annotation of transcriptomes

Trinity (http://trinityrnaseq.sourceforge.net/) software was applied to perform the de novo assembly for the filtered high-quality data. The used parameters were as follow: min_glue = 2, V = 10, edge-thr = 0.05, min_kmer_cov = 2, kmer size = 25, path_reinforcement_distance = 80, and group_pairs_distance = 250. The other parameters were set as the default. The data samples were merged and assembled, and contigs were obtained through overlap of the assembled sequences. Then the contigs were clustered according to the paired-end information of sequences and similarity of contigs. Local assembly was conducted to generate transcripts. The longest transcript in each local region was selected for use as a unigene. The following parameters were used to ensure a high quality of assembly: a minimum of 95% identity, a minimum of 35 overlapping bases, a minimum of 35 scores and a maximum of 25 unmatched overhanging bases at sequence ends. The consensus cluster sequences and singletons make up the final unigene dataset. For functional annotations, we first searched all unigene sequences against various protein databases such as Nr, SwissProt, COG, and KEGG using BLASTX, and then searched nucleotide database Nt using BLASTN, with an E-value cut-off of 10⁻⁵. For inconsistencies between the unigene alignment results in different databases, priority was given to the nr, Swiss-Prot, KEGG, and COG data sequentially. Unigenes with no matched data in any of these databases were analyzed with ESTScan software to predict the coding region and direction of the sequence.

Analysis on the digital gene expression profile

According to data in the constructed cDNA library, differences in mRNA expression between samples #2 and #3 were identified using the digital gene expression (DGE) profile. Gene expression was quantified as reads per kb per million reads (RPKM) [21]. According on Audic et al.'s detection method for differentially expressed genes based on sequencing [22], a strict algorithm was adopted to screen the differentially expressed genes. Those with a false discovery rate (FDR) ≤0.001 and |fold change| ≥2 were considered to be differentially expressed. Combined with the function annotation of differentially expressed genes, the pathway (KEGG) and gene ontology (GO) enrichment analysis as well as the pattern clustering of differentially expressed genes were performed. Fisher's exact test was applied in the enrichment analysis. The results were adjusted using the Bonferroni correction method. Thus, the pathways with obviously enriched differentially expressed genes and GO functional categories were obtained for further analysis.

qRT-PCR veritification on DGE

D. citri larvae of the same age were collected and sprayed with the spore suspension of L. attenuatum at the concentration of 1×10⁸ spores/mL. The infection period lasted 5 days. A proper amount of D. citri was picked out every day, frozen with liquid nitrogen, and preserved at -80°C. D. citri larvae treated with sterile water were used as the control. The primers for differentially expressed genes and reference genes were designed using Primer Premier 5 software (S1 Table). The total RNA from the six treatments above was extracted separately, and 1 μg was used in the reverse transcription with HiScript II Q RT SuperMix for qPCR with gDNA wiper (Vazyme, China), followed by real-time fluorescent quantitative PCR with ChamQ SYBR qPCR Master Mix (Vazyme). The reaction was carried out on an ABI 7500 Fast Real-Time PCR System (Applied Biosystems, USA). Relative gene expression was calculated using the Pfaffl method [23]. Both D. citri α-tubulin and β-actin were used as reference genes. The analytic software for qPCR was 7500 softwear v2.0.6. SPSS Statistics v19.0.0 software was used to perform independent sample t-tests (P<0.05).

Results

Transcriptome sequencing and assembly

Through transcriptome sequencing (runs accession number: SRR2980521), we found a total of 24,071,734 reads. The mean CycleQ20 value of the samples reached 100.00%, and the base Q30 was 81.41%, suggesting that the sequencing was reliable. All high-quality reads were assembled de novo into 2,318,371 contigs (Table 1) with a mean length of 82 nt. Contigs were clustered according to the paired-end information of sequences and similarity of contigs and then assembled into 119,428 transcripts with the mean length of 952 nt. Among them, 31,748 transcripts (26.59%) had a length greater than 1,000 nt. These transcripts were finally assembled into 84,733 unigenes, with a mean length of 772 nt. Among them, 16,602 unigenes (19.6%) had a length greater than 1,000 nt. The N50 lengths of the contigs, transcripts, and unigenes were 98 nt, 1,809 nt and 1,320 nt, respectively. Regarding the length distribution of transcripts and unigenes of D. citri, most were 200–300 nt, followed by 300–500 nt, with the smallest group being sequences longer than 2,000 nt.

[Figure omitted. See PDF.]

Table 1. cDNA library assembly for D. citri.

https://doi.org/10.1371/journal.pone.0162659.t001

Function annotation of unigenes

Among the 84,733 assembled unigenes, 28,190 unigenes had been annotated, as found through the sequence alignment on the NCBI (National Center for Biotechnology Information database) website. Specifically, 11,516 unigenes were annotated in the COG database, 17,081 in the GO database, 19,537 in the Swiss-Prot database, 9,387 in the KEGG database, and 27,746 in the Nr database (Table 2).

[Figure omitted. See PDF.]

Table 2. Functional annotation of the D. citri transcriptome.

https://doi.org/10.1371/journal.pone.0162659.t002

It was observed in GO analysis that 17,081 unigenes were named as 148,314 GO terms, of which most genes included more than one GO term. GO analysis was used mainly to predict the functions of D. citri proteins. Generally, GO terms are classified into three categories (Fig 1): biological process, molecular function, and cellular component. In our analysis, most of the GO terms participated in biological process (90,520, accounting for 61.03% of the total), followed by molecular function (32,827, 22.13%), and cellular component (24,967, 16.83%). The three largest sub-categories of GO terms were cellular process (10,798 GO terms) and metabolic process (10,357 GO terms) in biological process and binding (8,959 GO terms) in the molecular function.

[Figure omitted. See PDF.]

Fig 1. Functional annotation of assembled sequences based on gene ontology (GO) categorization.

GO analysis was performed at level two for three main categories (cellular component, molecular function, and biological process).

https://doi.org/10.1371/journal.pone.0162659.g001

We also used COG classifications to analyze the putative protein functions. In total, 11,516 unigenes were functionally classified into 25 COG categories (Fig 2). The largest category was “General function prediction only” (3,976, 35%), followed by “Translation, ribosomal structure and biogenesis” (1,134, 9.9%), “Carbohydrate transport and metabolism” (1,097, 9.5%), “Replication, recombination and repair” (1,067, 9.3%), “Amino acid transport, metabolism” (1,057, 9.2%). There were only a few of unigenes taking part in “Nuclear structure” (7 unigenes), “Cell motility” (20 unigenes), “RNA processing and modification” (81 unigenes), and “Chromatin structure and dynamics” (94 unigenes). No unigenes took part in Extracellular structures.

[Figure omitted. See PDF.]

Fig 2. Histogram representation of Clusters of Orthologous Groups (COG) classification.

https://doi.org/10.1371/journal.pone.0162659.g002

There were 9,387 unigenes that had been annotated and classified in KEGG database, distributed over 208 pathways (Fig 3). The numbers of unigenes involved in different pathways varied. Those pathways with less than 100 unigenes were classified into one category (others), and the remaining 3,416 pieces of unigenes were classified into 19 categories. The ribosome pathway (ko03010) was the largest category with 456 pieces of unigenes. The pathways of protein processing in the endoplasmic reticulum (ko04141, 398 pieces), spliceosome (ko03040, 358 pieces), and RNA transportation (ko03013, 336 pieces) also involved many unigenes.

[Figure omitted. See PDF.]

Fig 3. Pie chart representation of the distribution of functional annotations of Kyoto Encyclopedia of Genes and Genomes (KEGG).

ko03010: Ribosome; ko04141: Protein processing in endoplasmic reticulum; ko03040: Spliceosome; ko03013: RNA transport; ko04142: Lysosome; ko00190: Oxidative phosphorylation; ko00230: Purine metabolism; ko04120: Ubiquitin mediated proteolysis; ko04144: Endocytosis; ko03015: mRNA surveillance pathway; ko04146: Peroxisome; ko03008: Ribosome biogenesis in eukaryotes; ko04310: Wnt signaling pathway; ko00500: Starch and sucrose metabolism; ko04145: Phagosome; ko00240: Pyrimidine metabolism; ko00010: Glycolysis / Gluconeogenesis; ko03018: RNA degradation; ko00020: Citrate cycle (TCA cycle).

https://doi.org/10.1371/journal.pone.0162659.g003

Identification of immunity-related genes

By searching the transcriptome, we preliminarily identified 129 pieces of immunity-related unigenes, involving pattern recognition receptors, signal transduction factors, response factors, and so on (Table 3). These genes play important roles in endogenous and exogenous immunoreactions of D. citri. For example, prophenoloxidase (PPO) can participate in melanization; lysozyme can be capable of dissolving cell walls of fungi and Gram-positive bacteria; and the antimicrobial peptide and clip domain serine protease also perform functions related to humoral immunity when a pathogen invasion occurs. In terms of the pattern recognition receptors, we totally identified a total of 49 unigenes including 1 PGRP, 6 GNBPs, 19 βGRPs, 2 TEPs, 8 SCRs, 3 CTLs, 5 GALEs, 3 hemocytins, and 2 integrins. These genes play an important role for the insect to recognize external pathogens, triggering the downstream reaction. PGRPs mainly identify the special ingredient in the bacterial cell walls—peptidoglycan—and then trigger the transcription of antibacterial peptide or PPO activation cascade. The protein was originally found in the silkworm, and until now, 12 PGRP genes have been identified [24]. However, in the transcriptome of D. citri, we identified only one PGRP gene (PGRP—S2). In dealing with the infection of gram negative bacteria and fungi, GNBP/βGRP can identify and combine with the beta-1,3-glucan and trigger PPO cascade. We performed phylogenetic analyses on the complete six genes, DcβGRP-5,7,8,16,17,19, in the sequence of D. citri, and found they have a close genetic relationship with AtβGRP-4a and MhGNBP-2 in Hemiptera (Fig 4).

[Figure omitted. See PDF.]

Table 3. Summary of the immunity-related unigenes identified in D. citri transcriptome.

https://doi.org/10.1371/journal.pone.0162659.t003

[Figure omitted. See PDF.]

Fig 4. Phylogenetic analysis of β-1,3-glucan recognition proteins (βGRPs) from D. citri and other insect species.

The used amino acid sequences are from Dialeurodes citri (Dc), Triatoma infestans (Ti), Anasa tristis (At), Bombyx mori (Bm), Maconellicoccus hirsutus (Mh), Anopheles gambiae (Ag), Tribolium castaneum (Tc), Tenebrio molitor (Tm), Ostrinia furnacalis (Of), Spodoptera litura (Sl), Manduca sexta (Ms), Galleria mellonella (Gm), Helicoverpa armigera (Ha), Culex quinquefasciatus (Cq).

https://doi.org/10.1371/journal.pone.0162659.g004

Four signal transduction pathways, Toll, Imd, JNK, and JAK/STAT are known to be involved ininsect immunity [25]. However, in D. citri, we only identified unigenes related to the Toll pathway, including 16 Tolls, 1 ECSIT, 6 Pelle, 1 pellino, and 3 NF-κB genes, and did not find genetic information in the other three pathways. We perform phylogenetic analysis on the Toll protein and found DcToll-7 has high similarity with other species (Fig 5). In the cascade mediated by serine protease, D. citri has 19 clip-domain serine proteases (clip-domain SP), 2 serine protease homologs (SPHs), and 6 serpins. Among them, DcSPH-1 and DcSPH-2 are prophenoloxidase-activating factors, which can catalyze conversion of PPO into PO (phenoloxidase). Serpin can negatively control the PPO activation, competing for PPO with SP.

[Figure omitted. See PDF.]

Fig 5. Phylogenetic analysis of Tolls from D. citri and other insect species.

The used amino acid sequences are from Dialeurodes citri (Dc), Nilaparvata lugens (Nl), Acyrthosiphon pisum (Ap), Graminella nigrifrons (Gn), Halyomorpha halys (Hh), Cimex lectularius (Cl), Lygus hesperus (Lh), Triatoma infestans (Ti), Drosophila melanogaster (Dm), Bombyx mori (Bm), Tribolium castaneum (Tc), Anopheles gambiae (Ag), Pediculus humanus corporis (Phc), Ceratitis capitata (Cc).

https://doi.org/10.1371/journal.pone.0162659.g005

As immune factors of D. citri, we identified 18 PPO unigenes, 6 lysozyme unigenes, and 3 AMPs. Most of the insects contained two kinds of PPOs, in the form of dimmers. Through clustering analysis, we found DcPPO-6 and DcPPO-7 with PPO of Hemiptera, and DcPPO1 and DcPPO-17 had a more distant relationship with other species (Fig 6). This may suggest that D. citri has more than two kinds of PPO genes. AMP includes cecropin, defensin, attacin, gloverin, and so on. D. citri has 3 AMP genes with sequences that are too short to permit identification of the type of AMP.

[Figure omitted. See PDF.]

Fig 6. Phylogenetic analysis of prophenoloxidases (PPOs) from D. citri and other insect species.

The used amino acid sequences are from Dialeurodes citri (Dc), Tenebrio molitor (Tm), Anasa tristis (At), Riptortus pedestris (Rp), Acyrthosiphon pisum (Ap), Drosophila melanogaster (Dm), Bombyx mori (Bm), Anopheles gambiae (Ag), Tribolium castaneum (Tc), Spodoptera litura (Sl), Sarcophaga bullata (Sb), Ostrinia furnacalis (Of), Manduca sexta (Ms), Holotrichia diomphalia (Hd), Hyphantria cunea (Hc), Helicoverpa armigera (Ha), Galleria mellonella (Gm), Anopheles stephensi (Ans), Armigeres subalbatus (Ars), Nilaparvata lugens (Nl).

https://doi.org/10.1371/journal.pone.0162659.g006

Expression profile analysis

We obtained the RPKM values for unigenes, and through further calculation of Log₂ fold changes [Log₂Ratio(treatment group RPKM/control group RPKM)], we found 441 differentially expressed unigenes in control group and treatment group, among which 313 unigenes had been annotated (Data Accessibility) in the transcriptome database. The clustering analysis of differentially expressed genes suggested that Log₂FC values of most genes in the treatment group were positive (Fig 7A). Namely, these genes were upregulated in comparison to the control group. It was found in statistical analysis of differentially expressed genes that the genes with Log₂FC values in the range of -2 to 2 were in the majority (210 pieces, 47.62%; Fig 7B). Genes with Log₂Ratio values in the range of 2 to 4, ranked second with 148 pieces (33.56%). Generally, only small numbers of differentially expressed genes were downregulated. There were 38 genes (8.62%) with Log₂FC values of -4 to -2, and 9 genes (2.04%) with Log₂FC values of -6 to -4. In our two experimental groups, only one gene was upregulated by more than 10-fold, and its gene ID was unigene78104_c1. The gene ID of the most downregulated gene (by -6.38-fold) was unigene80382_c2.

[Figure omitted. See PDF.]

Fig 7. Clustering analysis (A) and quantitative analysis (B) of differentially expressed genes in D. citri after treatment with L. attenuatum.

The numbers beside the bar mean the values of Log₂FC.

https://doi.org/10.1371/journal.pone.0162659.g007

We analyzed 24 genes that may participate in the response against infection by L. attenuatum (Table 4). Among them, function annotation identified that the cuticle protein was the highest upregulated gene with a Log₂Ratio value of 10.04, followed by the vitellogenin gene. This implies that the two proteins may play roles in the response of D. citri against infection by L. attenuatum directly or indirectly. Regarding the melanization, we found three PPO genes and three clip-domain SP genes. Among them, the expression difference of DcSP-3 was 2.51-fold, and those of the two PPO genes were 1–2-fold. In addition, we found two lysozyme unigenes with fold changes of 2.47 and 2.60. Regarding the metabolic reaction process of insects, we examined the expression of cathepsin B.

[Figure omitted. See PDF.]

Table 4. Differentially expressed unigenes likely involved in the antifungal response of D. citri after treatment with L. attenuatum.

https://doi.org/10.1371/journal.pone.0162659.t004

Pathway analysis of upregulated genes of interest suggested that the differentially expressed genes were mainly involved in two pathways, the lysosome (ko04142) and MAPK signaling pathways (ko04010). In the lysosome pathway, 9 pieces of unigenes exhibited differential expression, with 8 upregulated and 1 downregulated. All of these genes are crucial to lysosome phagocytosis and immunoreactions. The MAPK signaling pathway plays the role of signal transduction in the processes of stress adaptation and inflammatory response. For this pathway, there were two pieces of differentially expressed genes, both of which were related to heat shock cognate 70.

qRT-PCR analysis of differentially expressed genes

To verify the accuracy of fold changes in the digital gene expression profile and investigate the dynamic tendency of key gene expression with infection by L. attenuatum, we used qRT-PCR to analyze key genes in the response of D. citri to infection for 1–5 days.

qRT-PCR analysis was conducted for the two pieces of cuticle protein unigenes most upregulated. After 1 day of infection, no significant change (P<0.05) in their expression was observed compared with that in the control group (CG) treated with clear water. However, by the second day of infection, DcCP-1 exhibited increased expression (Fig 8). The relative expression change after 2 days was 2.20-fold, and that after 3 days was 3.60-fold. Unexpectedly, the gene expression increased sharply on the fourth day, with 267.83-fold upregulation and 865.54-fold upregulation on the fifth day. The gene DcCP-2 exhibited obvious downregulation (P<0.05) on the second day and rebounded on the third day. On the fourth day, expression of the gene increased sharply as for DcCP-1, to 92.87-fold, and further increased to 202.63-fold on the fifth day.

[Figure omitted. See PDF.]

Fig 8. qRT-PCR–based verification and analysis of the differential expression of the identified genes.

Data are shown as means of three replicates ± standard deviation (SD). Asterisks above bars indicate significance differences between the treatment and control groups (t-test, P<0.05). The expression levels of unigenes in the control group (CG) are marked with a dashed line at Y = 1.0. The abscissa values 1, 2, 3, 4, 5 mean the treatment times (day) of D.citri infected by L. attenuatum.

https://doi.org/10.1371/journal.pone.0162659.g008

We verified the two upregulated vitellogenin genes, DcVg-1 and DcVg-2, and found that they showed the same expression tendency (Fig 8). During the first 2 days of infection, their expression was inhibited compared to that in the control group (P<0.05). On the third day, their expression levels were comparable to those in the control group, and after infection for 4 days, obvious upregulation of the two genes occurred (215.90-fold and 199.98-fold respectively). On the fifth day, these genes were upregulated by 327.19-fold and 287.84-fold, respectively. Thus, expression of the two vitellogenin genes was positively correlated with the infection time.

Among the genes participating in the lysosome pathway, two cathepsin B genes, DcCatB-1 and DcCatB-2, were analyzed through qPCR. DcCatB-1 was upregulated significantly (P<0.05) by 2.17-fold by the first day of infection and continued to be upregulated on days 2–5, with a peak fold change value of 3.69-fold on the second day (Fig 8). Expression of DcCatB-2 was inhibited on the first day, with a 0.43-fold difference from that in the control group. However, it was upregulated obviously on the second day and reached maximum upregulation by 3.93-fold on the third day. On the fourth and fifth days, its expression decreased, which might be caused by depletion of the enzyme in the latter stage.

For the melanization pathway, we tested the expression levels of DcPPO-7 and DcSP-3 (Fig 8). The results indicated that these two genes were significantly upregulated (P<0.05) on the first day of infection by 4.40-fold and 4.12-fold, respectively, and by 40.88-fold and 35.71-fold, respectively, by the fifth day. During melanization, serine protease can activate PPO and convert it to PO, which participates in the generation of melanin. This is an immediate immune response to the invasion of fungi. The two genes exhibited obvious upregulation on the first day of infection, which supports this theory.

We measured the expression level of the lysozyme gene DcLys-4 (Fig 8) and found that, like PPO, it was significantly upregulated (P<0.05) by 2.81-fold on the first day of infection, 17.37-fold on the third day, and 27.00-fold on the fifth day.

Heat shock cognate 70 (DcH70-1, Fig 8) acts in the MAPK signaling pathway and responds to exogenous stimuli. Expression of this gene did not change significantly in the first 2 days of infection. But during the last 3 days, it was significantly upregulated (P<0.05) with a 7.34-fold change on the fifth day.

In the expression profile, the hypothetical protein gene DcHP-1 also showed a large fold change of 8.13-fold on the fourth day and by 22.82-fold (Fig 8) on the fifth day. Further investigation of the importance of this gene is warranted.

Discussion

We obtained a transcriptome database containing 84,733 unigenes of D. citri using RNA-Seq technology. This amount of data was greater than that obtained previously [26]. By searching and screening the transcriptome database, we identified 129 immunity-related unigenes. These unigenes are related to pattern recognition receptors, information transduction factors and response factors. Using the digital gene expression profile, we identified 441 differentially expressed genes in D. citri infected with L. attenuatum. Among these genes, we found the response factors such as the PPO, lysozyme, and clip-domain SP participated in cascade. However, we did not find genes related to the pattern recognition receptors, suggesting that as the final effector, response factors show fluctuating expression

The cuticle is the first barrier for insects to defend against pathogen infection, in addition to being indispensible for maintaining the shape and mobility of insects [27, 28]. The major components of the insect cuticle are chitin and cuticle protein [29]. Cuticle protein contributes much to the stress resistance, drug resistance, and immunity of insects. When an insect is suffering from adverse environment conditions, cuticle protein genes are induced to strengthen or stabilize the cuticular structure, resist the effects of adverse factors, and maintain the insect's survival [30–32]. In researching the aphid's insecticide resistance mechanism, Silva et al. found that two RR2-type cuticle protein genes are upregulated [33]. Asano et al. found that when B. mori larvae are subjected to bacterial infection, the cuticle protein gene BmCb10 is significantly upregulated. They speculated that the gene could transmit the exogenous adverse stimulation to activate melanization [34]. He et al. proposed that cuticle protein may play a role in wound healing in Anopheles gambiae larvae and adults [35]. In this study, two cuticle protein genes showed significant upregulation upon D. citri infection by L. attenuatum. This finding provides some evidence for the role of cuticle protein in the immune defense of D. citri.

Melanization, an important aspect of the insect defense system, involves the regulation of the melanin cascade mediated by PPO [36]. Upon pathogen invasion, it activates PPO and transforms it into PO, which can transform phenolic substances into quinone intermediates, that aggregate and form melanin before enclosing, isolating, and killing pathogens. In addition, PO takes part in the processes of wound healing and skin hardening. The activation and transformation of PPO into PO is generally considered to be completed through the cascade of clip-domain SPs [37]. Gillespie et al. found that when Schistocerca gregaria was infected by Metarhizium anisoplia, the PPO level in the body increased while PO activity decreased, and the lysozyme level exhibited a significant decrease in comparison with the control group [11]. By knocking out the PPO gene in the mosquito Armigeres subalbatus, it was found that the melanization function was influenced greatly [38]. When fungal spores were injected into the insect body, the PO expression level was increased significantly [39]. We measured the expression levels of clip-domain SP and PPO in D. citri infected by L. attenuatum. On the first day of infection, both were significantly upregulated, and the fold changes in expression exceeded 35-fold by the fifth day. This indicates that these enzymes contribute much to the response against exogenous pathogenic fungi.

In the Digital Gene Expression Profiling, vitellogenin and cathepsin B with high expression levels received our attention. vitellogenin is an important source of energy for D. citri, and cathepsin B is an indispensable metabolic enzyme of lysosome. However, research on their roles in the immune response remains incomplete. Guo et al. showed that when Bemisia tabaci were fed a virus-infected plant, vitellogenin expression increased significantly [40]. Soderhall et al. demonstrated that a clotting protein belonging to the vitellogenin superfamily participates in the autologous immune defense in freshwater crayfish [41]. Raikhel et al. found that the upstream regulatory region of the vitellogenin gene of Aedes aegypti participates in the immune defense against pathogens [42]. Shi et al. demonstrated that vitellogenin can agglutinate erythrocytes of toad and chicken and has an inhibitory effect on various bacteria [43]. Amdma et al. showed that vitellogenin participates in the regulation of immune function and life in bees [44]. Cathepsin B enzymes can part in the processes of immune evasion [45, 46]. Futahashi et al. found upregulated expression of the cathepsin gene in many tissues of Burkholderia-infected Riptortus pedestris [47]. Zhang et al. cloned the cathepsin O gene of B. mori and detected its expression in the hemolymph of B. mori treated with Escherichia coli, they speculated that the enzyme participates in the immune response in the body [48]. Kocks et al. found that after immune stimulation, cysteine protease L was activated and highly expressed in lysosomes [49]. Wang et al. investigated two Cathepsin L genes of channel catfish and proved that the two genes acted in mucosal immunity [50].

In future studies, we may take advantage of the RNAi technique to verify the effect of these genes in the response of D. citri to infection of entomopathogenic fungi, and improve the fungal toxicity and field application effectiveness.

Supporting Information

[Figure omitted. See PDF.]

S1 Fig. Total RNA extraction from D. citri.

1: Total RNA for transcriptome sequencing; 2,3: Total RNA for digital gene expression profiling.

https://doi.org/10.1371/journal.pone.0162659.s001

(TIF)

S1 Table. Primer pairs for real time quantitative PCR in D. itri.

https://doi.org/10.1371/journal.pone.0162659.s002

(DOC)

Acknowledgments

This work was supported by the National Science and Technology Support Program (2012BAD19B06), the National Sparking Plan Project (2012GA811001), and Technological Projects of Chongqing (cstc2014yykfB80007, cstc2014fazktjcsf80033).

Author Contributions

1. Conceptualization: CR SY.

2. Formal analysis: SY LD.

3. Funding acquisition: CR LC HL.

4. Investigation: JY.

5. Methodology: SY XL.

6. Software: RL.

7. Validation: SY.

8. Writing – original draft: SY.

9. Writing – review & editing: SY RL.

Citation: Yu S, Ding L, Luo R, Li X, Yang J, Liu H, et al. (2016) Identification of Immunity-Related Genes in Dialeurodes citri against Entomopathogenic Fungus Lecanicillium attenuatum by RNA-Seq Analysis. PLoS ONE 11(9): e0162659. https://doi.org/10.1371/journal.pone.0162659

References

1. Mound LA, Halsey SH. Whitefly of the world. A systematic catalogue of the Aleyrodidae (Homoptera) with host plant and natural enemy data. John Wiley and Sons; 1978.

2. Nguyen R, Sailer RI, Hamon AB. Catalog of Aleyrodidae on Citrus and their natural enemies (Homoptera-Aleyrodidae). Contribution/Bureau of Entomology (USA); 1993.

3. Fasulo TR, Weems HV. Citrus Whitefly, Dialeurodes citri (Ashmead) (Insecta: Hemiptera: Aleyrodidae). UF/IFAS Extension. 1999; 128: 1–6.

4. Gill G. Citrus whitefly blows in. N Z Biosecurity. 2001; 26: 19.

5. Pyle K. The current status and economic impact of Australian citrus whitefly on New Zealand citrus. Report to New Zealand Citrus Growers Inc., Pyle Orchards and Consulting, Katikati, Tauranga, New Zealand; 2009. pp. 27.

6. Ma W, Li X, Dennehy TJ, Lei C, Wang M, Degain BA, et al. Pyriproxyfen resistance of Bemisia tabaci (Homoptera: Aleyrodidae) biotype B: metabolic mechanism. J Econ Entomol. 2010; 103: 158–165. pmid:20214381

7. Karatolos N, Denholm I, Williamson M, Nauen R, Gorman K. Incidence and characterisation of resistance to neonicotinoid insecticides and pymetrozine in the greenhouse whitefly, Trialeurodes vaporariorum Westwood (Hemiptera: Aleyrodidae). Pest Manag Sci. 2010; 66: 1304–1307. pmid:20799247

8. Karatolos N, Williamson MS, Denholm I, Gorman K, Ffrench-Constant RH, Bass C. Over-expression of a cytochrome P450 is associated with resistance to pyriproxyfen in the greenhouse whitefly Trialeurodes vaporariorum. Plos One. 2012; 7: e31077. pmid:22347432

9. Hall RA. Control of whitefly, Trialeurodes vaporariorum and cotton aphid, Aphis gossypii in glasshouses by two isolates of the fungus, Verticillium lecanii. Ann Appl Biol. 1982; 101: 1–11.

10. Korolev N, Gindin G. Vegetative compatibility in the entomopathogen Verticillium lecanii. Mycol Res. 1999; 103: 833–840.

11. Gillespie JP, Bailey AM, Cobb B, Vilcinskas A. Fungi as elicitors of insect immune responses. Arch Insect Biochem. 2000; 44: 49–68.

12. Sultana K, Batra LR, Stavely JR, Nasir MA. Hyperparasitism of Verticillium lecanii and Cladosporium cladosporioides on Uromyces appendiculatus, the causal organism for soybean rust. Pakistan J Phytopathol. 2000; 12: 42–45.

13. Kim JJ, Goettel MS, Gillespie DR. Potential of Lecanicillium species for dual microbial control of aphids and the cucumber powdery mildew fungus, Sphaerotheca fuliginea. Biol Control. 2007; 40: 327–332.

14. Zare R, Gams W. A revision of Verticillium section Prostrata. IV. The genera Lecanicillium and Simplicillium gen. nov.°). Nova Hedwigia. 2001; 73: 1–50.

15. Kim JJ. Influence of Lecanicillium attenuatum on the development and reproduction of the cotton aphid, Aphis gossypii. Biocontrol. 2007; 52: 789–799.

16. Gao Y, Xie YP, Xiong Q, Liu WM, Xue JL. Ultrastructural Exploration on the Histopathological Change in Phenacoccus fraxinus Infected with Lecanicillium lecanii. Plos One. 2015; 10: e0117428. pmid:25629309

17. Kim JJ, Kim KC. Selection of a highly virulent isolate of Lecanicillium attenuatum against cotton aphid. J Asia-Pac Entomol. 2008; 11: 1–4.

18. Hillyer JF, Christensen BM. Characterization of hemocytes from the yellow fever mosquito, Aedes aegypti. Histochem Cell Biol. 2002; 5: 431–440.

19. Söderhäll K, Cerenius L. Role of the prophenoloxidase-activating system in invertebrate immunity. Curr Opin Immunol. 1998; 1: 23–28.

20. Imler JL, Hoffmann JA. Signaling mechanisms in the antimicrobial host defense of Drosophila. Curr Opin Microbiol. 2000; 3: 16–22. pmid:10679426

21. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008; 5: 621–628. pmid:18516045

22. Audic S, Claverie JM. The significance of digital gene expression profiles. Genome Res. 1997; 7: 986–995. pmid:9331369

23. Pfaffl MW. A new mathematical model for relative quantification in real-time RT–PCR. Nucleic Acids Res. 2001; 29: e45–e45. pmid:11328886

24. Tanaka H, Ishibashi J, Fujita K, Nakajima Y, Sagisaka A, Tomimoto K, Sunagawa T. A genome-wide analysis of genes and gene families involved in innate immunity of Bombyx mori. Insect Biochem Molec. 2008; 38: 1087–1110.

25. Evans JD, Aronstein K, Chen YP, Hetru C, Imler JL, Jiang H, et al. Immune pathways and defence mechanisms in honey bees Apis mellifera. Insect Mol Biol. 2006; 15: 645–656. pmid:17069638

26. Chen EH, Wei DD, Shen GM, Yuan GR, Bai PP, Wang JJ. De novo characterization of the Dialeurodes citri transcriptome: mining genes involved in stress resistance and simple sequence repeats (SSRs) discovery. Insect Mol Biol. 2014; 23: 52–66. pmid:24164346

27. Delon I, Payre F. Evolution of larval morphology in flies: get in shape with shavenbaby. Trends Genet. 2004; 20: 305–313. pmid:15219395

28. Moussian B, Schwarz H, Bartoszewski S, Nüsslein-Volhard C. Involvement of chitin in exoskeleton morphogenesis in Drosophila melanogaster. J Morphol. 2005; 264: 117–130. pmid:15747378

29. Moussian B. Recent advances in understanding mechanisms of insect cuticle differentiation. Insect Biochem Molec. 2010; 40: 363–375.

30. Zhang J, Goyer C, Pelletier Y. Environmental stresses induce the expression of putative glycine-rich insect cuticular protein genes in adult Leptinotarsa decemlineata (Say). Insect Mol Biol. 2008; 17: 209–216. pmid:18477239

31. Purać J, Burns G, Thorne MA, Grubor-Lajšić G, Worland MR, Clark MS. Cold hardening processes in the Antarctic springtail, Cryptopygus antarcticus: clues from a microarray. J Insect Physiol. 2008; 54: 1356–1362. pmid:18703067

32. Teets NM, Denlinger DL. Surviving in a frozen desert: environmental stress physiology of terrestrial Antarctic arthropods. J Exp Biol. 2014; 217: 84–93. pmid:24353207

33. Silva AX, Jander G, Samaniego H, Ramse JS, Figueroa CC. Insecticide resistance mechanisms in the green peach aphid Myzus persicae (Hemiptera: Aphididae) I: a transcriptomic survey. Plos One. 2012; 7: e36366. pmid:22685538

34. Asano T, Taoka M, Shinkawa T, Yamauchi Y, Isobe T, Sato D. Identification of a cuticle protein with unique repeated motifs in the silkworm, Bombyx mori. Insect Biochem Molec. 2013; 43: 344–351.

35. He N, Botelho JM, McNall RJ, Belozerov V, Dunn WA, Mize T, et al. Proteomic analysis of cast cuticles from Anopheles gambiae by tandem mass spectrometry. Insect Biochem Molec. 2007; 37: 135–146.

36. Taft AS, Chen CC, Li J, Christensen BM. Molecular cloning of two prophenoloxidase genes from the mosquito Aedes aegypti. Insect Mol Biol. 2001; 10: 97–103. pmid:11240641

37. Gillespie JP, Kanost MR, Trenczek T. Biological mediators of insect immunity. Annu Rev Entomol. 1997; 42: 611–643. pmid:9017902

38. Ahmed A, Martin D, Manetti AGO, Han SJ, Lee WJ, Mathiopoulos KD et al. Genomic structure and ecdysone regulation of the prophenoloxidase 1 gene in the malaria vector Anopheles gambiae. P Natl Acad Sci USA. 1999; 96: 14795–14800.

39. Söderhäll K, Cerenius L. Role of the prophenoloxidase-activating system in invertebrate immunity. Curr Opin Immunol. 1998; 10: 23–28. pmid:9523106

40. Guo JY, Dong SZ, Yang XL, Cheng L, Wan FH, Liu SS et al. Enhanced vitellogenesis in a whitefly via feeding on a begomovirus-infected plant. Plos One. 2012; 7: e43567. pmid:22937062

41. Söderhäll K. Defence reactions in a crustacean. Dev Comp Immunol. 1997; 2: 137.

42. Raikhel AS, Kokoza VA, Zhu J, Martin D, Wang SF, Li C et al. Molecular biology of mosquito vitellogenesis: from basic studies to genetic engineering of antipathogen immunity. Insect Biochem Molec. 2002; 32: 1275–1286.

43. Shi X, Zhang S, Pang Q. Vitellogenin is a novel player in defense reactions. Fish Shellfish Immun. 2006; 20: 769–772.

44. Amdam GV, Simões ZL, Hagen A, Norberg K, Schrøder K, Mikkelsen Ø et al. Hormonal control of the yolk precursor vitellogenin regulates immune function and longevity in honeybees. Exp Gerontol. 2004; 39: 767–773. pmid:15130671

45. Ford L, Zhang J, Liu J, Hashmi S, Fuhrman JA, Oksov Y et al. Functional analysis of the cathepsin-like cysteine protease genes in adult Brugia malayi using RNA interference. Plos Neglect Trop D. 2009; 3: e377.

46. Rispe C, Kutsukake M, Doublet V, Hudaverdian S, Legeai F, Simon JC et al. Large gene family expansion and variable selective pressures for cathepsin B in aphids. Mol Biol Evol. 2008; 25: 5–17. pmid:17934209

47. Futahashi R, Tanaka K, Tanahashi M, Nikoh N, Kikuchi Y, Lee BL et al. Gene expression in gut symbiotic organ of stinkbug affected by extracellular bacterial symbiont. Plos One. 2013; 8: e64557. pmid:23691247

48. Zhang K, Su J, Chen S, Yu S, Tan J, Xu Met al. Molecular cloning, characterization and expression analysis of cathepsin O in silkworm Bombyx mori related to bacterial response. Mol Immunol. 2015; 66: 409–417. pmid:25996894

49. Kocks C, Maehr R, Overkleeft HS, Wang EW, Iyer LK, Lennon-Dumenil AM et al. Functional proteomics of the active cysteine protease content in Drosophila S2 cells. Mol Cell Proteomics. 2003; 2: 1188–1197. pmid:13130081

50. Wang R, Song L, Su B, Zhao H, Zhang D, Peatman E et al. Mucosal expression signatures of two Cathepsin L in channel catfish (Ictalurus punctatus) following bacterial challenge. Fish & shellfish immun. 2015; 47: 582–589.