Down-regulation of aminopeptidase N and ABC transporter subfamily G transcripts in Cry1Ab and Cry1Ac resistant Asian corn borer, Ostrinia furnacalis (Lepidoptera: Crambidae).

The Asian corn borer (ACB), Ostrinia furnacalis (Lepidoptera: Crambidae), is a highly destructive pest of cultivated maize throughout East Asia. Bacillus thuringiensis (Bt) crystalline protein (Cry) toxins cause mortality by a mechanism involving pore formation or signal transduction following toxin binding to receptors along the midgut lumen of susceptible insects, but this mechanism and mutations therein that lead to resistance are not fully understood. In the current study, quantitative comparisons were made among midgut expressed transcripts from O. furnacalis susceptible (ACB-BtS) and laboratory selected strains resistant to Cry1Ab (ACB-AbR) and Cry1Ac toxins (ACB-AcR) when feeding on non-Bt diet. From a combined de novo transcriptome assembly of 83,370 transcripts, ORFs of ≥ 100 amino acids were predicted and annotated for 28,940 unique isoforms derived from 12,288 transcripts. Transcriptome-wide expression estimated from RNA-seq read depths predicted significant down-regulation of transcripts for previously known Bt resistance genes, aminopeptidase N1 (apn1) and apn3, as well as a putative ATP binding cassette transporter group G (abcg) gene in both ACB-AbR and -AcR (log2[fold-change] ≥ 1.36; P < 0.0001). The transcripts that were most highly differentially regulated in both ACB-AbR and -AcR compared to ACB-BtS (log2[fold-change] ≥ 2.0; P < 0.0001) included up- and down-regulation of serine proteases, storage proteins and cytochrome P450 monooxygenases, as well as up-regulation of genes with predicted transport function. This study predicted the significant down-regulation of transcripts for previously known Bt resistance genes, aminopeptidase N1 (apn1) and apn3, as well as abccg gene in both ACB-AbR and -AcR. These data are important for the understanding of systemic differences between Bt resistant and susceptible genotypes.

Introduction

The productivity of cultivated crop plants has increased in part due to technologies that effectively suppress the feeding damage caused by herbivorous insect pests. This has involved the widespread use of chemical pesticides and genetically modified crops that express transgenic insecticidal protein toxins derived from the soil bacterium, Bacillus thuringiensis (Bt; 1). The widespread applications of chemical pesticides are implicated as likely causes for the increased incidence of functional resistance in insect populations 2, and has analogously occurred through selection pressures imposed by Bt crops 3. The control of several insect species has been compromised due to the evolution of resistance to Bt toxins, and include the maize pests the African stem borer, Busseola fusca 4, the Fall armyworm, Spodoptera frugiperda 5-8, and the corn rootworm, Diabrotica virgifera virgifera 9, 10, as well as the cotton pests Helicoverpa armigera, H. zea, and Pectinophora gossypiella 11-15. Other reports have documented a significant reduction in H. zea susceptibility towards the Cry1Ab toxin levels expressed by transgenic maize in the field 16. Analogously, over the past decade the survivorship of field-collected western bean cutworm (WBC), Striacosta albicosta, on Cry1F toxin in bioassay laboratory has increased 17, where larvae show little growth inhibition even at high toxin exposures 18.

Investigations into the biochemical, molecular and genetic mechanism by which Bt resistance evolve has identified a great variety of different potential mechanisms 19-21, which are thought to inhibit the normal Bt mode of action 22, 23. Membrane bound glycosylphosphatidylinositol (GPI) anchored aminopeptidase N (APN) was among the first proteins identified in the midgut of Lepidoptera capable of binding to Cry1A toxins in vitro 24, 25. The involvement of apn genes in the Cry1A and Cry9C toxin modes of action was suggested following knockdown by RNA interference (RNAi) that resulted in reduced susceptibilities among lepidopteran larvae 26-28. Implication of endogenous apn1 mutations in the evolution of Bt toxin resistance was shown to occur via the suppression apn1 transcription in Cry1Ca resistant S. exigua 29, Cry1Ac resistant T. ni 30, and Cry1Ab tolerant O. nubilalis 31. Similarly, levels of apn3 transcripts were reduced in Cry1Ab resistant O. nubilalis 31 and Cry1Ac resistant O. furnacalis 32. Initial indication that an ATP binding cassette (ABC) transporter may be involved in Bt toxin resistance traits came from QTL mapping results that linked insertion of a premature stop codon in abcc2 to Cry1Ac resistance in H. virescens 33. Different sets of deletions and point mutations in abcc2 orthologs were subsequently associated with resistance of P. xylostella and T. ni to Cry1Ac 34, B. mori resistance to Cry1Ab 35, and O. nubilalis survival on transgenic Cry1Fa maize 36. Furthermore, abcc2 and abcc3 transcripts were linked to Cry1Ac and Cry1Ca toxins resistance in S. exigua 37, as well as Cry1Ac resistance in the P. xylostella strain DBM1Ac-R 38. The repertoire of ABC transporters likely involved in Bt toxin resistance has expanded to implicate the pigment transporting ABCG1, Pxwhite, that shows down-regulated transcription in Cry1Ac resistant P. xylostella 39. Furthermore, a truncation of the abca2 gene of H. armigera was linked to Cry2Ab resistance 40, and spliceosome variants in H. armigera abcc2 are associated with Cry1Ac resistance 41.

A cell adhesion protein, cadherin (cad), was identified as an in vitro Cry1A toxin receptor in the midgut of M. sexta 42, and an ortholog from O. nubilalis analogously is capable of binding Cry1Ab, Cry1Ac and Cry1F toxins in vitro 43, 44 as well as Cry1A toxins in vivo 45. A mutation caused by transposon insertion into the H. virescens cad provided the initial evidence for involvement in resistance mechanisms 46, and additional mutations were associated with Bt resistance traits among field-derived Pectinophora possypiella 12, H. armigera 47, and O. furnacalis 48. Despite these prior findings, the cadherin locus is not linked to Cry1Ab or Cry1F resistance traits in O. nubilalis 31, 36, or Cry1Ac resistance in T. ni 49. Analogously, isoforms of membrane-bound alkaline phosphatase (ALP) in the midgut of M. sexta and H. virescens bind the Cry1Ac toxin 50, 51, and the down-regulation of alp was subsequently associated with resistance to Cry1Ac in H. armigera and Cry1Fa in S. frugiperda 52 and Cry1Ac in P. xylostella 38. Interestingly, Guo et al. demonstrated that alp down-regulation may be controlled by pathways involving a mitogen-activated protein kinase, MAPK4 38. Experiments on O. furnacalis show that although ALP is a putative receptor for Cry1Ac in the midgut of susceptible individuals, but neither toxin binding nor expression of alp transcripts were reduced among Cry1Ac resistant larvae 53.

