Genome-wide Identification of Tebufenozide Resistant Genes in the smaller tea tortrix, Adoxophyes honmai (Lepidoptera: Tortricidae)

The smaller tea tortrix, Adoxophyes honmai, has developed strong resistance to tebufenozide, a diacylhydrazine-type (DAH) insecticide. Here, we investigated its mechanism by identifying genes responsible for the tebufenozide resistance using various next generation sequencing techniques. First, double-digest restriction site-associated DNA sequencing (ddRAD-seq) identified two candidate loci. Then, synteny analyses using A. honmai draft genome sequences revealed that one locus contained the ecdysone receptor gene (EcR) and the other multiple CYP9A subfamily P450 genes. RNA-seq and direct sequencing of EcR cDNAs found a single nucleotide polymorphism (SNP), which was tightly linked to tebufenozide resistance and generated an amino acid substitution in the ligand-binding domain. The binding affinity to tebufenozide was about 4 times lower in in vitro translated EcR of the resistant strain than in the susceptible strain. RNA-seq analyses identified commonly up-regulated genes in resistant strains, including CYP9A and choline/carboxylesterase (CCE) genes. RT-qPCR analysis and bioassays showed that the expression levels of several CYP9A and CCE genes were moderately correlated with tebufenozide resistance. Collectively, these results suggest that the reduced binding affinity of EcR is the main factor and the enhanced detoxification activity by some CYP9As and CCEs plays a supplementary role in tebufenozide resistance in A. honmai.


Results
Whole-genome sequencing of A. honmai and de novo genome assembly. A draft genome sequence of A. honmai was constructed using the following NGS data on the Kanaya1960-S strain (Supplementary Table S1). Raw sequence data with 45.8 Gb, 36.0 Gb, 36.4 Gb, and 39.7 Gb were obtained from paired-end libraries with insert sizes of 180 bp and 300 bp, and mate-pair libraries with insert sizes of 3 kb and 8 kb, respectively. De novo assembly of these data using the Platanus assmbler 16 Table S2). The number of scaffolds is 72,695 with average length of 6,460 bp, and 35.26 Mb of gaps (7.5% of the total assembly) were contained in the scaffolds. Assessment of the draft genome assembly by BUSCO2 17 using the insecta_odb9 dataset showed that it covers 94.2% of complete BUSCO genes and 3.3% of fragmented BUSCO genes, whereas only 2.5% of BUSCO genes were missing (Supplementary Table S3), indicating a good quality of the draft genome assembly.

