Midgut Transcriptome Analysis of Clostera anachoreta Treated with Cry1Ac Toxin

Clostera anachoreta is one of the important Lepidoptera insect pests in forestry, especially in poplars woods in China, Europe, Japan and India, et al, and also the target insect of Cry1Ac toxin and Bt plants. In this study, by using the different dosages of Btcry1Ac toxin to feed larvae and analyzing the transcriptome data, we found six genes, HSC70, GNB2L/RACK1, PNLIP, BI1-like, arylphorin type 2 and PKM, might be associated with Cry1Ac toxin. And PI3K-Akt pathway, which was highly enriched in DEGs and linked to several crucial pathways, including the B cell receptor signaling pathway, toll-like receptor pathway, and MAPK signaling pathway, might be involved in the recovery stage when response to sub-lethal Cry1Ac toxin. This is the first study conducted to specifically investigate C. anachoreta response to Cry toxin stress using large-scale sequencing technologies, and the results highlighted some important genes and pathway that could be involved in Btcry1Ac resistance development or could serve as targets for biologically-based control mechanisms of this insect pest.


Introduction
Insect pest control in forestry and agriculture is increasingly achieved using biological, rather than chemical, insecticides, and the most successful of these is Bacillus thuringiensis (Bt).
However, abuse of Bt pesticides and large-scale planting of Bt plants has caused enormous challenges and threats, including insect resistance and tolerance. Laboratory-selected resistant insect populations have shown that resistance can develop through multiple mechanisms, including alteration of Cry toxin activation (Oppert et al., 1997), toxin sequestration by lipophorin (Ma et al., 2005) or esterases (Gunning et al., 2005), elevation of immune response (Hernández-Martínez et al., 2010), and alteration of toxin receptors, resulting in reduced binding processing and expression, and the metabolic pathway of chronic myelogenous leukemia are highly involved in Bt response and metabolism (Yang Y, 2018).
In this study, the transcriptome data of Clostera anachoreta larvae treated with a continuous sub-lethal dose of Bt Cry1Ac toxin were used to identify Bt-related genes. The data were then combined with growth and development indices to analyze the regulation of related genes involved in digestion, metabolism, and fitness costs, to explore the key genes involved in the response to Cry1Ac treatment and to clarify the mechanisms of the response of target insects to Bt or Bt plants.

Insect rearing
C. anachoreta larvae were collected in wild poplar woods near a highway with no biological control in Baoding, Hebei, China, in May 2016. The larvae were continuously fed without toxins for five generations prior to experimentation. The moths were placed in paper boxes with nets and fed a 10% honey solution. Eggs laid on the mesh were removed and transferred to glass bottles containing mesh and detached leaves. The larvae that hatched from eggs were cultured on the detached leaves, and then used in the experiments. All insect cultures were kept at 26±1°C with 70-80% relative humidity (RH) and a natural photoperiod.

Treatment with Cry1Ac toxin
Third instar larvae were selected for the feeding experiments. The larvae were reared on poplar leaves dipped in different concentrations of Cry1Ac activated toxin solutions (Beijing Meiyan Agricultural Technology), resulting in growth inhibition and increased mortality of the larvae within 5 days. The larvae were treated with three different concentrations of Cry1Ac toxin: 15 μL/mL (Sample A), 7.5 μL/mL (Sample B), and 1.5 μL/mL (Sample C). Larvae treated with 5 Mm/L sodium carbonate solution were used as the negative control treatment. All larval treatments were incubated for 5 days at 26±1°C with 70-80% RH and a natural photoperiod, followed by feeding with control leaves under the same conditions until fourth instar for dissection.
Each treatment group included three parallel samples, with approximately 10-15 larvae per bottle and 10-15 total bottles per treatment. Larvae were starved for 1 day prior to dissection. Dissected midguts were washed with diethyl pyrocarbonate water and blotted dry, then placed into centrifuge tubes (four to five midguts per tube), frozen with liquid nitrogen, and preserved at −80°C until use.