The Asian corn borer, O. furnacalis (Guenée), is a highly destructive pest insect of cultivated maize and is endemic to parts of East Asia and Australia, including China 54, 55. Although chemical insecticides and natural enemies are used for control 56, 57, these methods tend to show limited effectiveness due to difficulties in the timing of applications or environmental variation. Similar to other global regions, Bt transgenic crops hold great promise for the improvement of agricultural production in China, but concerns remain regarding potential evolution of resistance voiced nearly two decades ago 58 and now being observed amongst several insect pest populations 3, 59. Laboratory strains of O. furnacalis have been selected for high levels of resistance to Cry1Ab 32, 60 and Cry1Ac 61, but a complete understanding of the genetic mechanisms of these resistance traits remain incomplete. Recent advances in next generation DNA sequencing technologies allow the analysis of variation in transcript abundance (gene expression) through the application of RNA sequencing (RNA-seq) 62, and have allowed the comparison gene regulation between Bt resistant and susceptible phenotypes 63, 64. The following paper describes construction of a reference midgut transcriptome to which mapped RNA-seq read data were used to estimate gene expression differences among Bt resistant (Cry1Ab and Cry1Ac) and susceptible stains of O. furnacalis. These methods, in conjunction with secondary validation, demonstrate the capability of RNA-seq methods to identify candidate genes putatively involved in the evolution of Bt resistance traits. These data are important for the understanding of systemic differences between Bt resistant and susceptible genotypes, and how selection may result in fixed changes of expression level among numerous genes.

Materials and Methods

Insect resistance levels

Three O. furnacalis laboratory colonies were used in the current study; a Bt toxin susceptible colony (ACB-BtS), and two laboratory colonies selected for increased resistance to Cry1Ab (ACB-AbR) and a Cry1Ac toxins (ACB-AcR). Selection schemes were described previously for ACB-AbR 32, 60, and ACB-AcR 61. In brief, O. furnacalis larvae were collected from maize fields throughout central China, and reared on artificial diet with a photoperiod of L:D=16:8 h and relative humidity at 70~80%. Larvae in ACB-AcR were initially exposed to 2.5 ng Cry1Ac toxin/g of artificial diet, and exposures steadily increased to levels that caused 40 (LC40) to 70% mortality (LC70) in subsequent generations. ACB-AcR was reared in subsequent generations a final concentration of 400 ng Cry1Ac/g artificial diet as described previously 65. Selection was analogously replicated in order to generate the ACB-AbR colony 32, 60.

Resistance to Cry1Ab and Cry1Ac was respectively estimated among ACB-AbR and ACB-AcR larvae as the proportion of moribund larvae observed following a 7 d exposure to diet incorporated dose-response assays, and reported as a resistance ratio (RR) that is compared to the mortality among ACB-BtS larvae. Specifically, neonates from each strain were exposed to 9 concentrations of Cry1Ac (0, 1, 50, 100, 200, 400, 600, 800 and 1000 μg/g) and 9 concentrations of Cry1Ab (0, 0.1, 0.5, 1, 5, 10, 50, 100 and 200 μg/g) for 7 days to obtain 5% to 100% mortality. Bioassays at each toxin level were performed in triplicate (n = 48 neonates per replication), and tap water was used in the control treatment. The lethal concentration resulting in mortality among 50% of individuals (LC50) was calculated for each strain and each toxin using a Probit statistical model in Polo Plus (LeOra Software, USA). For each Bt toxin, a RR was calculated by dividing LC50 of the ACB-AbR or ACB-AcR strain by the LC50 estimated from the susceptible strain, ACB-BtS.

Midgut dissection and RNA extraction

Larvae from each of the ACB-AbR, ACB-AcR, and ACB-BtS colonies were fed on non-Bt toxin diet under standard rearing conditions, and 10 live 5th instars were randomly selected from each and placed on ice. The midgut tissues were dissected from each and the lumen was rapid cleaned with a solution of 0.7% NaCl (w/v) to remove debris, followed by submersion in liquid nitrogen. Dissected midguts were pooled by colony (n = 10 per pool), and total RNA was extracted separately from each the ACB-AbR, ACB-AcR and ACB-BtS pool using the Trizol reagent following manufacturer instructions (Invitrogen, CA). Resulting RNA samples were treated with DNase I (Promega, Madison, WI, USA) following manufacturer's instructions. Quantity and quality of RNA was assessed by denaturing gel electrophoresis and spectrophotometry on a Nanodrop 2000 (Thermo Scientific, Wilmington, DE, USA).

Midgut transcriptome sequencing, assembly and functional annotation

Three single-end cDNA libraries were constructed, one for each of the three midgut total RNA samples derived from ACB-BtS, ACB-AbR, and ACB-AcR. Purified total RNA from each sample was separately converted to complementary DNA (cDNA) using the ScriptSeq™ mRNA-seq Library Preparation Kit (Illumina, San Diego, CA) following the manufacturer's protocol. In brief, mRNA was isolated from 20 µg of purified total RNA using oligo (dT) magnetic beads, and the resulting mRNA fragmented by incubation with divalent cations at 94 °C for 5 min. First strand cDNA was generated by reverse transcription primed by random hexamers, followed by synthesis of the second-strand cDNA using RNase H and DNA polymerase I. After end-repair and ligation of Illumina sequencing and flow cell adaptors, the products were amplified by PCR and purified with QIAGEN MiniElute PCR Purification Kit (Qiagen, Venlo, Netherlands), and single-end 100 bp sequence read data was generated on an Illumina HiSeqTM 2000 platform at Novogen Co., Ltd. Beijing, China. Raw Illumina read data were submitted to the National Center for Biotechnology Information (NCBI) short read archive (SRA) database (Table 1; ACB-BtS: SRX751199, ACB-AbR: SRX751197, ACB-AcR: SRX750576) under BioProject PRJNA266299.

Table 1

Lethal concentration (LC) and resistance ratio (RR) estimates for Cry1Ab and CryAc toxins among Ostrinia furnacalis larvae from ACB-AbR and ACB-AcR strains compared to those in the susceptible ACB-BtS strain.

Bt toxin	ACB strain	LC50 (95% CI)* μg/g	RR
Cry1Ab	ACB-BtS	0.36 (0.13-0.51)	/
	ACB-AbR	69.77 (49.10-98.66)	193.80
Cry1Ac	ACB-BtS	0.27 (0.19-0.34)	/
	ACB-AcR	967.76 (842.75-1092.81)	3584.30

* 95% confidence interval