Construction of a reference transcriptome assembly for A. honmai. A reference transcriptome
assembly for A. honmai was constructed from pooled RNA sequencing (RNA-seq) data of a susceptible strain (Kanaya1960-S) and a resistant strain (Yui2012-R) (Supplementary Table S1). Basic statistics of this assembly are shown in Supplementary Table S4. In total, 105,001 contigs (Trinity transcripts) were obtained in this assembly and used as reference transcripts throughout the study. These contigs were further clustered into 63,454 genes (Trinity genes). Trinity genes were used for calculating gene expression levels for differentially expressed genes (DEGs) analysis.
Double digest restriction-site associated DNA analysis and construction of a linkage map. The loci responsible for tebufenozide resistance were mapped onto the genome using double digest restriction-site associated DNA (ddRAD-seq) analysis 18 . After cleaning the 555.2 million raw reads (56.0 Gb) obtained from a parental female (Yui2012-R), a parental male (Kanaya1960-S), and 92 F 2 individuals that survived the tebufenozide treatment, 526.7 million clean reads were assigned to each individual and used for further analyses. In total, 20,313 single nucleotide polymorphism (SNP) markers were identified between the parental female and male using Stacks modules 19 , and 945 SNP markers were genotyped on 85 or more F 2 survivors (Supplementary Excel File S1).
Using JoinMap 20 , a linkage map containing 30 linkage groups (LGs) with a total size of 1859.2 cM was generated using 927 of the 945 SNP markers (Supplementary Fig. S1 and Excel File S1). These 927 SNP markers were mapped on 579 unique positions in the linkage map. Although the chromosome number of A. honmai has not yet been determined, the number of LGs was consistent with that of an allied species, Adoxophyes orana (n = 30) 21 . Furthermore, corresponding chromosomes in Bombyx mori were identified for all the LGs by mapping the SNP markers on B. mori chromosomes using the basic local alignment search tool (BLAST) (Supplementary Table S5). Of these, two corresponding A. honmai LGs were identified for three B. mori chromosomes, 11 (LG19 and LG23), 23 (LG6 and LG28), and 24 (LG13 and LG14), which is consistent with previous reports that these three chromosomes were generated by ancestral chromosomal fusion events of lepidopteran species with n = 31 karyotype [22][23][24] . On the other hand, two corresponding B. mori chromosomes, 1 (Z) and 15, were identified for A. honmai LG1 (Supplementary Table S5). The first half region of the LG1 corresponds with B. mori chromosome 15 and the second half region corresponds with B. mori chromosome 1 (Z), which is consistent with a Z-autosome fusion identified in another allied species, Cydia pomonella 25 . One-to-one correspondence was identified between remaining A. honmai LGs and B. mori chromosomes. These results verified that the linkage map of A. honmai was successfully established.
Identification of candidate loci for tebufenozide-resistance by ddRAD-seq. Using this linkage map, we searched for the region responsible for tebufenozide-resistance. Because only F 2 progenies that survived to high-dose tebufenozide exposure were used for this analysis, these survivors were expected to be either resistant homozygotes (R/R) or heterozygotes (R/S) for tebufenozide-resistance genes. Therefore, we expected www.nature.com/scientificreports www.nature.com/scientificreports/ genotypes of SNP markers tightly linked to resistant loci in the F 2 progenies (i.e., resistant-type SNP markers) were either identical to that of the resistant mother (R/R) or heterogenic to that of the parents (R/S). To identify regions where resistant-type SNP markers were significantly enriched, we performed a chi-squared test for the 579 unique SNP markers mapped on the constructed linkage map. We identified SNP markers on 20 LGs with significant deviation from the Mendelian expected ratio (1:2:1) in F 2 survivors (Supplementary Table S6). Of these, the SNP markers 29620 in LG17 and 110923 in LG12 showed significantly high genotype frequency of "R/R or R/S": 100% and 96.7%, respectively. In particular, the genotype frequency of R/R in 29620 was 80.9%, which is much higher than that of marker 110923 (52.94%). We selected the genomic regions encompassing the most significant marker and its neighboring markers in the two LGs as candidate responsible regions for tebufenozide-resistance: 13.83 to 14.96 cM region on LG17 (neighboring markers are 127043 and 48478; Fig. 1) and 18.375 to 18.997 cM region on LG12 (neighboring marker is 14044; Fig. 2).

