Journal of Invertebrate Pathology 129 (2015) 1–6

Contents lists available at ScienceDirect

Journal of Invertebrate Pathology journal homepage: www.elsevier.com/locate/jip

Homologs to Cry toxin receptor genes in a de novo transcriptome and their altered expression in resistant Spodoptera litura larvae Liang Gong a,b, Huidong Wang b, Jiangwei Qi b, Lanzhi Han c, Meiying Hu b,⇑, Juan Luis Jurat-Fuentes d,⇑ a Key Laboratory of Plant Resource Conservation and Sustainable Utilization, South China Botanical Garden, Chinese Academy of Sciences, Guangzhou 510650, Guangdong Province, China b Key Laboratory of Pesticide and Chemical Biology, College of Natural Resources and Environment, South China Agricultural University, Guangzhou 510642, Guangdong Province, China c State Key Laboratory for Biology of Plant Diseases and Insect Pests, Institute of Plant Protection, Chinese Academy of Agricultural Sciences, Beijing 100193, China d Department of Entomology and Plant Pathology, University of Tennessee, Knoxville, TN 37996, USA

a r t i c l e

i n f o

Article history: Received 6 November 2014 Revised 4 May 2015 Accepted 8 May 2015 Available online 14 May 2015 Keywords: Spodoptera litura Bacillus thuringiensis Transcriptome Cry1Ca receptor Resistance

a b s t r a c t Insect resistance threatens sustainability of insecticides based on Cry proteins from the bacterium Bacillus thuringiensis (Bt). Since high levels of resistance to Cry proteins involve alterations in Cry-binding midgut receptors, their identification is needed to develop resistance management strategies. Through Illumina sequencing we generated a transcriptome containing 16,161 annotated unigenes for the Oriental leafworm (Spodoptera litura). Transcriptome mining identified 6 contigs with identity to reported lepidopteran Cry toxin receptors. Using PCR we confirmed their expression during the larval stage and compared their quantitative expression in larvae from susceptible and a field-derived Cry1Ca resistant strain of S. litura. Among reduced transcript levels detected for most tested contigs in the Cry1Ca-resistant S. litura larvae, the most dramatic reduction (up to 99%) was detected for alkaline phosphatase contigs. This study significantly expands S. litura transcriptomic resources and provides preliminary identification of putative receptor genes with altered expression in S. litura resistant to Cry1Ca toxin. Ó 2015 Elsevier Inc. All rights reserved.

1. Introduction Insecticidal Cry proteins from the bacterium Bacillus thuringiensis (Bt) are used in pesticides and produced by transgenic crops against agricultural pests and insect vectors of human disease (Sanchis, 2011), yet evolution of insect resistance threatens the sustainability of these technologies. While populations of three lepidopteran pests have evolved practical resistance to Bt crops (Dhurua and Gujar, 2011; Storer et al., 2010; van Rensburg, 2007), the mechanisms involved are unknown. Diverse research

Abbreviations: Bt, Bacillus thuringiensis; APN, aminopeptidase-N; GPI, glycosylphosphatidyl-inositol; ALP, alkaline phosphatase; ABCC2, ATP-binding cassette transporter subfamily C2; GO, gene ontology; COG, cluster of orthologous groups of proteins; KO, KEGG orthology; CR, cadherin repeats; BBMV, brush border membrane vesicles; NGS, next-generation sequencing; RPKM, reads per kilobase per million mapped reads. ⇑ Corresponding authors at: Tianhe No. 483, Guangzhou 510642, Guangdong Province, China. Tel.: +86 (020) 85280308 (M. Hu). 370 Plant Biotechnology Building, 2505 E.J. Chapman Drive, Knoxville, TN 37996, USA. Tel.: +1 (865) 974 5931 (J.L. Jurat-Fuentes). E-mail addresses: [email protected] (M. Hu), [email protected] (J.L. Jurat-Fuentes). http://dx.doi.org/10.1016/j.jip.2015.05.008 0022-2011/Ó 2015 Elsevier Inc. All rights reserved.