Read statistics and quality assessment was conducted using FastQC 66. Contaminating Illumina adapters and nucleotides within a 20 bp window with a mean Phred 33 quality score less than 20 (q < 20) were trimmed from the raw data using Trimmomatic 0.32 67. The filtered high quality fastq read data from all libraries were separately used as input for the Trinity in silico read normalization script 68 using default parameters except a kmer size (--KMER_SIZE) was set to 25 and maximum coverage per kmer (--max_cov) limited to 30. Trinity r20140414 69 was used to assemble contigs (Inchworm), cluster contigs (Chrysalis) and reconstruct transcript isoforms (Butterfly) according to Haas et al. 70. The resulting contigs were submitted to the National Center for Biotechnology Information (NCBI) transcriptome shotgun assembly (TSA) database under RNA-seq Assembly project (r20131110). Open reading frames (ORFs) of ≥100 amino acids were predicted from Trinity fasta output using TransDecoder 1.0, and the resulting .cds file used as the reference for subsequent RNA-seq analyses. Putative functional annotations were obtained for all O. furnacalis derived proteins sequences obtained from the TransDecoder 1.0 .pep files using Blast2GO 71 as described by Garczynski et al. 72, except that BLASTp 73 was used to query the NCBI non-redundant (nr) protein database.

Estimates of differential-expression among midgut transcripts

Non-normalized trimmed single-end fastq read data derived from each Illumina RNA-seq library were each mapped separately to the ORF-containing transcripts identified by TransDecoder (.cds file) using Bowtie 2.2.1 74 with default parameters and alignments output in sequence alignment map (SAM) file format. Each SAM file was converted to a binary sequence alignment map (BAM) format and counts of the number of reads aligned to each ORF-containing contig within indexed BAM files were output a tab-delimited text using SamTools 75. The level of expression was estimated for each transcript among replicated ACB-BtS and ACB-AbR and ACB-AcR alignments using read counts as a proxy for transcript abundance 76, and implemented using the DESeq package 77 with the false discovery rate (FDR) set at 0.001 and quantile for removing insignificant estimates at 0.4. MA-plots comparing log2[fold-change] to mean of normalized read counts were generated, as well as a plot showing the distribution of p-values. Significance threshold were set at two levels 1) a moderate threshold of log2[fold-change] ≥ |1.0| and p-values ≤ 0.001, and 2) a high stringency threshold of log2[fold-change] ≥ |2.0| and p-values ≤ 0.001.

Nucleotide and derived protein sequences corresponding to each Trinity-defined transcript isoform (component; comp from Butterfly output) within putatively differentially expressed transcripts [log2[fold-change] ≥ |1.0| and p-values ≤ 0.001] were obtained from the reference transcriptome using the PERL script select Seqs.pl. Transcripts putatively annotated within our Blast2GO analyses as candidate Bt resistance genes were used as queries of all translated amino acid sequences using the tBLASTn algorithm (identity cutoffs ≥ 98% over ≥25 amino acids), and consensus transcripts were reconstructed used CAP3 78. Derived protein sequence of each CAP3 assembled contig was obtained using the Virtual Ribosome 79 and reading frame adjusted based on BLASTx search of the NCBI nr protein database. In all cases, the nearest ortholog from the model lepidopteran B. mori were downloaded in FASTA format. Additionally, homologous sequences from the O. nubilalis midgut transcriptome were identified by BLASTp query of a local database composed of translated fasta sequences obtained from Vellichirammal et al., 64. Multiple sequence alignments performed among B. mori, O. nubilalis and O. furnacalis sequences with ClustalW v2.0 using default parameters 80. Conserved protein primary and secondary structures were identified by search of the Conserved Domains Database (CDD) 81, and use of the Smart Modular Architecture Research Tool (SMART) 82, and Transmembrane protein topology Hidden Markov Model (TMHMM) prediction algorithm 83.

Three closely related aminopeptidase N paralogs are known within the genome of O. nubilalis (OnAPN3a, OnAPN3b, and OnAPN3c) that show high degrees of amino acid sequence similarity 31, 84. Differentially expressed O. furnacalis transcripts annotated as encoding APN3-like proteins were aligned with Onapn3a, Onapn3b, and Onapn3c orthologs using the ClustalW algorithm (default parameters), and imported into the program Molecular Evolutionary and Genetic Analysis (MEGA) 6.0 85. Gluzincin aminopeptidase, GAMEN, and zinc binding motifs (HEX2HX18E) were identified from CDD and SMART output as described above. Putative GPI-anchors signals were predicted using GPI-SOM 86. Phylogenetic inferences were made using the genetic distances with the Neighbor-joining (NJ) methods as well as by Maximum likelihood, and node support provided by 1,000 pseudoreplicates.

Quantitative RT-PCR (qRT-PCR) analysis was conducted to validate the abundance of O. furnacalis midgut transcripts that were annotated as candidate Bt resistance genes or that showed significant differential expression in RNA-seq experiments. The candidate Bt resistance genes alkaline phosphatase and cadherin which were not predicted to show significant up- or down-regulation from RNA-seq analyses were also tested by qRT-PCR. The O. furnacalis β-actin gene (accession number: GU301782) was chosen as a reference for normalizing qRT-PCR data. Oligonucleotide primers for qRT-PCR were designed with Primer Express 3.0 (ABI, USA) (). Total RNA was extracted from midgut of ACB-AcR, ACB-AbR and ACB-BtS strains as described previously and 2.0 µg used in 1st strand cDNA synthesis. Each strain including 30 larvae, and each 10 larvae were treated as biological duplicates. Reactions were setup in triplicate using 20 µL reaction volumes containing 10 µL of SYBR Green Universal Master Mix II (Takara, Japan), 1.0 µL of PCR Forward Primer (10 μM), 1.0 µL of PCR Reverse Primer (10 μM), 0.5 µL of ROX reference dye II (50X), 1.0 µL of 1st strand cDNA as template, and 6.5 µL sterilized ultra-pure grade water. Amplification was performed using the ABI 7500 Fast Real-Time PCR System, at 95 °C (30 s), followed by 40 cycles at 95 °C (10 s) and 60 °C (20 s) (2 steps per cycle), and then melt-curve analysis to evaluate the specificity of the amplifications. Normalization of the experimental cycle threshold (CT) data relative to the β-actin reference gene was done as previously described 87, 88, and relative expression of normalized data was made to values from the ACB-BtS strain by the comparative 2-∆∆CT method 89. The significance of any differences in estimated CT values between experimental unknowns was assessed by Nested ANOVA, and thresholds set at α = 0.05.

Structural variant predictions