Synteny analysis on
LG17. To obtain further information on the location of the loci responsible for tebufenozide resistance, we carried out synteny analysis on LG17. A. honmai and B. mori genomes revealed highly conserved synteny between the candidate genomic region on the LG17 of A. honmai (13.83 to 14.96 cM) and the 4.778 to 7.295 Mb region on chromosome 10 of B. mori (Fig. 3). Notably, EcR, which encodes a target molecule of tebufenozide 26 , is located on the middle of this syntenic region in B. mori, suggesting that the EcR of A. honmai (AhEcR) is also located on the corresponding genomic region on LG17 (13.83 to 14.39 cM) and is likely close to marker 29620 (Fig. 3).
Comparison of EcR cDNA sequences between resistant and susceptible strains. Because the candidate responsible region in LG17 seems to contain EcR, we compared EcR cDNA sequences between susceptible and resistant strains based on RNA-seq data. The reference transcriptome assembly mentioned in section "Construction of a reference transcriptome assembly for A. honmai" was used for identifying the EcR cDNA sequence of A. honmai by tblastn search (https://blast.ncbi.nlm.hih.gov/) using the EcR peptide sequence of B. mori as query. A full-length cDNA sequence encoding the EcR B1 isoform in A. honmai (AhEcRB1) was identified. Mapping the clean reads of individual RNA-seq data obtained from two resistant and two susceptible strains to the AhEcRB1 sequence allowed identifying SNPs causing three non-synonymous amino acid substitutions (D26N, A415V, and K546R) in the deduced EcR protein sequence (Table 1). Among them, only SNP C1244T corresponding to A415V is consistent with DAH-resistance phenotypes in the four A. honmai strains, suggesting that this substitution is the most significant.
To validate the linkage of C1244T substitution in AhEcR with DAH resistance, the distribution of C1244T genotypes was compared with that of the SNP markers 29620 and 115778, which were tightly linked to tebufenozide-resistance, using the ddRAD-seq samples ( Supplementary Fig. S2). Accordingly, the genotypes associated with these three SNP markers were exclusively correspondent in 89 F 2 individuals with only one of Adoxyophyes honmai. Genotype frequency indicates the rate of F 2 survivors after tebufenozide treatment with the same genotype to homozygous resistant strain (R/R; diamonds) or with the genotype of homozygous resistant strain or heterozygous strain (R/R or R/S; Squares) in each marker. The ID and position of SNP markers on LG17 are shown at the top. The p-value (circles) indicates the significance of deviation from Mendelian inheritance. The region between the two blue-dotted lines, encompassing the SNP marker 29620, indicates the candidate region responsible for tebufenozide resistance. EcR (bold red font) was estimated to be localized in this region by synteny analysis. See further details in Fig. 3. www.nature.com/scientificreports www.nature.com/scientificreports/ exception (98.8%), further supporting the idea that the SNP C1244T, which engenders A415V substitution in AhEcR protein ( Supplementary Fig. S3), is critical for tebufenozide-resistance.

Functional characterization of the A415 mutation in AhEcR.
To verify the biochemical significance of AhEcR A415V substitution on tebufenozide-resistance, we cloned two AhEcRB1 cDNAs, namely AhEcRB1_S and AhEcRB1_R, from the susceptible strain Kanaya1960-S and the resistant strain Yui2014-R, respectively. The nucleotide sequences of the two cDNAs showed three SNPs, namely T1149C and C1244T, and A1637G; however, because T1149C caused a synonymous substitution, A415V (C1244T) and K546R (A1637G) were the only difference at the protein level. The competition assays using in vitro translated proteins from the AhEcRB1/AhUSP1 heterodimeric complex and [ 3 H]ponasterone A (PonA) revealed that this ecdysteroid, was equally competitive in AhEcRB1_S (pIC 50 = 7.88; IC 50 = 13.1 nM) and AhEcRB1_R (pIC 50 = 7.94; IC 50 = 11.5 nM) (Fig. 4A,B). In contrast, tebufenozide was about five-fold less competitive in AhEcRB1_R (pIC 50 = 7.18; IC 50 = 66.1 nM) than in AhEcRB1_S (pIC 50 = 7.82; IC 50 = 15.1 nM) (Fig. 4C,D). Thus, the binding affinity of PonA was similar between wild and resistant EcR, but the binding affinity of tebufenozide against resistant species was significantly lower (1/4-1/5) than that against wild type. These results revealed that, in the resistant strain, the A415V substitution decreased the affinity to tebufenozide without influencing the affinity to a natural ecdysteroid.
Evaluation of candidate responsible genes in LG12. Next, we searched genes related to tebufenozide resistance in the second candidate responsible region in LG12. Genes of A. honmai on 110 kb genomic region (scaffold8255; 1 b-110 kb) including the significant SNP marker 110923, located at 18.977 cM in LG12, showed apparent conserved synteny with genes of B. mori on 60 kb genomic region (17.51-17.57 Mb) in chromosome 17 (chr17) (Fig. 5). Three P450 genes, namely CYP9A20, CYP9A19, and CYP9A21, were located on the 60 kb genomic region in chr17 of B. mori, while eight homologous CYP9A genes (Supplementary Table S7) were located in the corresponding genomic region of A. honmai. Notably, the SNP marker 110923 was located on an exon of gene CYP9A169 (Fig. 5). These results suggested that one or more CYP9A homologs encoded in scaffold8255 might be relevant to the tebufenozide resistance in the LG12 region.
Analysis of DEGs using RNA-seq data. Generally, the increase of detoxification activity is one of the main causes for insecticide resistance. To identify tebufenozide-metabolyzing enzyme genes that are supposed to be commonly and highly expressed in resistant strains, we performed DEGs analysis. According to the strategy shown in Fig. 6, we compared the gene expression level calculated from RNA-seq data on two tebufenozide resistant strains (Yui2014-R and Met-Sel-R) and two susceptible strains (Kanaya1960-S and Haruno2014-S). First, we compared Yui2014-R to Kanaya1960-S and to Haruno2014-S, and identified 1,648/1,104 and 154/60 up-/ down-regulated DEGs, respectively ( Supplementary Fig. S4). By comparing these DEGs, 114 and 23 common up-and down-regulated DEGs were extracted, respectively (Fig. 6). Similarly, we compared Met-sel-R (LC 50 to tebufenozide = 1,600 ppm) to Kanaya1960-S and to Haruno2014-S, and identified 2,514/4,003 and 688/1,255 up-/down-regulated DEGs, respectively ( Supplementary Fig. S5). Of these, 453 and 870 common up-and down-regulated DEGs were in common between the two comparisons (Fig. 6). Finally, we compared the common DEGs obtained from the two comparisons and identified 44 and 11 overall common up-and down-regulated   www.nature.com/scientificreports www.nature.com/scientificreports/ No other P450, CCE, or GST genes were found in the final common down-regulated DEGs. Notably, CYP9A169 and CYP9A171, likely localized near the tebufenozide-responsible region in LG12 (Fig. 5, Supplementary  Table S7), were found in the overall common up-regulated DEGs. In contrast, no homologous coding genes of B. mori in chromosome 10 were found among these DEGs. These results further support that CYP9A homologs encoded in the candidate responsible region in LG12 might contribute to the tebufenozide resistance.   Table S7) was evaluated by RT-qPCR in 13 regional strains with different magnitude of tebufenozide-resistance (Supplementary Table S8). Each strain showed unique expression profiles in the four CYP9A genes. However, in general, the expression levels of CYP9A genes tended to be higher in the strains with lower corrected mortalities to tebufenozide (Fig. 7). Significant correlations were observed between the corrected mortality and the mRNA expression levels in CYP9A167 and CYP9A168 (r = −0.768 and −0.663, p < 0.05; Spearman correlation) (Supplementary Table S9). Furthermore, a high correlation was observed in the sum of four CYP9A genes (r = −0.7845, p = 0.0015). These results suggested that high expression of CYP9A family genes in the LG12 region is involved in the tebufenozide-resistance in A. honmai.