has identified binding of Cry toxin to receptor proteins in the insect midgut as critical for susceptibility, with practical resistance to Bt crops and pesticides associated with alterations in toxin-receptor interactions through mutations or down-regulation of receptor genes (Adang et al., 2014). Consequently, the identification of Cry toxin receptor genes is critical to developing effective resistance management strategies. Traditional biochemical and genetic analyses have identified diverse putative Cry receptors in the midgut of lepidopteran larvae, including cadherins, glycosylphosphatidylinositol (GPI)–anchored aminopeptidase-N (APN) and alkaline phosphatase (ALP), and ATP-binding cassette transporter subfamily C2 (ABCC2) proteins (Heckel, 2012; Pigott and Ellar, 2007). Recent advances in nucleic acid sequencing have allowed the identification of genes with homology to putative Cry toxin receptors in non-model insects (Pauchet et al., 2009) and testing their alterations in resistant insects (Lei et al., 2014; Tetreau et al., 2012). Larvae of the Oriental leafworm moth, Spodoptera litura (Fabricius) (Lepidoptera: Noctuidae), are devastating pests of diverse crops in the Indo-Australian tropics and are susceptible to selected Cry proteins, most notably Cry1Ca (van Frankenhuyzen, 2009). The only reported Cry toxin receptor in S. litura larvae is a

2

L. Gong et al. / Journal of Invertebrate Pathology 129 (2015) 1–6

midgut APN (SlAPN) that binds (Agrawal et al., 2002; Kaur et al., 2007) and serves as a functional receptor for Cry1Ca toxin (Rajagopal et al., 2002). The APN epitope interacting with SlAPN has been recently identified and reported using phage display (Kaur et al., 2014). While there are no data available on mechanisms of resistance to Cry1Ca in S. litura, resistance to Cry1Ca in the closely-related species Spodoptera exigua involved lack of expression of an APN1 gene (Herrero et al., 2005) with low identity (34%) to SlAPN. In contrast, Spodoptera frugiperda resistance to Cry1Fa was associated with reduced ALP levels (Jurat-Fuentes et al., 2011). The goal of the present study was to use Illumina sequencing to produce a S. litura transcriptome that would assist in identifying Cry1 toxin receptor homologs and their alterations in resistant populations. We identified 6 contigs expressed in larvae that encoded homologs to Cry toxin receptors, including cadherin, APN, ALP and ABCC2 genes. We also report the isolation of field-collected Cry1Ca-resistant S. litura and the comparative expression of Cry receptor homologs in susceptible and Cry1Ca-resistant S. litura larvae. 2. Materials and methods 2.1. Insect strains and bioassays A laboratory (LS) S. litura strain reared in the laboratory for >10 years and a field-derived population (FS) were maintained as previously reported (Gong et al., 2014). The FS strain was generated from eggs collected at vegetable fields in Guangzhou Tianhe district in China (23.18°N, 113.39°E) in 2012. Purified activated Cry1Ca toxin was provided by Marianne Pusztai-Carey (Case Western Reserve University, Cleveland, OH, USA). Neonates from the LS and FS strains were used in bioassays with Cry1Ca diluted with distilled water containing 0.1% Triton X-100 on the surface of meridic diet (Hinks and Byers, 1976). Ten insects and three biological replicates were used per toxin dose. Mortality was assessed 7 days after neonate inoculation and mortality parameters were calculated by probit analysis (Finney, 1971). 2.2. RNA isolation and cDNA library construction Samples from eggs, larvae (1st–4th instar), pupae and adults (male and female) of the LS strain were collected during 24 h after metamorphosis or molt for RNA isolation. Each sample (80 mg) was kept at 80 °C until used for total RNA extraction with Trizol reagent (Life Technologies, CA, USA) according to manufacturer’s instructions. Genomic DNA was eliminated by treatment with DNAse I (New England Biolabs, MA, USA). Total RNA concentration was measured using the QubitÒ RNA Assay Kit in a QubitÒ 2.0 Flurometer (Life Technologies), and integrity assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). All samples had RIN (RNA Integrated Number) values >8. Equal amounts (0.75 g) of RNA from each developmental stage were pooled, and the resulting RNA (3 g) was used for transcriptome sequencing. Sequencing libraries were generated using the Illumina TruSeq™ RNA Sample Preparation Kit (Illumina, San Diego, USA) following manufacturer’s recommendations. 2.3. Sequencing and de novo transcriptome assembly Clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer’s instructions. After