The individual .BAM alignment generated for ACB-BtS, -AbR and AcR above were merged into as single file using the SamTools mpileup command. The resulting merged BAM file was interrogated for variant nucleotides with respect to the reference transcriptome using SamTools and BCFtools as described previously 90, and vcfutils.pl VarFilter options filtered the variant call format (.vcf) output file for read depths ≥ 20 and for variants ≥ 10 bp from any alignment gap. The vcf v4.1 file was further filtered for variant sites fixed differently between ACB-BtS and -AbR or -AcR, and mutations causing non-synonymous (amino acid) changes were predicted from annotated CDS data.

Results

After 142 generations of selection for survival on Cry1Ab toxin diet, the ACB-AbR laboratory colony has an estimated RR of 193.8 compared to the susceptible ACB-BtS strain (Table ), which is an increase compared to the RR of 106.8 previously reported after 51 generations of selection 60. Analogously, the ACB-AcR strain shows an estimated 3,584.3-fold increase in Cry1Ac resistance compared to ACB-BtS after 144 generations of selection on Cry1Ac toxin, diet (Table ).

A reference O. furnacalis midgut transcriptome was assembled from paired end mRNA-seq data derived from ACB-AbR, ACB-AcR and ACB-BtS strains reared on control artificial diet (without Bt toxin), and therefore was useful to investigate fixed (constitutive) differences in gene expression between strains. The combined set of libraries for each strain contained ≥ 41.2 million reads and ≥ 4,117 Mbp of sequence data (Table A). Following removal of adaptor, ambiguous and low-quality sequence data with a Phred 33 quality score < 20 (q < 20), ≥ 35.2 million reads and ≥ 3,400 Mbp of sequence data remained from each strain. Following read normalization, assembly of contigs, contig clustering, and isoform reconstruction using the Trinity pipeline, a combined midgut reference assembly of 83,370 O. furnacalis midgut-expressed transcripts was obtained. The assembled sequences were deposited in the NCBI Transcriptome Shotgun Assembly (TSA) database and linked to the BioProject PRJNA266299. These sequences do not represent unique genes, but instead contain a total count of transcript isoforms (alternative splice variants or closely-related paralogs) defined within 41,088 unique components (comp; defined by the program Chrysalis) within a set of transcripts by Butterfly (Table B; see Haas et al., 70 for additional details). From this reference assembly, a total of 12,288 unique coding sequences were identified from 28,940 components (isoforms groups; Table B).

A transcriptome was previously assembled from whole O. furnacalis larvae and used to predict putative DEGs between ACB-AbR and -BtS 32, but the combined reference assembly in this study represents the first midgut-specific genomic resource for this species, and provided the means by which to assess expression-level differences within midgut tissues. Nearly 60% of the 28,940 coding sequences predicted within the O. furnacalis midgut transcriptome did not receive any functional gene annotations. Although problematic for efforts to obtain putative homology-derived functions for our dataset, these data may represent potential genomic novelty what could be applied in future efforts to control crop damage by this insect. Regardless, annotation, GO enrichment and pathway analyses were possible for a number of transcripts (Table ).

The goals of this study were to determine the transcripts and functional gene categories that are putatively involved in Cry1Ab and Cry1Ac resistance mechanisms. This was addressed by estimating transcript levels in the midgut of Cry1Ab and Cry1Ac resistant O. furnacalis based on counts of mapped RNA-seq read data, followed by application of statistical models to estimate quantities that differ significantly when compared to the midgut of susceptible larvae. Statistically significant fold-changes in transcript levels remain difficult to interpret within a biological context, but we follow the precedent of a large number of prior RNA-seq studies that assume that there is a biological consequence at the cellular or functional level. In addition, our data estimated differential expression among two colonies independently selected for Cry1Ab and Cry1Ac toxin resistance, wherein we interpreted constitutive up- or down-regulation in both colonies as strong evidence for involvement of those genes in the evolution of Cry1A toxin resistance.

To assess the significance of any differences in estimated midgut transcript levels between resistant and susceptible larvae, the current study applied two criteria within our DESeq analysis; Threshold 1 of log2[fold-change] ≥ |1.0| and p-values ≤ 0.001 that identified ≤ 1,701 significantly up- or down-regulated transcripts between ACB-AbR or -AcR and ACB-BtS phenotypes (Fig. ), and ≤ 392 transcripts at a more stringent criteria of log2[fold-change] ≥ |2.0| and p-values ≤ 0.001 (Fig. ). This comparison showed that a larger proportion of genes were significantly down-regulated compared to up-regulated, and down-regulated transcripts included those encoding aminopeptidase N paralogs and a member of the ABC transporter family (ABCG).