RT-qPCR analysis of CCEs and GST homologs. We also examined the expression level of AhCCE1-6 and
AhGST1 by RT-qPCR. The expression level of AhCCEs tended to be higher in resistant strains ( Supplementary  Fig. S6). In particular, the expression levels of AhCCE2, AhCCE3, AhCCE4, and the sum of six AhCCE genes   Supplementary Table S8. www.nature.com/scientificreports www.nature.com/scientificreports/ showed significant correlations with the corrected mortality (Supplementary Table S10). By contrast, no significant correlation was found between the expression level of the AhGST1 and the corrected mortality ( Supplementary Fig. S7, Table S10). These results suggested that higher expression of several CCE genes are involved in the tebufenozide resistance in A. honmai, wherease the involvement of GST was unclear.

Discussion
With the aid of a series of NGS analyses, we have conducted genome-wide searching of the genes responsible for tebufenozide resistance in A. honmai and identified candidate genes that implicate two resistant mechanisms. One is the decreased sensitivity to tebufenozide caused by an amino acid substitution in the target molecule EcR, and the other is increased detoxification activity caused by the overexpression of multiple CYP9A and CCE genes.
First, ddRAD-seq analysis identified a candidate locus in LG17. The following synteny analysis between A. honmai and B. mori using the newly-generated draft genome sequences and linkage map of A. honmai revealed that the LG17 corresponds to chromosome 10 of B. mori and the locus contained the EcR gene. Successively, the RNA-seq analysis identified the non-synonymous SNP causing the A415V amino acid substitution in the ligand-binding domain of the EcR protein, and direct cDNA sequencing analysis clarified its tight linkage with tebufenozide resistance. Finally, biochemical analysis using in vitro-translated EcR protein suggested that the affinity to tebufenozide was significantly decreased in the EcR of resistant strains, whereas the affinity to the endogenous ligand, 20E, was not. Collectively, we conclude that the EcR A415V substitution is the major factor accounting for the expression of tebufenozide resistance in A. honmai.
This A415V substitution may reasonably explain the decline of susceptibility to DAH-IGRs without changing the vital function of EcR. The crystal structure of EcR in tobacco budworm, Heliothis virescens, revealed that the amino acid residues critical for binding to a natural, plant-derived ecdysteroid (PonA) are widely conserved among insects 26 . In contrast, different amino acid residues were used for binding to a lepidopteran specific DAH (BYI06830) with similar structure to tebufenozide 27 . The A415 residue in A. honmai EcR (equivalent to A386 residue in H. virescens EcR), neighbours R413, which directly interacts with ecdysteroid in the ligand binding pocket, and V414, which is only conserved among lepidopteran EcRs 28 ( Supplementary Fig. S3). It is thus plausible that A415V substitution modifies the local structure of the ligand binding pocket to selectively alleviate the interaction with DAH-IGRs but not with natural ligands. Further studies on molecular dynamics and crystal structure analysis are required to elucidate the interaction between A415V residues and DAH-IGRs.
The ddRAD-seq analysis identified the second candidate locus in LG12 and the synteny analysis revealed that it corresponds to the region in B. mori chr 17 containing three CYP9A genes. This LG12 region contains at least eight CYP9A genes including two CYP9A19-like (CYP9A164 and CYP9A165), three CYP9A20-like (CYP9A166, CYP9A167, and CYP9A168), and three CYP9A21-like (CYP9A169, CYP9A170, and CYP9A171) genes, suggesting the occurrence of gene duplication in the three CYP9A genes (Fig. 5). Similar extensive duplication events in the CYP9A gene cluster were reported in the three noctuid moths, Spodoptera litura 29 , Spodoptera frugiperda, and Helicoverpa armigera 30 , which have developed high resistance to many insecticides.
The DEG and RT-qPCR analyses further confirmed the involvement of these CYP9A genes in tebufenozide resistance. In general, enzymatic detoxification is one of the major mechanisms for insecticide resistance and P450s play the most important roles. Increased levels of P450 gene expression in resistant insect strains have been reported in numerous cases 31 . In this study, DEG analyses showed that CYP9A166 was commonly upregulated in the resistant strains. Furthermore, RT-qPCR analysis revealed that the expression levels of two CYP9A genes (CYP9A167 and CYP9A168), and sum of four CYP9A family genes encoded in the LG12 region were significantly correlated with the magnitude of resistance (Supplementary Table S9). Altogether, the tebufenozide resistance in A. honmai is, at least in part, attributable to the elevated detoxification activity of the enzymes abundantly produced from the highly expressed CYP9A genes.
In addition to CYP9A genes, the DEG and RT-qPCR analyses indicated that higher expression of CCE genes are also involved in the tebufenozide resistance. The relative importance of each CYP9A and CCE gene on tebufenozide resistance remains unclear. We must note that some field populations of A. honmai in Shizuoka prefecture had acquired strong resistance to various insecticides, such as carbamates, organophosphates, and pyrethroids 14 . In addition, development against diamide insecticides, including flubendiamide and chlorantraniliprole, have been recently reported 32 . The Yui2014-R strain used for the DEG analysis performed in the present study was also resistant to both DAH-IGRs and diamide insecticides (Uchiyama et al., unpublished). P450s and CCEs are known to play a substantial role in the degradation of a number of insecticides, such as carbamates, organophosphates, pyrethroids, and diamides [33][34][35][36] . Although we compared two tebufenozide-resistant strains with two tebufenozide-susceptible strains having different genetic backgrounds to reduce the possibility to pick-up genes unrelated to tebufenozide resistance, it is not completely deniable that some of the highly expressed P450s and CCEs were related to resistance to other insecticides. For further determining the role of each CYP9A and CCE gene in tebufenozide resistance, substrate specificity and catalytic activity of individual enzymes will need to be analysed using recombinant proteins.
The results obtained in the present study are consistent with the mode of inheritance of tebufenozide resistance reported in a previous study: autosomal, incompletely dominant, and controlled by polygenic factors 15 . Both EcR and CYP9As were found to be located in autosomes. It is likely that the decline of binding activity in EcR by the amino acid substitution is a recessive trait, while the increase of detoxification activities by the overexpression of CYP9A and CCE enzymes is a dominant trait. Thus overlapping of these two traits likely appears incomplete dominance.
We believe that the modification of EcR is the main factor while CYP9A overexpression is a supplemental factor, which enhances the former action, based on the following reasons. The correlations between the CYP9A and CCE gene expression levels and the magnitude of tebufenozide-resistance were low to moderate (Supplementary  Table S9 and S10). In addition, the expression levels of CYP9A and CCE genes in some resistant strains (e.g. strains 6,7,8,11, and 12 for CYP9As; strains 4, 5, and 8 for CCEs) were as low as in a susceptible strain (Strain 1-3) (Fig. 7), suggesting that higher expression of CYP9A and CCE genes is dispensable in some resistant strains. In contrast, our preliminary analysis revealed that the frequency of A415V substitution was very highly correlated with the magnitude of tebufenozide-resistance in these populations (r = −0.9049, P = 0.0003, n = 11, Spearman correlation coefficient). Further confirmation of the importance of the EcR substitution in tebufenozide-resistance using more local strains will be reported elsewhere (Uchibori-Asano et al., in preparation).
In conclusion, using NGS techniques, we clarified molecular mechanisms of tebufenozide resistance in A. honmai. The combination of ddRAD-seq, RNA-seq, and genome sequence analyses was highly efficient to reveal tebufenozide resistance mechanisms. The discovery that A415V substitution in EcR causes tebufenozide resistance will be valuable for developing management systems not only for A. honmai but also for other lepidopteran pests.