cluster generation, libraries were sequenced on an Illumina Hiseq 2000 platform to generate 100 bp paired-end reads. Raw reads in fastq format were processed through custom Perl scripts to remove adapter sequences and reads containing poly-N and of low quality, resulting in clean reads. For quality assessment of the sequencing, Q20 (Qphred = 10log 10(e), where e is the error rate of sequencing) and GC-content values were calculated. Unigenes were assembled from the cleaned sequence data using Trinity (Grabherr et al., 2011), and all the obtained contigs were used in the downstream analysis. The raw data generated from this study were deposited in the Sequence Read Archive (SRA) at NCBI with accession No. SRX469308. 2.4. Sequencing analysis The RSEM software (Li and Dewey, 2011) was used to estimate the number of reads mapping to each unigene as RPKM (Reads Per Kilobase per Million mapped reads). All the unigenes were used for sequence similarity searches of the NCBI nr (non-redundant) protein and nucleotide databases, and the Swiss-Prot and Pfam (http://pfam.sanger.ac.uk/) databases using BLAST with an E-value cutoff of 10 5 (Korf, 2003). Functional annotation by gene ontology (GO) was analyzed by Blast2GO (Conesa et al., 2005). The COG and KEGG pathway annotations were performed using Blastall software by searching the Cluster of Orthologous Groups (http:// www.ncbi.nlm.nih.gov/COG/) and Kyoto Encyclopedia of Genes and Genomes databases (http://www.genome.jp/kegg/), respectively. The transcriptome shotgun assembly project was deposited at GenBank under accession GBBY00000000. The version described in this manuscript is the first version, GBBY01000000. 2.5. Semi-quantitative and quantitative PCR Expression of selected genes during the egg, larva (pooled 1st, 2nd, 3th and 4th instar larvae), pupa and adult (pooled females and males) stages of S. litura was determined by semi-quantitative PCR. Total RNA extracted as above was used for first-strand cDNA synthesis using an oligo (dT)18 primer and AMV reverse transcriptase (TaKaRa). Semi-quantitative PCR reactions were carried out using cDNA (600 ng), 0.4 lM of each primer (Table S1), and LA Taq (Takara Dalian) in a final reaction volume of 25 ll. The S. litura b-actin gene (GenBank number DQ494753) was used as internal reference for normalization. The PCR amplification was performed using the following conditions: one cycle: (94 °C, 3 min); 27 cycles (94 °C, 30 s; 50 °C, 45 s; 72 °C, 1 min) and a last cycle 72 °C for 10 min. The presence of amplicons was tested using 1.5% agarose gel electrophoresis and directly sequenced. Quantitative RT-PCR (qRT-PCR) was performed in midguts (three 5th instar larval guts per biological replicate) from larvae of the LS and FS strains. Reactions were carried out on a BioRad iQ5 real-time PCR detection system using 100 ng of cDNA, 0.2 lM of primers (Table S1) and SYBR Premix Ex Taq (TaKaRa). Amplification consisted of an initial denaturation step at 95 °C for 30 s followed by 40 cycles of 95 °C for 5 s, 55 °C for 30 s and a dissociation step added at the end. Analysis of the amplification and melting curves was performed according to the manufacturer’s instructions. The relative amounts of transcript were first normalized to the endogenous reference gene (same as for semi-quantitative PCR), and then normalized relative to the transcript levels in the laboratory (LS) strain of S. litura according to the 2 DDCt method (Livak and Schmittgen, 2001). Data shown are the means and corresponding standard errors calculated from three biological replicates each tested in triplicate. Significance of expression differences was tested with an ANOVA test on rank-transformed data followed by a post hoc multiple comparison versus a control (expression in the susceptible LS strain) procedure

L. Gong et al. / Journal of Invertebrate Pathology 129 (2015) 1–6

3

(Dunett’s Method) with significance at the p < 0.01 level in Sigma Plot v.11.0 (Systat Software Inc.). 3. Results and discussion 3.1. Rationale, sequencing and functional annotation Selected Bt pesticides (Chatterjee and Mondal, 2012) and transgenic crops expressing Bt toxin genes (Zaidi et al., 2009) have been proposed for effective control of S. litura larvae. However, there is no information on potential resistance mechanisms in S. litura against Bt toxins, which is necessary for developing effective resistance management practices. This knowledge gap is further exacerbated by the unavailability of significant genetic data for S. litura in public databases. Toward resolving this issue, we performed Illumina sequencing of the S. litura transcriptome to generate 53,063,813 raw reads, which after removal of adapters (591,263), sequences containing ‘‘N’’ (39,404) and low quality regions (1,601,445) resulted in 50,831,701 clean reads (95.6% of raw reads). These clean reads were further assembled into 56,927 unigenes with the N50 value of 3479 and mean size of 1695 bp (ranging from 201 bp to 37,991 bp, Fig. S1). The assembled S. litura transcriptome had call accuracy Q20 values of 97.45% and the GC (%) content was 47.45. Using BLAST against the NCBInr database we annotated 16,161 genes (44.67% of all unigenes). Addition of data from supplementary searches against alternative databases resulted in a total of 19,055 contigs, 15,900 of them longer than 200 bp, which have been deposited in NCBI as BioProject ID: PRJNA232395. 3.2. Functional classification and pathway analysis We used gene ontology assignment programs to assign 15,130 annotated contigs to three gene ontology classes (biological process, cellular component and molecular function), with 2163 functional terms assigned (Fig. 1). Cellular process, metabolic process, cell, cell part, binding and catalytic activity were the most common categories in each GO domain. Among categories with fewest members were cell proliferation, cell junction and receptor regulator activity. Using clusters of eukaryotic orthologous groups of proteins (KOG) we further classified the annotated contigs into 26 groups, representing major phylogenetic lineages (Fig. S2A). The results showed that the largest group (2145 contigs, >20% of total) contained only a general functional prediction, and the following group was signal transduction with 1331 contigs (approximately 15% of total). Mapping the annotated sequences to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database resulted in 7685 contigs matching in the database and assigned to 238 KEGG pathways (Fig. S2B). The most common pathways observed were signal transduction (about 9% of genes) followed by transport and catabolism, and carbohydrate metabolism (about 7% of genes). 3.3. Identification of Cry toxin receptor homologs in the S. litura transcriptome Most currently available functional data identify cadherins, APNs, ALPs, and ABCC2s as Cry toxin receptors in Lepidoptera (Adang et al., 2014). Predicted cadherin, APN, ALP, and ABCC2 contigs in the S. litura transcriptome were compared through phylogenetic analyses to proposed lepidopteran Cry1 toxin receptors (See Supplementary Figure Legends). Out of the 33 detected cadherin contigs in the transcriptome only comp20295 displayed relevant identity (76.6%) and clustered (Fig. S3) with a cadherin from S. exigua (AFH96949.1). That S. exigua cadherin is the only cadherin reported as Cry1Ca receptor (Ren et al., 2013). However, reduced expression (about 80%) of

Fig. 1. Gene function classification (Gene Ontology or GO terms) for unigenes annotated in the S. litura transcriptome. Shown are the identified functions and the corresponding percent and total number of unigenes for each GO category.

the S. exigua cadherin by gene silencing only resulted in a 50% reduction in susceptibility to Cry1Ca, indicating that additional Cry1Ca receptors exist in that insect (Ren et al., 2013). In BLASTx searches of the NCBInr database (data not shown), comp20295 was identified (99% identity) as a S. litura cadherin-like receptor (AFJ04291.1). The only other contigs displaying significant identity to cadherin proteins in the database were comp11997 and comp12002_c2, which displayed 54% (E-value 2e-29) and 92% (E-value 1e-69) sequence identity to a mutated cadherin from Cry1Ac-resistant Helicoverpa armigera (ACY69027.1) and a protocadherin from Bombyx mori (XP_004924479.1), respectively. Based on these observations, we selected comp20295, as a putative Cry toxin receptor, and contigs comp11997 and comp12002_c2 as having sequence identity but not clustering with Cry receptor cadherins, for expression analyses. Sequence annotation identified 24 predicted APN contigs in the S. litura transcriptome, of which only 7 clustered with Cry toxin

4

L. Gong et al. / Journal of Invertebrate Pathology 129 (2015) 1–6

receptor APNs (Pigott and Ellar, 2007) in phylograms (Fig. S4). Contig comp19222 clustered with APNs in Class 4, and was identified as representing (>90% identity) the SlAPN protein (AAK69605.1) previously described as a receptor for Cry1Ca toxin in S. litura (Rajagopal et al., 2002). Interestingly, the Cry1Ca-binding epitope in SlAPN (Kaur et al., 2014) is not found in any other S. litura APN or APN contig in our transcriptome (Fig. S5). The next contig displaying highest identity to Bt toxin receptors was contig comp18626_c1, which displays 75% identity to a Class 1 APN not expressed in a Cry1C-resistant strain of S. exigua (Herrero et al., 2005). Other contigs displaying lower identity to Cry toxin receptor APNs included contig comp19608 grouped with Bt receptor APNs in Class 2, contig comp13857 clustered with Ostrinia nubilalis ACV04931.1 in APN Class 8, contig comp19217 clustered with Bt toxin receptors in APN Class 3, and comp19516 clustering with Plutella xylostella CAA10950.1 in APN Class 5. Contig comp22241 also clustered in APN Class 5, but displayed lower identity (