Additional genes highly down-regulated among Cry1Ab and Cry1Ac resistant larvae included those in putative functional categories involved in lipid metabolism and transport, binding, and stress response (Fig. . Some of the annotations that were made among the down-regulated transcripts for lipid storage and transport proteins (arylphorin, methionoe-rich protein), lipid metabolizing enzymes (e.g. lipases), and stress-response associated genes glutathione S-transferase (GST), cytochrome P450 monooxygenases, and heat shock proteins (HSP) 68- and 70-like (). Several uridine diphosphates (UDP)-glucuronosyltransferase-like transcripts were also significantly repressed in the midgut of both Cry1Ab and Cry1Ac resistant O. furnacalis.

Also among the significantly down-regulated transcripts in the high stringency threshold group (log2[fold-change] ≤ -2.0 and p-values ≤ 0.001), several were annotated as previously known Bt resistance genes (see Introduction) including members of the aminopeptidase N gene family, apn3 paralogs and apn8, an ABC transporter in subgroup G, abcg, as well 10 serine protease genes. Sequence alignment of derived APN3 paralogs from species of Ostrinia indicated ≥ 96% intraspecific amino acid identity (Fig. ) and subsequent phylogenetic analyses resolved the association between paralogs with APN3a and 3c predicted to be the most recently diverged duplicates (Fig. ). Estimates of significantly reduced expression in RNA-seq data (Table ) were corroborated by qRT-PCR results for Ofapn1 and Ofapn8, as well as for Ofapn3b, in the midgut of both ACB-AbR and -AcR (Fig. ). Candidate Bt resistance genes were also annotated within the lower stringency threshold group of O. furnacalis transcripts (≥ 1.0 log2[fold-change] ≤ -1.0), and included those encoding ABC transporters from three different subgroups (A, C, and G), ABCA2, ABCC2, and ABCG4. The number of significantly up-regulated transcripts were smaller in number for both Cry1Ab and Cry1Ac resistant O. furnacalis compared to down-regulated genes (Fig. ), and interestingly included two up-regulated ABC transporter subgroup G members (comp15850 and comp9932 (scarlet); Table ).

The putative ABCG1 encoded by O. furnacalis comp3109_c0_seq1:140-646(+) is the ABC transporter that was the highly down-regulated in both Cry1Ab and Cry1Ac resistant strains based on RNA-seq and qRT-PCR. The derived 168 amino acid sequence from the 646 bp O. furnacalis abcg1-like transcript (comp3109_c0_seq1:140-646(+)) shows ≥ 71% amino acid sequence identity with the putative ACBG1-like transporters from B. mori (XP_004929563.1; 667 aa; aligned positions 510 to 667) and P. xylostella (XP_011551601.1; 669 aa; aligned positions 502 to 669). These orthologs correspond to gene models KAIKOGA014788 on B. mori chromosome 12 and Px012058 on P. xylostella scaffold 47 (Fig. ). Interestingly, the 386 aa sequence encoded by the B. mori gene model BGIBGA010555 spans the C-terminal half of XP_004929563.1 (positions 281 to 667) as well as Px012058 (positions 283 to 669). Similarly, the 2,289 bp the comp48645_c0_seq2 from an O. nubilalis midgut transcriptome 64 encodes a 345 aa long variant located at the C-terminus of the ortholog (). Annotation of the O. nubilalis OnBAC 35G22 full insert sequence (GenBank HTGS KX793137) using these data defined a putative 1,773 bp CDS fragment encoding a 590 aa partial OnABCG1 peptide, with the missing the N-terminus residing on another BAC (Fig. ). Alignments to this genomic CDS showed that homology with truncated transcripts began at exon splice junctions, thus providing evidence that truncated copies are indeed splice variants (Fig. ; ). Analogously, mapping of O. furnacalis RNA-seq reads to the 1,773 bp partial Ofabcg1 transcriptional unit on OnBAC 35G22 provided evidence that the Ofabcg1 comp3109_c0_seq1:140-646(+) transcript is the major isoform in the midgut. The above RNA-seq evidence and qRT-PCR data confirm that this Ofabcg1 isoform is down-regulated in Cry1Ab and Cry1Ac resistant larvae (Fig. ). SMART and TMHMM secondary structural predictions indicated that the 168 aa OfABCG1 contains three C-terminal transmembrane domains (TMD4, TMD5, and TMD6) whereas the OnABCG1 homolog retains all six, but both have complete NBD deletions.

OfABCG4-encoding transcripts are also predicted to be down-regulated in the midgut of Cry1Ab and Cry1Ac resistant larvae (Table ; comp17719_co_seq4:1-1962(+)), which has the full OfABCG4 CDS (). In contrast, qRT-PCR results suggest that Ofabcg4 transcript levels are not significantly different in ACB-AbR or -AcR compared to ACB-BtS (Fig. ). No Ofabcg4 transcription isoforms (variants) were identified. Four transcripts from a single OfABCA2-like CDS were differentially regulated among resistant strains, and most significantly reduced among Cry1Ab resistant larvae (Table ). No O. furnacalis transcripts were found that encoded a full-length OfABCA2 peptide, but the entire sequence was reconstructed from fragments using CAP3 and shown to share > 97% identity with the translated sequence from O. nubilalis Oncomp46431_c0_seq2 (). Analogously, CAP3 reconstructed the complete CDS for the Ofabcc3 transcript that was significantly down-regulated in the ACB-AbR colony (Table ; ). ABC transporter functional domains were identified within all aligned peptides. Transcripts for abcc3 was the only ABC transporter family C member predicted to show significant differential regulation within RNA-seq analyses and only occurred via repression within the ACB-AbR colony (Table ), which was supported by qRT-PCR data (Fig. ). In contrast to our RNA-seq results, qRT-PCR indicated that Ofabcc1 and Ofabcc2 are up-regulated in both resistant strains, with the exception of Ofabcc2 in ACB-AbR.

Although the numbers of up-regulated transcripts are comparatively less in number, several genes in the same functional categories were shared with those that were down-regulated (Fig. . Specifically, these include a number of serine proteinases (Table ) as well as and lipid storage and metabolizing proteins, and GST (). Interestingly, transcripts encoding candidate Bt binding receptors a membrane-bound alkaline phosphatase (comp13277_c0_seq1:1-414(+)) was not significantly down-regulated in RNA-seq and corroborated by qRT-PCR (Fig. ), in contrast cadherin was significantly down-regulated in qRT-PCR but not within RNA-seq results. Cry1Ab and Cry1Ac resistant larvae share the putative significant up-regulation of the ABCG-like scarlet gene, but was barely over the threshold (log2[fold-change] ≤ +1.08).

A total of 12,487 variant sites with genotype quality score of 999 were predicted within transcripts across all RNA-seq alignments to the reference transcriptome (data not shown). Filtering of this .vcf output for homozygous variants resulted in 51 sites (50 SNPs and 1 indel), of which 30 were fixed between ACB-BtS and both -AbR and -AcR. Among the other sites, 20 were fixed between ACB-AbR compared to both -BtS and -AcR, and one SNP was fixed within ACB-AcR (). Annotations defined 48 putative protein coding regions, in which SNPs predicted 5 non-synonymous changes. None of the fixed variant sites were predicted within transcripts from DEGs for any known candidate Bt resistance gene (Table ).

Discussion

Laboratory selected colonies remain effective models for the study of field-evolved resistance, in that mutations identified in selected colonies have also been detected among resistant phenotypes in the field 47, 91. Furthermore, lab selected strains have been instrumental in the identification of genes involved in the Bt toxin mode of action as well as resistance mechanisms. Potential field relevance of this O. furnacalis trait conferring a 193.8-fold increase in resistance to Cry1Ab was demonstrated by neonate survival and development to pupation when feeding on silk tissues from the Cry1Ab-expressing transgenic event Mon810 60, which compares to 100% mortality of an unselected field strain 92. The ACB-AcR strain has a higher estimated Cry1Ac resistance ratio (3,584.3-fold increase compared to ACB-BtS) after 144 generations of selection (Table ), but commercial maize hybrids expressing Cry1Ac toxin are not available in order to perform on-plant assays for survivorship. High fold-increases in resistance have analogously been estimated for these ACB-AbR and -AcR strains (Table 32, 61, which agrees with prior evidence that chronic Bt toxin exposures can lead to incremental increases in the tolerance towards Bt toxins 93, 94. Granted chronic low-level exposures such as those used for O. furnacalis and O. nubilalis may favor the selection of resistance traits that are polygenic and potentially influenced by epistasis (gene-gene interactions) 95); a premise that was partially supported by the detection of more than one quantitative trait locus (QTL) influencing Cry1Ab tolerance in O. nubilalis 31 as well as maybe even for Cry1Ac resistance in H. virescens 46. Moreover, these low-dose exposures may be pertinent for modeling Bt resistance traits that develop among lepidopteran species that feed on transgenic plant tissues with low Bt titres 96.

The Cry1Ab resistant ACB-AbR strain was reported to have a high degree of cross-resistance to Cry1Ac 60, 61, which is analogous to levels of cross-resistance to Cry1Ac among Cry1Ab selected O. nubilalis 97, 98 and Heliothis virescens 99. These observations may be due to Cry1Ab and Cry1Ac sharing common midgut receptors 100, 101, that could be supported by results of competitive toxin binding assays 102 as well as immunoblotting data from O. nubilalis that demonstrate binding to similar midgut protein receptors 44, 103. Surprisingly, ACB-AcR shows low cross-resistance towards Cry1Ab 104, but is analogous to the lack of cross-resistance to Cry1Ab among Cry1Ac selected P. xylostella 105. This might suggest that although Cry1Ab and Cry1Ac may share some of the same midgut receptors, the molecular mechanism(s) involved in the development of these two O. furnacalis resistance traits may not be identical.

Sequencing of cDNAs and subsequent assembly of these data into transcripts is an effective means to discover novel genes 106, and has been greatly facilitated among non-model organisms by recent advances in high throughput RNA sequencing (RNA-seq) 107. Due to complexities of assembling a heterogeneous mixture for differential splicing and allelic variants within RNA-seq data, the number of unique transcripts assembled in our O. furnacalis reference transcriptome (n = 83,370) is likely greater than the number of genes expressed in the midgut. Despite of this, analogous numbers of transcripts have been assembled and used successfully to estimate levels of differential expression among the non-model lepidopteran species Lymantria dispar 108, P. xylostella 63, and O. nubilalis 64. Additional difficulties pertaining to gene annotation are often encountered with the study of non-model organisms. Specifically, our results showed that nearly 60% of the O. furnacalis midgut transcripts failed to receive any functional annotation, but are somewhat typical among non-model organisms 93. The results may be influenced by the high proportion of transcripts belonging to the class of long non-coding RNAs (lncRNAs) that are transcribed as mRNA but do not encode open reading frames 109. The assignment of lncRNA origins from numerous and widespread portions of mammalian genomes suggest that most regions show transcriptional activity 110, and are highly prevalent among RNA-seq data from insects 111. In this era of high throughput sequencing computational annotations are relied upon to assign putative functional roles 112, and undoubtedly hinders the interpretation of the biological consequence of differentially transcription estimated from RNA-seq data 113. Despite these shortfalls, the current study did not address any of the above stated gaps in knowledge, but instead used available genomic resources and empirical data to infer the possible functions of putative differentially expressed genes.

RNA-seq experiments provide the opportunity to take a systems approach for evaluating global gene regulation, and can estimate significant genetic differences between phenotypes 114, 115. Our statistical analysis of RNA-seq data and subsequent annotation of significantly up- and down-regulated genes demonstrate the over-representation of certain functional categories including binding, metabolism and transport in Bt resistant midgut tissues (Fig. ). Closer interrogation of functional annotations show the presence of transcripts encoding peptides likely involved in lipid binding and transport, as well as lipid and protein metabolism. Since it remains difficult to investigate all of these differentially-regulated genes within the context of Bt resistance mechanisms, we used criteria that assumed genes and gene functional categories that were dysregulated between the two Cry1A toxin resistance traits to indicate the likely presence of a biochemical basis, or contrastingly, unique to once trait or the other. Moreover, the dearth of empirical evidence from non-model organisms in regards to defining gene function led both o the inability to derive putative functions for many transcripts as well as our confidence that all annotations are absolutely correct. Regardless, several transcripts estimated to be significant down-regulation in the midgut of both Cry1Ab and Cry1Ac resistant O. furnacalis, included members of large gene families previously identified as candidate Bt resistance genes; aminopeptidase N, ABC transporter subgroups, and serine proteases family members (Table ). Upon the addition of qRT-PCR data, our results also suggested that the O. furnacalis cadherin (Ofcad) is also repressed in ACB-AbR and -AcR strains (Fig. ) and agrees with previous qRT-PCR results from the ACB-AcR train 48. Results from our study and prior research on ACB-AcR 48 are difficult to compare due to use of whole larval RNA as template in the latter, suggesting any conclusions might be confounded by tissue-specific differences. The lack of any reduction in the expression of O. furnacalis alkaline phosphatase is analogous to prior results using the ACB-AcR strain 53.

Down-regulation of aminopeptidase N transcripts in the midgut of Cry1Ab and Cry1Ac larvae: The involvement of APN receptors in the Bt mode of action and mechanisms of resistance is well documented (see Introduction). There are eight orthologous apn genes, apn1 thru apn8, which resulted from ancestral tandem duplication in the lepidopteran lineage 116-118, and confirmed to occur in O. nubilalis by genetic linkage and physical mapping 31. Of these orthologs, apn1 and apn3 are identified as Cry1A toxin receptors in O. nubilalis 44 and O. furnacalis 43. Experimental evidence also shows that the O. nubilalis apn1 (Onapn1) gene product expressed in the surface of cultured cells bind Cry1Ab toxins in vivo, and OnAPN3a and OnAPN8 peptides analogously bind Cry1Fa 84. Results of the current study showed that Ofapn1, Ofapn8, as well as Ofapn3b are down-regulated in the midgut of larvae from both ACB-AbR and -AcR strains (Table ). Association of reduced Ofapn1 transcript levels with Cry1Ab and Cry1Ac resistance may not be surprising given linkage to Cry1A resistance traits in O. nubilalis 31 and T. ni 30, or even our implication of Onapn3 due to prior evidence from O. nubilalis 31 and whole larval RNA-seq data from ACB-AcR 32. Comparatively unique to all prior studies, our analysis allowed partitioning of RNA-seq data among the Ofapn3a, Ofapn3b, and Ofapn3c paralogs (Table ; Fig. ). Furthermore, our results from transcript-specific qRT-PCR suggest that Ofapn3b may be the paralog that is highly down-regulated among Cry1A toxin resistant Ostrinia. What remains interesting is that Onapn3 and Onapn8 bound Cry1Fa toxin but not Cry1Ab when expressed singularly in cell culture, whereas Onapn1 bound Cry1Ab and not Cry1Fa 84. This presents a potential discordance between binding assay results and genetic linkage or association studies. Specifically, considering that Crava et al., 84 identified OnAPN3a and OnAPN8 protein products as likely receptors for Cry1Fa, but the involvement of APN8 and APN3a in Cry1A toxin resistance mechanisms remains uncertain. Since expression of APN8 was not quantified in the midgut of Cry1A resistant T. ni 30 or O. nubilalis 31, the novelty of this association in O. furnacalis is unknown.

Our evidence shows that several APN gene family members are down-regulated in Cry1Ab and Cry1Ac resistant larvae. This might suggest that > 1 APN receptor might affect toxin binding wherein sequential binding could expose novel receptor epitopes 22, or alternatively that binding is not always indicative of a role in the toxin modes of action 119. Additionally, it was previously proposed that Onapn1 and Onapn3 may be part of a co-regulated gene network 31, wherein similar cis- and trans-acting factors may result in the coordinated transcription of gene family members. This premise may be partially supported by tissue specific expression of seven apns in B. mori 118, suggesting that tissue expression patterns may have been retained among gene duplicates. Thus is may be conceivable that multiple apn transcripts might be down-regulated as a result of a shared co-regulatory mechanism. Additional research is undoubtedly required in order to unravel the mechanisms that regulate apn expression in Lepidoptera, as well as how any perturbation in this transcriptional regulation can lead to the evolution of Bt resistance.

Down-regulation of transcripts for ABC transporters in the midgut of Cry1Ab and Cry1Ac resistant larvae: ATP-binding cassette (ABC) proteins are membrane-bound transporters implicated in the movement of solutes across lipid membranes and into the extracellular space 23, 120, 121. Gene coding regions are comprised of full transporter (FT) versions that have paired nucleotide binding domain (NBD) and transmembrane domain (TMD) or a half transporter (HT) that have one each of the NDB and TMD, and hetero- or homodimerization of the latter is required for functionality 122. Fifty one ABC transporter protein coding regions were identified in the B. mori genome and classified into within 8 subgroups by phylogenetic analysis (A-H) 123. Bt resistance mechanisms have implicated a number of ABC transporters from subgroups A 124, C 34-36, 46 and G 39 (see Introduction for additional details), and interestingly members from these three subgroups were down-regulated in Cry1Ab and Cry1Ac resistant O. furnacalis (Table ). Our qRT-PCR results confirm that the abundance of Ofabcg1 transcripts show the most significant proportional reduction in the midgut of Cry1Ab and Cry1Ac resistant larvae (Fig. ). The ABCG subfamily is composed of HTs that have an inverted arrangement of NBD and TMDs, and the functions of a vast majority of insect orthologs still remains unknown 125.

Regardless, the Drosophila melanogaster ABCG1 white that functions in pigment transport into the eye 126, was linked to P. xylostella Cry1Ac resistance 39. Also, regulation of 20-OH ecdysone response is regulated by the D. melanogaster ABCG1 E23 gene 127 and its ortholog from Tribolium castaneum 104. As suggested by Guo et al., 39 ABCG1 HTs may form a heterodimers with other subfamily G members, such that the down-regulation of the Ofabcg1 may have a large number of functional consequences. Although extrapolated through vast evolutionary time periods, several human ABCG members are involved in lipid transport 128, and the prediction that putative lipid metabolic and transport genes are differentially regulated in Cry1A resistant O. furnacalis might possibly suggest an analogous role of OfABCG1 and G4 in species of insects. If found to be true the down-regulation of ABCG transporters could likely affect an array of cellular functions, such that the expression of genes that are functionally related may also be affected. Specifically, putative reduced lipid influx or efflux from cells may result in a decreased requirement for lipid metabolism and intracellular transport, and further suggest that interpretation of RNA-seq data could be confounded by effects that are downstream of causal genetic factor - in this case reduced ABCG1 and G4 transporter levels.

Even though ABCC2 can functional as a Cry1A toxin receptor 129, specific roles within Bt resistance mechanisms remain elusive. This could involve their direct binding with Bt toxins 23, where it is hypothesized that the extracellular regions of ABC transporters may be sufficient to initiate toxin oligimerization, membrane insertion and pre-pore formation 40, but could also require the involvement one or more additional receptors 33. Alternatively, it was proposed that multiple potential Bt modes of action might exist all capable of forming pre-pores that can insert into the midgut epithelium, such that a number of receptors could be involved. This further suggests that mutations or changes in expression of one or more of these receptors (ALP, APNs, cadherin, and/or ABC transporters) may be capable of influencing Cry toxin resistance traits to different degrees. Indeed, additional research is required in order to decipher the individual roles as well as the interactions between Bt receptors within the framework of toxin modes of actions.

Mutations that cause differences in the level of transcript of candidate Bt resistance genes have been associated with resistant strains of Lepidoptera (see Introduction; Table ). Analogously, structural changes within several of these genes are predicted to cause frameshifts, translation of truncated proteins, or changes in functional amino acids that are proposed to impact resistance. There are 30 SNPs that are fixed differently within midgut expressed transcripts between resistant (ACB-AbR and -AcR) compared to ACB-BtS, and 21 others that are fixed in either ACB-AbR or -AcR. At the functional level, five SNPs lead to predicted changes in the protein sequence and do not exist within any previously identified candidate Bt resistance gene. These results might suggest that structural variation within proteins encoded by Cry1Ab or Cry1Ac resistant O. furnacalis, including candidate Bt resistance genes (ABC transporters, cadherin, alkaline phosphatase, amino peptidases, or proteases), may not contribute to the resistance traits within the colonies used in this study. Any potential linkage of the predicted fixed non-synonymous changes to Bt resistance remains unknown, but will likely be explored within future research.

Conclusions

This study quantitative compared the midgut transcripts among two resistance strain (ACB-AbR and ACB-AcR) and susceptible strain (ACB-BtS). The RNA-seq read depths predicted significant down-regulation of transcripts for previously known Bt resistance genes, aminopeptidase N1 (apn1) and apn3, as well as the abcg gene in both ACB-AbR and -AcR. In contrast, mutations impacting the structure of encoded candidate Bt resistance genes were not detected. These data are important for the understanding of systemic differences between Bt resistant and susceptible genotypes, and modes of adaption to Bt toxins.

Supplementary Data 1.

Click here for additional data file.

Supplementary Data 2.

Click here for additional data file.

Supplementary Data 3.

Click here for additional data file.

Supplementary Data 4.

Click here for additional data file.

Supplementary Data 5.

Click here for additional data file.

Table 2

Summary for the Ostrinia furnacalis midgut transcriptome A) Illumina HiSeq2000 raw and processed read data, and B) the assembled reference sequences.

		Raw read data			Trimmed read data