Sample collection and preparation
2.3.1 RNA isolation, quantification and qualification Total RNA was isolated from midgut of C. anachoreta larvae using TRIzol reagent (Invitro-gen,Carlsbad,CA,USA) and RNAprep pure Tissue kit (TIANGEN, DP431) following the manufacturer's protocol. RNA degradation and contamination was monitored on 1% agarose gels.
RNA concentration was measured using Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

Library preparation for Transcriptome sequencing
A total amount of 1.5 µg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) following manufacturer's recommendations and index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer( 5X) . First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase ( RNase H-) . Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3' ends of DNA fragments, NEBNext Adaptor with hairpin loop structure were ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 250-300 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA). Then 3 µl USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37°C for 15 min followed by 5 min at 95 °C before PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. At last, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system.

Clustering and sequencing (Novogene Experimental Department)
The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq platform and paired-end reads were generated.

Data analysis 2.4.1 Quality control
Raw data (raw reads) of fastq format were firstly processed through in-house perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from raw data. At the same time, Q20, Q30, GC-content and sequence duplication level of the clean data were calculated. All the downstream analyses were based on clean data with high quality.

Transcriptome assembly
The left files (read1 files) from all libraries/samples were pooled into one big left.fq file, and right files (read2 files) into one big right.fq file. Transcriptome assembly was accomplished based on the left.fq and right.fq using Trinity  with min_kmer_cov set to 2 by default and all other parameters set default.

Gene functional annotation
Gene function was annotated based on the following databases: Nr (NCBI non-redundant protein sequences); Nt (NCBI non-redundant nucleotide sequences); Pfam (Protein family); KOG/COG (Clusters of Orthologous Groups of proteins); Swiss-Prot (A manually annotated and reviewed protein sequence database); KO (KEGG Ortholog database); GO (Gene Ontology).

Quantification of gene expression levels
Gene expression levels were estimated by RSEM (Li et al, 2011) for each sample: 1. Clean data were mapped back onto the assembled transcriptome. 2. Readcount for each gene was obtained from the mapping results.

Differential expression analysis
For the samples with biological replicates: Differential expression analysis of two conditions/groups was performed using the DESeq R package (1.10.1). DESeq provide statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting P values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted P-value <0.05 found by DESeq were assigned as differentially expressed.
For the samples without biological replicates: Prior to differential gene expression analysis, for each sequenced library, the read counts

Analysis of growth and development indices
Continuous bioassays were conducted using C. anachoreta larvae in the same growth stage.
Larval mortality and the developmental duration of third generation larvae were significantly affected by increased dosages of Cry1Ac. Insect and pupa weights were highest at 7.5 μg/mL(Sample B), and were almost consistent with the control weights ( Figure 1).