Methods
Insects. Smaller tea tortrix, A. honmai, were collected from tea fields in different regions in Japan, and maintained in the laboratory on an artificial diet (Insecta LFS, NOSAN) at 25 °C under 16L8D conditions 14 . The strains used for genome sequencing, ddRAD-seq, and RNA-seq analyses are listed in Supplementary Table S1. Strain Kanaya1960-S was the representative tebufenozide-susceptible strain, and it was collected in the 1960s and maintained in the laboratory since then without exposure to any insecticides. The Met-sel-R was established in the laboratory by selection on an artificial diet dipped into methoxyfenozide for 22 generations, as indicated in Supplementary Table S11. Strains used for RT-qPCR of cytochrome P450 genes are listed in Supplementary  Table S8. Bioassays for determining corrected mortality of the strains are described in Supplementary Methods. Whole genome sequencing and de novo genome assembly. Genomic DNA was extracted from 24 adult males of the Kanaya1960-S strain using the DNeasy Blood and Tissue kit (QIAGEN). Because the W chromosome of Lepidoptera usually contains highly repetitive sequences 37 , we avoided extracting DNA from females (W/Z), and only used males (Z/Z) to increase assemble quality. Moreover, two-rounds of inbreeding, starting from a single pair, were carried out to improve genomic DNA homogeneity. Two paired-end (insert size is 180 bp and 500 bp) and two mate-pair (insert size is 3 kb and 8 kb) libraries were generated based on the manufacturer's protocol (https://www.illumina.com/). Each library was sequenced on one lane of the Illumina HiSeq2000, generating 101-bp paired-end and mate-pair reads, respectively (Macrogen Japan Corp.). The paired-end reads and mate-pair reads were cleaned by platanus_trim version 1.0.7 and platanus_internal_trim version 1.0.7 (http://platanus.bio.titech.ac.jp/pltanus_trim), respectively. Platanus version 1.2.4 16 was used for generating a draft genome assembly of A. honmai. First, all the cleaned paired-reads were merged and de novo assembled using the assemble command in Platanus. Second, scaffolding of all the cleaned paired-end and mate-pair reads was performed on the assembled contigs using the scaffold command in Platanus. Finally, gap-closing using all the cleaned reads was performed on the scaffolds using the gap_close command in Platanus, and scaffolds with lengths of 200 bp or more were extracted. The completeness of the generated genome assembly was assessed by BUSCO2 17 using insecta_odb9 (1658 BUSCO genes).
Transcriptome assembly by RNA-seq. Second-instar larvae of Yui2012-R and Kanaya1960-S strains were fed with tea leaves dipped in tebufenozide solution (LC 10 = 46.9 ppm for Yui2012-R and 2.67 ppm for Kanaya1960-S) for 24 h, and surviving larvae were collected. Total RNA was extracted from the larvae using the RNeasy Plus Mini kit (QIAGEN). Two larvae were used in each extraction, which was performed in triplicate. A cDNA library was constructed for each Total RNA sample and then sequenced on the Illumina HiSeq2000 platform by Hokkaido System Science Co. Ltd using the 101-bp paired-end mode. The obtained reads were cleaned using Trimmomatic 0.32 38 . All the clean reads were merged and de novo assembled by Trinity r20140717 39 . Each generated RNA-seq contig was compared with the National Center for Biotechnology Information non-redundant (NCBI-nr) protein database by blastx search (e-value < 1e-3) and top-hit information was extracted. The coding sequence (CDS) was predicted for each RNA-seq contig by TransDecoder bundled with the Trinity r20140717. The predicted CDSs were compared with Pfam domain database by HMMER3 40 and homologous conserved domains were identified. The annotated RNA-seq contigs were used as a reference transcriptome assembly for DEG analysis and linkage analyses. The details of "Analysis of DEGs based on RNA-seq" are described in Supplementary Methods.