A) Illumina read data
	Libraries	Count	Length (bp)		Count	Length (bp)
	ACB-AbR.fastq	52,243,872	100		44,502,745	96.4
	ACB-AcR.fastq	47,658,725	100		40,650,855	96.4
	ACB-BtS.fastq	41,167,042	100		35,229,832	96.5
B) Assemblies
	Midgut transcriptome	Count	Length (bp)		Components	Isoforms/comp
	Reference, total1	83,370	617.9±554.2		41,088	2.00±3.45
	Reference, cds2,3	28,940	638.3±423.8		12,288	2.24±4.18

1. Assembly of pooled normalized Illumina read data from ACB-AbR, -AcR, and -BtS libraries. 2. Coding sequences in the reference assembly using Transdecoder (>100 amino acids). 3. Used as template for RNA-seq read mapping and estimation of differential expression.

Table 3

The statistics for the annotation of O. furnacalis midgut transcriptome across annotation resources. NA = no annotation(s) received.

	Database
	BlastX	BlastP	Swiss-Prot	eggNOG	GO	KEGG	Annotated	NA
Unigene Number(%)	28,641(28.64)	21,044(21.04)	32,683(32.68)	12,441(12.44)	2,364(2.36)	3,553(35.54)	39,869(39.87)	60,134(60.13)