Illumina Sequencing and de Novo Assembly
After sequencing using Illumina HiSeqTM 2000 system, we generated 57516756, 56506756, 45779646 and 46263752 raw reads in Sample A, B, C and K (Control) respectively. After removal of low-quality reads and reads containing adaptor or N (It means the bases in reads were not certain), the libraries yielded totals of 54563000, 53586232, 43332730 and 44893556 clean reads.
Their Q30 values (meaning that the base recognition accuracy is 99.9%) were 89.25%, 91.51%, 89.07% and 91.67%, respectively (Table 1) Table 1). The sequencing raw data was not uploaded.  Figure 4). This was consistent with the genetic relationships among these species, supporting the validity of our transcriptome data.
Figure 5 also showed that 205 and 320 unigenes were involved in the cluster for "Defense mechanisms" and "Cell wall/membrane/envelope biogenesis", respectively. And 5304 unigenes were involved in the cluster for Digestion and absorption consisted of "Energy production and conversion", "Lipid transport and metabolism", "Carbohydrate transport and metabolism" and "Amino acid transport and metabolism". In this experiment, the largest group was "General function prediction only" which the same as many other studies ( Tang which might be more important for C. anachoreta responding to BtCry1Ac in midgut. With this in mind, we paid more attention to the signal transduction mechanisms. Furthermore, many of the unigenes participated in defense mechanisms and digestion and absorption, which might imply that the organism would need much more energy to defend foreign substances after dieting Cry1Ac toxin.
In the categories of MF, the largest term is "binding" (GO:0005488) into which 11881 unigenes were classified, followed by the terms "catalytic activity" (15647, GO:0003824). In the BP category, we observed that 201 and 4936 unigenes belonged to the categories "immune system process" and "response to stimulus", respectively, 20493 unigenes were participated in the "metabolic process" (GO:0008152). It was noteworthy that in GO analysis, one unigene could be classified into more than one GO term due to its multiple functions. Hence, there were common unigenes in these three catergories "immune system process", "response to stimulus" and "metabolic process" (Figure 7). Nevertheless, 2913 unigenes (14.21%) in the cluster of "metabolic process" were related to the stress or defense response, only 0.04% and 3.5% unigenes were divided into in the cluster of "immune system process" and "response to stimulus", respectively (Tariq, M. et al., 2015) The midgut chosen for our study contained the certain receptor protein of Cry1Ac, which was dissected and extracted from larvae dieting with Cry1Ac toxin in different doses, might be an important site of resistance or defense response to Cry1Ac. It could explain the difference of the expression between each sample involved in resistance or defense system ( Figure   8).

KEGG analysis
Unigenes were queried against the KEGG database to identify the biological processes in which they are involved, and a total of 14,829 unigenes were assigned to 231 KEGG pathways.

Analysis of differentially expressed genes (DEGs)
In the midgut of C. anachoreta treated with Cry1Ac, the three treatments (A, B, and C)   Based GO enrichment analysis, most DEGs were involved in "Catalytic activity" and "Hydrolytic activity" under the "Molecular Function" category, except for Sample C versus the control; there was no significant difference between the control and the low-dose Cry1Ac treatment (Sample C).

Genes of interest and Bt-related genes in the transcriptome of C. anachoreta
We investigated DEGs potentially involved in insecticidal tolerance in C. anachoreta, including metabolism and Bt-related genes. DEG sequences encoding enzymes functioning in ingestion and detoxification as well as Cry1Ac receptors were extracted from the midgut RNA-Seq data. Cadherin, ABC transporter, ALP, aminopeptidase, P450, GST, carboxylesterase, and serine proteases including trypsin and chymotrypsin (FDR < 0.05, absolute value of log2FC ≥ 2) were significantly differentially expressed (Table 5). Arylphorin subunit beta and arylphorin type 2, related to regeneration of the midgut, also differed significantly between the Cry1Ac treatments and the control. There were many unigenes involved in metabolic systems related to resistance/tolerance or immunity, including GST, P450, CarE, and serine proteases (

Pathways significantly enriched in DEGs
Metabolic pathway enrichment analysis revealed 158 common DEGs involved in 102 pathways (Supplementary Table 7) potentially involved in digestion, absorption, and immunity.
DEGs were significantly enriched in ten functional categories; six categories, Fructose and mannose metabolism, Glycerolipid metabolism, Glycolysis/gluconeogenesis, Carbohydrate digestion and absorption, Pentose and glucuronate interconversions, and Pyruvate metabolism, are involved in digestion and absorption of food and energy metabolism; two categories, Viral myocarditis and Cardiac muscle contraction, are related to cardiac muscle and may affect blood circulation. The two DEGs involved in the Histidine metabolism pathway might influence production of diapause hormone B, thereby affecting insect development duration. Finally, the Tight junction pathway, a Cell-junction pathway, is common in monolayer columnar epithelium and monolayer cubic epithelium in the midgut, and plays an important role in closing the gap at the top of the cell and blocking entry of extraneous macromolecules into the tissue. All DEGs involved in these pathways might by affected by Cry1Ac toxin (P < 0.05) ( Table 8), suggesting that greater energy is needed to sustain life activities when organisms are threatened by drugs or other stressors.  Table 9) with important roles in the response of C. anachoreta to Cry1Ac toxin. These genes were significantly enriched in nine functional categories, four of which, Glutathione metabolism, Drug metabolism-cytochrome P450, Metabolism of xenobiotics by cytochrome P450, and Chemical carcinogenesis, are related to detoxification and immunity (Table 10). Based on the transcriptome data, the response of C. anachoreta to Cry1Ac toxin is related to its tolerance/resistance and immunity to Bt toxin. In this study, the larval midgut transcriptome was sequenced and analyzed by Novogene, and 151,090 high-quality unigenes were obtained, including 1660, 1175, and 562 DEGs between the three larvae treatments with Cry1Ac toxin and the control, respectively. In the P. xylostella midgut transcriptome, 2925 and 2967 unigenes were differently expressed between two resistant strains compared to susceptible strains, respectively (Lei Y, 2014). The DEGs were primarily associated with Bt receptors (e.g., ABC transporter, APN, ALP, and CAD), digestive function (α-amylase, lipase, protease, and carboxypeptidase), immune response (HSC70 and IL22), detoxication enzymes (P450, GST, and also primarily associated with digestion function, with enriched GO terms including monocarboxylic acid metabolic process, pyruvate metabolic process, lipid biosynthetic process, and serine-type endopeptidase activity. We also focused on genes of interest involved in insecticides and Bt resistance or tolerance. These findings may facilitate the understanding of interactions between insects and Bt toxin.

DEGs
There were many DEGs related to growth, development, and metabolism in the midgut of C. are related to carbohydrate degradation and energy metabolism, suggesting that larvae require more energy to stay alive when suffering from stress.
In summary, whether upregulated or downregulated, the DEGs participated a variety of pathways, such as the immune pathway, detoxification metabolism, etc., to regulate cells or life processes, allowing the larvae to remain alive. However, the data are insufficient to determine the specific mechanisms involved in pathways that play crucial roles in the midgut of C. anachoreta after treatment with Cry toxin.

DEGs in pathways of interest
The mechanism of action of Cry toxins has been studied for many years ( Here, the response to Bt Cry1Ac toxin might involve the PI3K-Akt pathway, which was highly enriched in DEGs and linked to several crucial pathways, including the B cell receptor signaling pathway, toll-like receptor pathway, and MAPK signaling pathway. The B cell receptor and toll-like receptor pathways are related to adaptive and innate immunity in the early response to Cry toxin, while the MAPK signaling pathway is involved in the mechanism of resistance in resistant insects. In our study, the PI3K-Akt pathway was activated by a FAK tyrosine kinase (also known as RTK), which was stimulated via the recognition of integrin beta-1 and Von Willebrand factor; the genes of all three were upregulated. FAK is required during development and might regulate tumor suppressor p53. The binding of growth factors to RTK stimulates class Ia PI3K isoform, and PI3K catalyzes the production of phosphatidylinositol-3,4,5-triphosphate (PIP3) at the cell membrane. PIP3 in turn serves as a second messenger that helps to activate Akt. Once active, Akt can control key cellular processes by phosphorylating substrates involved in apoptosis, protein synthesis, metabolism, and cell cycling. In our study, the downregulated AKT gene might have encoded less activated Akt, thereby reducing phosphorylation of downstream factors, and increasing downstream reactions such as BAD or Casp9 activation of the apoptosis process. It could also result in the accumulation of PEPCK, indirectly activating the process of glycolysis, which could generate more ATP (energy) and NADH (reducing power) for other life activities.
Our results differed from those from other resistant insects, possibly due to the timing of midgut dissection, which took place after a period of recovery. During the recovery period until growth to fourth instar, the larvae dieted on control leaves. However, due to damage to the midgut caused by the pore-forming toxin, the larvae required more energy to digest and absorb nutrition from food, and to heal midgut epithelial cells through proliferation and apoptosis.
Therefore, in the midgut of C. anachoreta larvae treated with Bt Cry1Ac and recovered to fourth instar, the glycolysis and oxidative phosphorylation pathways might be activated to generate more ATP (energy) and NADH (reducing power) for larval growth and development.
Meanwhile, apoptosis and regeneration pathways might be stimulated and upregulated for recovery of the larval midgut.