ddRAD-seq analysis and construction of linkage groups (LGs). A resistant strain (Yui2012-R) female
and a susceptible strain (Kanaya1960-S) male were crossed, and the resulting F 1 siblings were randomly crossed to obtain F 2 individuals. Second-instar F 2 larvae were fed tea leaves dipped in tebufenozide solution (150 ppm) and surviving larvae (92 individuals) were collected after 10 days (fifth-instar larvae). Genomic DNA was extracted individually from the female and male parents and from the 92 F 2 larvae using the DNA blood and tissue kit (QIAGEN). The ddRAD-seq library was constructed using all genomic DNAs according to the method described in Peterson et al. 18 . Two restriction enzymes, EcoRI and MspI, were used for digesting genomic DNAs. To recognize the genomic DNA of each individual, a unique pair of inline barcodes (ligated onto the EcoRI-associated DNA) and index barcodes (located within sequencing adapters) were assigned to each sample. The ddRAD-seq library was sequenced on two lanes of the Illumina HiSeq2000, generating 101-bp paired-end reads (Macrogen Japan Corp.). Sequence data were processed using Stacks 1.26 modules 19 . Reads obtained for each individual were first demultiplexed and cleaned using the process_radtags module. A catalog of SNP markers was then generated by processing the cleaned reads successively using ustacks (−m = 3 and −M = 3 options), cstacks (−n = 2 option), and sstacks. Finally, F 2 individual genotypes were calculated in the genotypes module. The DNA markers