Table 4

Estimates of transcripts levels (expression) among candidate Bacillus thuringiensis toxin resistance genes in the midgut of Ostrinia furnacalis based on RNA-seq read mapping to each Trinity component (comp). The transcripts significantly up- (+) or down-regulated (-) in two selected lines (ACB-AbR or -AcR) are based on log2[fold-change] estimates compared to the common susceptible stain (ACB-BtS). NS = not significant.

Trinity component (transcript)	Annotation	ACB-AbR	ACB-AcR	gene
comp3109_c0_seq1:140-646(+)	ABC transporter G1	-2.55(3.1e-21)	-3.54(1.8e-58)	OfABCG1d
comp9932_c0_seq1:509-967(+)	ABC transporter G1	NS	+1.32(5.08e-05)	OfABCG1u
comp15850_c0_seq2:228-680(-)	ABC transporter G1S	+1.00(7.85e-05)	+1.08(5.18e-05)	Ofscarlet
comp17719_c0_seq4:1-1962(+)	ABC transporter G4	-1.52(<0.0001)	-1.48(<0.0001)	OfABCG4
comp19185_c1_seq11:1-807(-)	ABC transporter A2	-1.30(4.12e-42)	NS	OfABCA2
comp19455_c0_seq5:1-5130(+)	ABC transporter A2	-1.26(<0.0001)	NS	OfABCA2
comp16920_c0_seq5:3-1388(-)	ABC transporter A2	-1.19(7.51e-89)	-1.01(1.81e-55)	OfABCA2
comp19185_c1_seq2:783-1433(-)	ABC transporter A2	-1.16(1.06e-33)	-1.04(2.97e-30)	OfABCA2
comp18111_c0_seq1:196-2433(+)	ABC transporter C3	-1.29(30.3e-130)	NS	OfABCC3
comp16077_c0_seq4:2-3469(-)	ABC transporter C3	-1.19(3.13e-162)	NS	OfABCC3
comp19374_c4_seq1:579-2912(+)	Aminopeptidase N1	-1.95(<0.0001)	-1.50(<0.0001)	OfAPN1
comp18815_c0_seq1:194-3178(-)	Aminopeptidase N1	-1.71(<0.0001)	-1.36(<0.0001)	OfAPN1
comp19522_c5_seq3:173-997(-)	Aminopeptidase N3a	-2.95(<0.0001)	-2.29(,0.00001)	OfAPN3a
comp19316_c0_seq9:676-2964(+)	Aminopeptidase N3b	-2.78(<0.0001)	-2.06(<0.00001)	OfAPN3b
comp19316_c0_seq2:479-2614(+)	Aminopeptidase N3b	-2.63(<0.0001)	NS	OfAPN3b
comp19522_c5_seq11:1-474(-)	Aminopeptidase N3c	-3.67(<0.0001)	-2.62(<0.00001)	OfAPN3c
comp19316_c0_seq8:676-1932(+)	Aminopeptidase N3c	-2.73(<0.0001)	-2.11(<0.00001)	OfAPN3c
comp19503_c0_seq3:292-801(-)	Aminopeptidase N8	-2.26(7.0E-87)	+2.07(4.3E-64)	OfAPN8
comp19411_c8_seq1:3-434(-)	chymotrypsin	-4.27(<0.0001)	-4.41(<0.0001)
comp12401_c0_seq1:1-579(-)	trypsin-like SP	-3.60(1.9e-299)	-5.11(<0.0001)
comp14870_c3_seq2:396-728(+)	trypsin alkaline c	-3.50(1.59e-32)	-2.14( <0.0001)
comp16028_c0_seq1:1-447(-)	trypsin-like SP	-3.11(<0.0001)	-3.70(<0.0001)
comp19588_c4_seq4:1-753(+)	trypsin-like SP	-2.34(<0.0001)	-2.46(<0.0001)
comp14501_c0_seq1:93-548	trypsin alkaline c	-2.21(4.57e-14)	NS
comp18273_c0_seq1:295-888(-)	trypsin-like SP	-2.05(<0.0001)	-2.42(<0.0001)
comp19320_c6_seq4:2-748(+)	trypsin-like SP	-2.02(<0.0001)	-2.37(<0.0001)
comp13807_c1_seq3:621-1112(+)	trypsin-like SP	-3.35(5.08e-32)	-5.07(<0.0001)
comp17412_c0_seq2:220-786(+)	trypsin-like SP	-3.27(<0.0001)	-3.80(<0.0001)
comp5178_c0_seq1:3-548(-)	trypsin-like SP	NS	-3.01(3.34e-23)
comp18630_c10_seq1:120-701(-)	trypsin-like SP	NS	-2.76(<0.0001)
comp15635_c1_seq3:174-1067(+)	chymotrypsin	-3.37(<0.0001)	-2.90(1.72e-23)
comp19379_c9_seq4:444-881(+)	chymotrypsin	-2.45(<0.0001)	-2.19(<0.0001)
comp9315_c0_seq2:3-521(+)	chymotrypsin	-2.24(6.84e-22)	NS
comp11509_c0_seq1:286-633(-)	chymotrypsin	+4.29(2.61e-68)	+4.67(3.69e-68)
comp11509_c0_seq2:504-851(-)	chymotrypsin	+3.84(3.87e-63)	+4.51(2.82e-67)
comp5693_c0_seq2:2-433(+)	trypsin-beta-like	+3.77(9.3e_193)	+4.53(4.26e-210)
comp16545_c0_seq2:933-1256(+)	chymotrypsin	+3.35(1.80e-39)	+4.35(1.73e-46)

* Differentially expressed genes have a FDR value of ≤ 0.0001 and the absolute value of log2 (fold change) ≥ 1 or log2 (fold change) ≤ -1 between ACB-AbR or ACB-AcR and ACB-BtS.