Evolution of Gustatory Receptor Gene Family Provides Insights into Adaptation to Diverse Host Plants in Nymphalid Butterflies

Evolution of Gustatory Receptor Gene Family Provides Insights into Adaptation to Diverse Host... The host plant range of herbivorous insects is a major aspect of insect–plant interaction, but the genetic basis of host range expansion in insects is poorly understood. In butterflies, gustatory receptor genes (GRs) play important roles in host plant selection by ovipositing females. Since several studies have shown associations between the repertoire sizes of chemosensory gene families and the diversity of resource use, we hypothesized that the increase in the number of genes in the GR family is associated with host range expansion in butterflies. Here, we analyzed the evolutionary dynamics of GRs among related species, including the host generalist Vanessa cardui and three specialists. Although the increase of the GR repertoire itself was not observed, we found that the gene birth rate of GRs was the highest in the lineage leading to V. cardui compared with other specialist lineages. We also identified two taxon-specific subfamilies of GRs, characterized by frequent lineage-specific duplications and higher non-synonymous sub- stitution rates. Together, our results suggest that frequent gene duplications in GRs, which might be involved in the detection of plant secondary metabolites, were associated with host range expansion in the V. cardui lineage. These evolutionary patterns imply that the capability to perceive various compounds during host selection was favored during adaptation to diverse host plants. Key words: gustatory receptor, gene duplication, chemoreception, host range expansion, host plant selection, butterfly. Introduction ancestral species would result in a high speciation rate in Herbivorous insects comprise more than one-third of de- the following lineage. However, the evolutionary processes scribed species, and many authors claim that host plant asso- of host range expansion in herbivorous insects are poorly un- ciations have been major forces driving insect species derstood, because the genetic basis of host plant adaptation diversification (Ehrlich and Raven 1964; Mitter et al. 1988; has been largely unexplored in most insect taxa. Wiens et al. 2015). In particular, the host plant range, which Host plant selection by ovipositing females is one of the key refers to the variety of plant species consumed or used as determinants of host range (Janz and Nylin 1997). Because oviposition sites, is considered an important factor in adapta- the larvae of herbivorous insects generally have low dispersal tion of herbivorous insects. Several studies have suggested abilities, adult females have to recognize and lay eggs on that host range evolution has facilitated species diversification plants that are suitable for larval diet. In butterflies, adult (Janz et al. 2006; Janz and Nylin 2008; Nylin et al. 2014). females can discriminate their host and non-host plants by According to the “oscillation hypothesis” detailed in Janz touching the surface of the plant with their foreleg tarsi and Nylin (2008), expansion of the host plant range in an (Renwick and Chew 1994). This drumming behavior allows The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non- commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com Genome Biol. Evol. 10(6):1351–1362. doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 1351 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Suzuki et al. GBE butterflies to perceive soluble secondary metabolites on the repertoire size of GRs. V. cardui is one of the most polypha- plant, and adult females decide to oviposit or not based on gous butterfly species, with more than 10 plant orders the chemical blend (Nishida et al. 1987; Pereyra and Bowers recorded as its hosts (Robinson et al. 2010). Since it has 1988; Huang et al. 1993). From these observations, chemo- been estimated that host range expansion in V. cardui oc- sensory gene families have been seen as the biggest candidate curred within the genus Vanessa (Nylin et al. 2014), genomic genes regulating host selection behavior of butterflies. analysis among V. cardui and related host specialist species Although the precise genetic mechanisms underlying the provides an opportunity to address relationships between preference for specific hosts for oviposition are yet to be elu- host plant ranges and GRs in butterflies. Specifically, we in- cidated, Ozaki et al. (2011) showed that a gustatory receptor vestigated the numbers of gene gains and losses of GRs gene (GR) encoded a receptor for a host plant-specific com- within a monophyletic species group, structures and clusters pound (i.e., synephrine), and regulated the oviposition behav- of the GR gene phylogeny, and amino acid substitution rates ior of the Asian swallowtail butterfly, Papilio xuthus. across the GR gene family. Moreover, Briscoe et al. (2013) revealed that many GR genes had female-biased expression patterns in the postman butter- Materials and Methods fly, Heliconius melpomene, suggesting the importance of GRs Butterflies and Host Plants for host selection by females. Because drumming behavior is widely observed among butterflies, GRs likely play essential Four butterfly species from the tribe Nymphalini (V. cardui, host selection roles in other butterfly taxa, which also influ- Vanessa indica, Polygonia c-aureum,and Araschnia burejana) ences host range expansion. were used in this study. V. cardui was chosen as the generalist GRs comprise one of the biggest gene families in insect sample because it is the only butterfly species that is both gen- genomes. This gene family evolved under the birth-and-death eralist at the level of plant order (i.e., using more than three model, in which gene repertoires are shaped by multiple gene plant orders, according to Nylin et al. 2014) and available in gains (duplications) and losses (deletions or pseudogeniza- Japan. Other three specialist species are also commonly found tion), and the gene number varies among lineages (Nei in Japan, and each has a different phylogenetic distance from et al. 2008; Sanchez-Gracia et al. 2009). Several studies V. cardui. These specialists rely on a single plant order, the have reported possible associations between the repertoire Rosales, as their hosts. V. cardui was also recorded to utilize sizes of the GR family and the varieties of resource use in Rosales (family Urticaceae), but its primary host plants are from insects. In Drosophila, for example, two specialist species the order Asterales. Butterfly host use data were obtained from have fewer numbers of GRs than their non-specialist sister an online database HOSTS (Robinson et al. 2010). Nylin et al. species because of accelerated gene losses, which were ap- (2014) suggested the host range expansion in the ancestor of parently caused by relaxed selection after specializations on V. cardui occurred after the divergence of the V. cardui lineage chemically homogeneous environments (McBride 2007; and the V. indica lineage, which was around 20 million years McBride et al. 2007). Differences in the number of GRs ago (Wahlberg and Rubinoff 2011). have also been observed in the other direction, that is, host generalists: recent studies on genome sequences of extreme Sampling and RNA Extraction polyphagous herbivores in the family Noctuidae, such as the cotton bollworm moth Helicoverpa armigera and the tobacco Adult females were captured in the wild around Sendai, cutworm Spodoptera litura, revealed that they have remark- Japan, from 2015 to 2016. They were placed in cages with ably larger GR repertoires in the genomes compared with their host plants for oviposition. After oviposition, the eggs other lepidopteran host specialists (Xu et al. 2016; Cheng were collected and reared under constant environmental con- et al. 2017; Gouin et al. 2017; Pearce et al. 2017). These ditions of 25 C, 16L/8D until eclosion. Three female individ- observations raise the possibility that evolutionary transition uals of each species had all their legs dissected within 2 days from host specialist to generalist in butterflies is accompanied of eclosion. The legs of each individual were separately ho- by the increase in the number of GRs by gene duplications. mogenized as soon as possible after dissection, and the total However, the evolutionary relationships between the GRs and RNA was extracted from the homogenates using a Maxwell host expansion remain uncertain, because previous works did 16 LEV Plant RNA kit (Promega Corporation, WI, USA) follow- not present genomes for specialist species in Noctuidae, and ing the manufacturer’s protocol. The qualities of total RNA specialists used for comparisons were all phylogenetically dis- samples were measured using an Agilent 2100 Bioanalyzer tant (e.g., Bombyx mori). Also, the timings of host range ex- (Agilent Technologies, CA, USA). pansion within Noctuidae were not clearly understood. In this article, we analyzed the evolutionary dynamics of RNA-Sequencing and De Novo Assembly GRs in four closely related butterfly species, including the painted lady Vanessa cardui, to test our hypothesis that host Although detection and comparison of the gene repertoires range expansion is associated with the increase in the should be based on whole genome sequences, they were not 1352 Genome Biol. Evol. 10(6):1351–1362 doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Gustatory Receptor Genes and Host Range Expansion GBE available for our study species. Instead, we applied RNA- as the query for TBLASTN to search for contigs that were not sequencing and de novo assembly to identify candidate detected in the first round. genes. As insect GRs generally have very low expression levels The amino acid sequences of the GR genes of our study (Clyne et al. 2000; Ozaki et al. 2011), we focused our se- species (V. cardui, V. indica, P. c-aureum,and A. burejana) quencing efforts on female legs, in which the most diverse and those of five other insect species (H. melpomene, Briscoe GRs are expressed among body parts of the butterfly (Briscoe et al. 2013; D. plexippus, Zhan et al. 2011; P. xuthus, Ozaki et al. et al. 2013), rather than sequencing all tissues evenly at low- 2011; B. mori, Guo et al. 2017;and Drosophila melanogaster, sequencing depths. Our strategy for RNA-seq was to combine Robertson et al. 2003) were aligned using MAFFT 7.273 (Katoh 100 bp and 300 bp paired-end reads obtained using HiSeq and Standley 2013) with the L-INS-i algorithm. An initial gene 2500 and MiSeq sequencers (Illumina, CA, USA), respectively, phylogeny was constructed using the maximum likelihood to obtain the longest possible gene sequences. Briefly, paired- method in RAxML 8.2.8 (Stamatakis 2014), followed by boot- end cDNA libraries were constructed using 500 ng of total strap analyses with 100 replications. We utilized “-m RNA from each butterfly obtained using the TruSeq RNA PROTGAMMAAUTO” option to identify the most suitable pro- Sample Prep Kit v2 (Illumina) according to the manufacturer’s tein substitution model, and JTT model was assigned. protocol. Three libraries were obtained for the HiSeq system At this point, however, we observed many cases where one and one library for the MiSeq system from each species (i.e., GR gene of the nymphalid species seemed to be fragmented RNA samples from one individual of each species was used for into several contigs in the gene phylogeny. To avoid overesti- both HiSeq and MiSeq libraries (see supplementary fig. 1, mation of the number of GRs, we searched for contigs that Supplementary Material online for a schematic explanation). could be integrated into a single gene model. We first identified All RNA-seq runs were conducted at NODAI genome research sets of orthologous genes occupying close positions in the phy- center (Tokyo, Japan) between 2015 and 2016, except for a logeny. Nucleotide sequences of each gene set were then MiSeq library of A. burejana that was sequenced at JT aligned using MAFFT. Here, we set specific criteria to determine Biohistory Research Hall (Osaka, Japan) in 2016. whether the multiple contigs originated from a single gene: 1) If Low-quality raw RNA-seq reads were removed using the multiple contigs of one species aligned to their orthologous fastq_quality_filter of FASTX-toolkit 0.0.13 (http://hannonlab. genes without overlaps between contigs, they were combined cshl.edu/fastx_toolkit/; last accessed May 19, 2018), with a set- into a single gene model. The contig identities were later con- ting of -v -Q 33 -q 20 -p 30. Processed reads were assembled firmed by PCR and electrophoresis (Supplementary Material using Trinity 2.1.1 (Grabherr et al. 2011) with the –trimmo- online). For a few gene models, sequence gaps were deter- matic option on the default settings. All the reads from both mined by Sanger sequencing. 2) If multiple contigs of the spe- HiSeq and MiSeq data were collectively assembled (supple- cies overlapped in alignments, we calculated synonymous mentary fig. 1, Supplementary Material online). Qualities of substitution rates (dS) in the overlapping regions, following de novo assemblies were listed in supplementary table 1, the method used in Duncan et al. (2014).If dS was <0.19, Supplementary Material online. we collapsed the contigs into a single gene model (i.e., contigs were considered as alleles). Otherwise, we counted those con- tigs as paralogs. As dS is often considered as an index for diver- Gene Homology Search, Phylogenetic Analysis, and Gene gence times between homologs (Lynch and Conery 2000), Model Construction what we tried here was collapsing contigs with shorter diver- To detect candidate GR genes from the Trinity assemblies, gence times. The cut-off value for dS (0.19) was average pair- TBLASTN searches (e-value ¼ 1e-05) (https://blast.ncbi.nlm. wise dS between V. cardui and V. indica, calculated from 30 nih.gov/Blast.cgi; last accessed May 19, 2018) were performed randomly chosen BUSCO genes using codeml in PAML v4.8 usingthe aminoacidsequences of GRsof B. mori (Guo et al. (Yang 2007)(supplementary table 2, Supplementary Material 2017), Danaus plexippus (Zhan et al. 2011), and H. melpomene online). Because of this methodology, we had to repeat the (Briscoe et al. 2013) as input queries and our Trinity assemblies phylogeny construction process and gene alignment tests until as databases. Contig profiles selected in TBLASTN were further no contigs remained that could be combined with others. tested with BLASTX (e-value ¼ 1e-02) against the NCBI non- Hereafter, we considered each gene model as a unit of GR gene. redundant protein database (nr). Contigs that matched GRs of After obtaining the final set of GRs, a final gene phylogeny other insect species in BLASTX were considered to be candidate was constructed using the maximum likelihood method in GR genes of the focal species. For contigs having multiple iso- RAxML, with JTT substitution model assigned by forms, we chose the longest contig as the representative se- “PROTGAMMAAUTO” option. Bootstrap analyses were fur- quence. Amino acid sequences of candidate GR contigs were ther carried out with 500 replications. The phylogeny was predicted using TransDecoder (https://transdecoder.github.io; visualized using iTOL v3 (Letunic and Bork 2016). Based on last accessed May 19, 2018) and Sequence Manipulation Suite ligand information of several insect GRs from previous studies, (Stothard 2000). Subsequently, we repeated the same proce- we inferred the ligand affinities of each phylogenetic dure using amino acid sequences of the candidate GR contigs subfamily. Genome Biol. Evol. 10(6):1351–1362 doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 1353 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Suzuki et al. GBE As chemoreception plays a major role in host selection of genes exclusively expressed in other tissues or at other devel- butterflies, we also explored other chemosensory gene fam- opmental stages would have been ignored from the analysis. ilies expressed in female legs. Homology searches and con- In this case, the number of gene gains would be underesti- structions of gene models were also performed for olfactory mated, and the number of losses would be overestimated. receptor (OR), ionotropic receptor (IR), odorant-binding pro- Nevertheless, we believe that the overall pattern of gains and tein (OBP), and chemosensory protein (CSP) gene families in losses was reliable for several reasons. First, legs are one of the thesamemanneras GRs. primary gustatory-sensing organs for butterflies, thus we as- sumed that the numbers of GRs expressed in female legs could be used as substitutes for the numbers of the whole Comparisons of the Number of GRs GR repertoires. In fact, it has been shown that female legs of The difference in the number of GRs (i.e., the number of GR H. melpomene expressed the most diverse set of GRs among gene models) among species was determined. First, we its body parts (50 out of 73 total GRs, 68%) (Briscoe et al. counted the numbers of insect-specific BUSCO (single-copy 2013). The positive relationships between the number of orthologs conserved among the lineage) genes (v3, Insecta) expressed genes at a tissue and the total repertoire size (Simao et al. 2015) in each assembly to consider the variation have been reported in ORs. For example, D. melanogaster in the number of genes due to the de novo assembly proce- expressed 44 out of 62 total ORs (70%), and female H. mel- dure (supplementary table 3, Supplementary Material online). pomene expressed 67 out of 70 total ORs (96%) only in the We calculated the ratio of the number of GRs to that of antennae, the primary olfactory organ of insects (Couto et al. BUSCO in each assembly and tested the heterogeneity of 2005; Briscoe et al. 2013). Although phylogenetically distant the frequency among species using the pairwise Fisher test from insects, similar patterns were found for vertebrate ORs with the Benjamini–Hochberg correction (Benjamini and expressed at olfactory epithelium, the most important olfac- Hochberg 1995). Other chemosensory genes detected were tory organ for mammals: 295 out of 366 ORs (80%) and 817 also analyzed in the same procedure. out of 1,154 ORs (70%) were detected from human and mouse, respectively (Zhang et al. 2004, 2007). In contrast, Estimation of GR Gain/Loss Events and Birth/Death Rates the relationship may not be straightforward for GRs. In S. We estimated the numbers of gains and losses of GRs along litura, the most GR-rich tissue was larval maxilla, but it the species tree using Notung 2.8.1.7 (Stolzer et al. 2012), expressed only 84 out of 237 total GRs (35%) (Cheng et al. which reconciles a gene phylogeny onto a species tree. Input 2017). Adult legs possessed only 20 GRs (8%) (Cheng et al. species tree was manually drawn based on phylogenetic in- 2017). However, these percentages might have been under- formation of species reported in previous studies (Wahlberg estimated, because only 123 GRs were detected in their tran- 2006; Wahlberg et al. 2009, 2013; Wahlberg and Rubinoff scriptome analysis across all tissues of larvae and adults 2011). We used the “–phylogenomics” command imple- (Cheng et al. 2017). Taken together, these observations imply mented in Notung to estimate counts of gains (duplications) that the number of expressed chemosensory genes at an or- and losses (deletions or pseudogenization) of GRs occurring in gan is generally influenced by the size of total gene reper- each phylogenetic branch. We then tested the difference in toires, but the extent to which it reflects the total repertoires the ratio of the number of gains and losses among the most would depend on the importance of chemoreception at the recent branches for each species using the pairwise Fisher test organ. At least, we expect the numbers of expressed GRs in with the Benjamini–Hochberg correction (Benjamini and butterfly legs reflect high proportions of the whole repertoire Hochberg 1995). The same analysis was also performed for sizes, as leg chemoreception is commonly important among OR and IR gene family. butterflies. Second, the gene gain/loss estimation by Notung We then estimated the gene birth rate (ß) and the gene is based on the topology of the input gene phylogeny (Stolzer death rate (@), defined as the number of gene gains or losses et al. 2012). Since we collected complete GR repertoires from per million years per gene, for each phylogenetic branch after several species, our final GR gene phylogeny, in general, was the common ancestor of Nymphalini. Divergence times of highly supported by bootstrap analysis. Therefore, we assume lineages were inferred using BEAST 2.4.7. (Bouckaert et al. that the topology of the phylogeny would not dramatically 2014), based on the nucleotide sequences of five gene change even if missing GRs of Nymphalini were included, and regions (COI, EF-a, wingless, RpS5, and NADPH) of 34 species the overall pattern of gene gains and losses would not be obtained from Genbank (supplementary table 4, supplemen- affected. Further analysis using complete GR gene repertoires tary fig. 2, Supplementary Material online). We calculated the needs to be performed to verify this assumption. gene birth and death rate using the formula developed by Niimura et al. (2014) (supplementary table 5, Supplementary Evolutionary Rate Analysis Material online). However, we admit our analysis of gene gains and losses We took two approaches to infer d /d ratios across the GR N S based on leg transcriptomic data may be inaccurate, because family. First, we chose pairs of orthologous GRs between 1354 Genome Biol. Evol. 10(6):1351–1362 doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Gustatory Receptor Genes and Host Range Expansion GBE Table 1 Number of chemosensory genes (i.e., gene models) and BUSCO genes (v3, Insecta) detected from female leg transcriptomes. Species Vanessa cardui Vanessa indica Polygonia c-aureum Araschnia burejana Host Range Generalist Specialist Specialist Specialist GR 50 28 17 45 OR 25 15 17 24 IR 20 21 27 35 OBP 29 27 32 30 CSP 31 36 33 39 BUSCO 1618 1586 1609 1642 V. cardui and V. indica, and between V. cardui and A. bure- Table 2 jana, based on phylogenetic positions on the GR gene phy- P-values for the Pairwise Fisher’s Test Using Benjamini–Hochberg Correction, Testing Differences in the Ratio of the Number of GRs to logeny. If a GR gene had a one-to-many relationship with its that of BUSCOs among four species orthologs, we considered each different combination as one V. cardui V. indica P. c-aureum pair. In total, we obtained 24 orthologous pairs between V. V. indica 0.04275 — — cardui and V. indica (supplementary table 6, Supplementary P. c-aureum 0.00037 0.12034 — Material online), and 26 pairs between V. cardui and A. bur- A. burejana 0.6034 0.11378 0.00217 ejana (supplementary table 7, Supplementary Material on- line). After nucleotide alignments using MAFFT, we calculated pairwise d /d ratios with codeml in PAML 4.8 N S (Yang 2007). Second, we searched for orthologous gene gene conversion were detected, we re-estimated pairwise groups consisting three species (V. cardui, V. indica,and A. d /d ratios after removing the region of putative conversion burejana as the outgroup), and aligned those nucleotide N S from the alignments. sequences within groups. When one-to-many genes were included in the group, each different combination was con- sidered as an original sample. We obtained 14 orthologous Results groups (OGs) in total (supplementary table 8, Supplementary Detection and Comparisons of the Numbers of GRs Material online). d /d ratios for each branch in the unrooted N S three-species tree were calculated using codeml under branch Although a variation of total assembled bases was observed, model (Yang 2007). The estimation based on branch model is our de novo assemblies exhibited generally similar qualities likely to be more accurate than pairwise, but a problem in our among four species in terms of N50, average contig length, data set was that phylogenetic relationships of Nymphalini and the numbers of BUSCO (supplementary tables 1 and 3, GRs were so complicated it was difficult to find clear ortholog Supplementary Material online). The number of GRs found in groups among three species. Therefore, we integrated pair- V. cardui, V. indica, P. c-aureum,and A. burejana were 50, 27, wise dn/ds estimation to cover across the GR phylogeny. Most 17, and 45, respectively (tables 1 and 2). Most of the anno- GRs in our data set were partial sequences, and consequently, tated GRs were partial sequences. The frequencies of GRs in the analysis was performed only in aligned regions. We then the assembly (i.e., GR/BUSCO ratio) were significantly differ- tested the variation of mean d /d ratios using one-way N S ent between V. cardui and V. indica (pairwise Fisher test; P < ANOVA for phylogenetic subfamilies and using t-test be- 0.05), and between V. cardui and P. c-aureum (pairwise Fisher tween one-to-one and one-to-many pairs, respectively. For test; P < 0.001). However, the difference was not observed three-species OGs, we also tested the difference in selective between V. cardui and A. burejana (pairwise Fisher test; P ¼ pressures between the V. cardui lineage and the V. indica 0.603) (tables 1 and 2). For the other chemosensory gene lineage using paired Wilcoxon test. Any genes showing dS families, we did not observe any significant differences in < 0.01 were excluded from the analysis. the number of genes between V. cardui and other specialists Additionally, we searched for signatures of gene conver- (tables 1 and 2). sion among homologous genes of Nymphalini using GENECONV (Sawyer 1989), as gene conversion can modify the estimation of evolutionary rates. Protein sequences align- GR Gains and Losses along the Nymphalini Phylogeny ments of more than three genes (paralogs or homologs in the It was estimated that losses of GRs occurred more frequently close positions) were created by MAFFT 2.723 (Katoh and than gains in most branches in the Nymphalini lineage (fig. 2 Standley 2013) and used as input data. For genes in which and table 3). In contrast, the most recent branch leading to Genome Biol. Evol. 10(6):1351–1362 doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 1355 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Suzuki et al. GBE Table 3 The Estimated Numbers of Gains and Losses in the GR Family, Along With Gene Birth Rates (b) and Death Rates (@), Among the Most Recent Branches of Four Species. Species Vanessa cardui Vanessa indica Polygonia c-aureum Araschnia burejana Host Range Generalist Specialist Specialist Specialist No. of Gains 7 1 0 3 No. of Losses 9 25 40 24 Gain/Loss ratio 0.778 0.04 0 0.125 Birth rate (b) 0.0536 0.00101 0 0.00099 Death rate (d) 0.00689 0.02519 0.02853 0.00793 Two major taxon-specific GR subfamilies were observed in Table 4 P-values for the Pairwise Fisher’s Test Using Benjamini-Hochberg the gene phylogeny (fig. 2). One comprised nearly half of all Correction, Testing Differences in Gain/Loss Ratios of GRs. the genes, exhibiting frequent lineage-specific expansions. V. cardui V. indica P. c-aureum This subfamily only included GRs of Lepidopteran species (i.e., no Drosophila GRs), and was thus named the V. indica 0.0079 — — P. c-aureum 0.0003 0.4727 — “Lepidoptera-specific” (LS) subfamily. The other smaller sub- A. burejana 0.0487 0.6104 0.0916 family only contained GRs of butterfly species and was des- ignated the “butterfly-specific” (BS) subfamily (fig. 2). This subfamily originated from the fructose receptor clade, but, unlike the sugar and fructose subfamilies, it was characterized by several lineage-specific gene expansions. Putative species- specific gene duplications in the most recent branches of Nymphalini, including seven putatively duplicated genes in V. cardui exhibited a relatively high GR gain/loss ratio (seven the V. cardui lineage, were mapped on either LS or BS sub- gains/nine losses) (fig. 2 and table 3). Comparison of GR gain/ families (fig. 2). loss ratios among the most recent branches for each species revealed that the gain/loss ratio in V. cardui was significantly higher than those of the other species (pairwise Fisher test; P Evolutionary Rate Analysis < 0.01 for V. cardui–V. indica; P < 0.001 for V. cardui–P. c- aureum; P < 0.05 for V. cardui–A. burejana)(table 4). We We estimated pairwise d /d ratios between orthologous N S then estimated the gene birth and death rate at each phylo- genes, and those in lineages of V. cardui and V. indica based genetic branch. The branch leading to V. cardui showed the on the three-species tree, to examine patterns of selective highest birth rate (ß ¼ 0.00536), whereas the death rate was pressures acting on GRs among phylogenetic subfamilies similarto thatofthe branch leading to A. burejana (V. cardui: (i.e., sugar, fructose, CO , LS, and BS). For the pairwise anal- @ ¼ 0.00689, A. burejana: @ ¼ 0.00793) (table 3 and sup- ysis, significant differences in d /d were observed among N S plementary table 5, Supplementary Material online). subfamilies in both combinations of species (one-way However, expansion of the repertoire itself during host range ANOVA; P < 0.05 for V. cardui–V. indica; P < 0.01 for V. expansion was not observed (from 52 to 50 GRs, fig. 1). The cardui–A. burejana)(fig. 3). Particularly, LS and BS subfamilies same analysis for OR and IR families did not indicate any showed higher d /d ratios, whereas sugar, fructose and CO N S 2 lineage-specific accelerations of gene gains or losses (supple- receptor genes exhibited very low evolutionary rates (fig. 3). mentary fig. 3, Supplementary Material online). Moreover, GRs in one-to-many relationships had higher d /d N S ratios compared with one-to-one genes in both combinations (t-test; P < 0.01 for V. cardui and V. indica; P < 0.05 for V. Phylogenetic Analysis cardui and A. burejana)(fig. 3). The variation of d /d esti- N S Ligand information of several GRs could be estimated from mated for lineages of V. cardui and V. indica showed similar the gene phylogeny because they were included in the same patterns in terms of mean values (fig. 3). However, the differ- subfamilies as GRs whose ligands had been reported in past ences in d /d among subfamilies and between one-to-one N S studies, namely, CO (Kwon et al. 2007), sugar (Dahanukar and one-to-many genes were not significant for the V. cardui et al. 2007), and fructose receptors (Sato et al. 2011)(fig. 2). lineage, mainly because of low sample coverage and a few These subfamilies were mainly characterized by one-to-one outlier genes with d /d > 2(fig. 3). Although the mean N S orthologous relationships among distantly related species d /d was higher in the V. cardui lineage than in V. indica, N S (fig. 2). there was no significant difference in evolutionary rates 1356 Genome Biol. Evol. 10(6):1351–1362 doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Gustatory Receptor Genes and Host Range Expansion GBE FIG.1.—Estimation of gene gains and losses in the GR family. The results below the H. melpomene lineage are shown in the figure (see supplementary fig. 3, Supplementary Material online for the results on the whole phylogeny). Nymphalini lineage is boxed in dashed line. Numbers at each tip show the current repertoire size of GRs, and numbers at each node indicate an inferred GR repertoire size of a common ancestor. Estimated numbers of gene gains (þ) and losses () in GRs are shown on branches. The generalist species (V. cardui) is labeled in bold. The branch where host range expansion occurred (according to Nylin et al. 2014) was colored in red. Divergence times were estimated with BEAST 2.4.7. between orthologs of these species (supplementary table 8, Nymphalini (fig. 1). Moreover, the gene birth rate was highest Supplementary Material online). at the V. cardui branch, whereas the death rate was not We detected a possible gene conversion between clearly different from those in the other branches (table 3). VindGR23a and VindGR23b (supplementary table 9, These findings imply that gene duplications in the GR family Supplementary Material online). To examine the effect of occurred frequently in the course of host range expansion, gene conversion, we calculated pairwise d /d ratios for these which is consistent with our hypothesis. However, as Notung N S genes after removing regions of the putative gene conversion. simply considers clusters of conspecific genes in the phylog- As a result, estimated d /d ratios were slightly higher without eny (i.e., genes in one-to-many relationships) to be species- N S conversion than original values (supplementary table 9, specific gene duplication events, it was uncertain whether Supplementary Material online). those one-to-many genes represented actual duplications. Nevertheless, the results showed that one-to-many ortholo- gous pairs had higher d /d ratios than one-to-one pairs N S Discussion (fig. 3), which is consistent with the observation that dupli- The genetic basis underlying the transition from specialist to cated chemoreceptor genes have higher evolutionary rates generalist has been unexplored in herbivorous insects. In but- than non-duplicated homologs (Gardiner et al. 2008; terflies, GRs play important roles in host plant selection by Almeida et al. 2014). Therefore, we expect one-to-many ovipositing females (Ozaki et al. 2011). It has been suggested genes in our data set repersented recent duplication events that the repertoire sizes of this gene family are associated with in the GR family. the diversity of resource use (McBride 2007; McBride et al. In contrast, the expansion of the GR repertoire size itself 2007; Xu et al. 2016; Cheng et al. 2017; Gouin et al. 2017; was not observed along with the increased rate of gene du- Pearce et al. 2017). Thus, we hypothesized that evolutionary plication (fig. 1). Although this result was not consistent with events in the GR family, particularly the increase of the reper- our hypothesis, we note that it was largely influenced by our toire by gene duplications, were associated with adaptation to study design using partial GR repertoires, in which many GRs diverse host plants by butterflies. To test this hypothesis, we were likely to be missing. The algorithm of Notung predicts investigated characteristics of GRs among four closely related the number of genes at each node as the number of OGs species from the tribe Nymphalini, including the generalist V. retained in at least one species below the node (Stolzeretal. cardui and three specialists. 2012). For example, the estimated gene number at the com- mon ancestor of V. cardui and V. indica is the sum of the numbers of OGs found only in V. cardui, those found only Frequent Gene Duplications of GRs Were Associated With in V. indica, and those found in both species. Therefore, if Host Range Expansion partial gene repertoires are placed at the tips, there would be a strong bias that the gene numbers at nodes tend to be We observed that the gain/loss ratio was particularly higher at larger than those at tips, regardless of the actual changes in the branch leading to V. cardui than the other branches within Genome Biol. Evol. 10(6):1351–1362 doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 1357 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Suzuki et al. GBE FIG.2.—Phylogenetic relationships of GRs from nine insect species. The maximum likelihood tree was constructed with RAxML based on amino acid sequences of GRs. Bootstrap analysis was carried out with 500 replicates. Black dots indicate bootstrap support >80%. Subfamilies with putative ligand information are colored in yellow (sugar), orange (fructose), and gray (CO ). Taxon-specific subfamilies are colored in green (Lepidoptera-specific, LS) and light green (Butterfly-specific, BS). Putative species-specific gene duplications are labeled with colored bars. Vcar, V. cardui; Vind, V. indica;Pcau, P. c-aureum; Abur, A. burejana;Hmel, H. melpomene;Dple, D. plexippus;Pxut, P. xuthus;Bmor, B. mori;Dmel, D. melanogaster. the GR gene repertoires. Rather, we would like to focus on proposed as a variable explaining the variation of chemore- the fact that the increase in gene birth rate was observed ceptor gene repertoire sizes among Drosophila (Gardiner et al. despite the bias toward gene loss. Changes in the GR reper- 2008), which could also be affecting nymphalid butterflies. To toire sizes before and after host range expansion should be our knowledge, however, there is no evidence showing that further addressed in whole genome-based studies. genome size is significantly different among our study species, Another concern is that the variation of the GR repertoire thus we could not reach any conclusion about the relatively sizes between specialist and generalist was not entirely con- large GR repertoire for A. burejana. At least, the fact that sistent with previous studies (Xu et al. 2016; Cheng et al. lineages for V. cardui and A. burejana can be clearly separated 2017; Gouin et al. 2017; Pearce et al. 2017). Specifically, by the gene birth rates supported the importance of GRs the number of GRs detected from the generalist V. cardui during host range expansion. was significantly larger than V. indica and P. c-aureum, but Because we used only a single generalist species, it was not different from A. burejana. Again, this could be caused by difficult to identify whether the higher rate of gene duplica- the use of partial GR repertoires: Even if there is a significant tions in GRs was specifically associated with host plant range. difference in the repertoire sizes at the whole genome scale, However, several studies suggest that there could be associ- that difference could be reduced when we extracted only ations between gene duplications and adaptation to diverse genes expressed in female legs. Another explanation is that environments in various organisms. For example, species that A. burejana lineage, the most distant branch from V. cardui, survive under a variety of climatic conditions tend to have was influenced by other factors that maintained its GR reper- higher proportions of duplicated genes in the genomes toire relatively large. Besides host range, genome size was among Drosophila (Makino and Kawata 2012)and mammals 1358 Genome Biol. Evol. 10(6):1351–1362 doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Gustatory Receptor Genes and Host Range Expansion GBE FIG.3.—Variations of d /d ratios across the GR family, estimated in two ways. The upper section shows variations among phylogenetic subfamilies, N S analyzed with one-way ANOVA. The lower section shows comparisons between one-to-one and one-to-many genes, analyzed with t-test. One-to-many genes were not found in the V. indica lineage among ortholog sets for the three-species branch model estimation. (Tamate et al. 2014). In the soil bacteria genus Frankia, (fig. 2), and higher evolutionary rates compared with sugar metabolism-related genes were expanded by gene duplica- and CO clades (fig. 3). Based on these results, we assume tions in a strain infecting diverse plant orders (Normand et al. that genes in LS and BS are involved in the perception of plant 2007). Collectively, these observations raise the possibility that secondary metabolites by nymphalid butterflies. In gene duplications and gene family expansions are common Drosophila, the GR subfamily responding to bitter tastants, processes facilitating adaptation to diverse ecological niches such as plant-derived secondary metabolites, evolved faster across various organisms. Therefore, it would be reasonable than sugar and CO subfamilies (McBride et al. 2007; Weiss to assume that the gene duplication of GRs in the V. cardui et al. 2011; Ling et al. 2014; Delventhal and Carlson 2016; lineage are associated with adaptation to diverse host plants. supplementary fig. 4 and supplementary table 10, However, a more comprehensive study including multiple Supplementary Material online), which resembles the estima- generalist lineages will be necessary to investigate whether tion for LS and BS in our data. Moreover, more than half of the pattern we observed is peculiar to generalist lineages or lost genes during host specialization in Drosophila belonged specific to V. cardui. to the putative bitter receptor clade (McBride et al. 2007; supplementary table 10, Supplementary Material online), which is consistent with the frequent turnover of GR gene Taxon-Specific GR Subfamilies and Perception of Plant repertoires in LS/BS. As secondary metabolite repertoires on Secondary Metabolites plants are likely to be quite variable at both intra- and inter- The phylogenetic analysis identified two taxon-specific GR species scales, divergent evolution of GRs in LS/BS in terms of subfamilies, LS and BS (fig. 2). These subfamilies were char- sequences and gene repertoires might reflect their role as acterized by frequent lineage-specific gene gains and losses, receptors for secondary metabolites. The fact that GRs for including genes duplicated in the recent V. cardui lineage synephrine receptor in P. xuthus (Ozaki et al. 2011) and those Genome Biol. Evol. 10(6):1351–1362 doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 1359 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Suzuki et al. GBE expressed exclusively in female legs in H. melpomene (Briscoe case, duplication of GRs can trigger host range expansion. et al. 2013) were included in BS also imply that the subfamily Under this assumption, the number of GRs would be in- plays some roles in host selection. Although GRs in LS do not creased under selection in species or populations that favor have empirical knowledge to support their functions, three a more extensive host plant range, whereas the number GRs of H. armigera, which could be classified in LS in our gene would be stable as long as specialized hosts are favored and phylogeny, responded to extracts of cotton leaves, suggesting maintained. their functions as bitter taste receptors (Xu et al. 2016). We are unable to determine which scenario is more plau- Overall, these facts imply that GRs for secondary metabolite sible at this point. In fact, the two scenarios of adaptation may receptors of nymphalid butterflies are most likely to exist in have occurred simultaneously in the course of host range ex- the taxon-specific subfamilies. We acknowledge that we can- pansion. To understand the evolutionary relationships be- not ultimately verify their roles without functional analysis, tween the number of GRs and host range expansion in thus candidate bitter GRs should be deorphanized in the fol- detail, we have to answer more questions, such as the timings lowing studies. for duplication and fixation of GRs, ligand affinities and be- havioral influences of each receptor, and distribution of ligands on candidate host plant taxa in nature. Expansion of Detectable Secondary Metabolites May Be an Adaptation to a Wider Host Plant Range Conclusion Our findings together suggest that frequent gene duplica- Although a substantial number of studies have focused on tions in GRs, which might be responsible for secondary oviposition behavior of butterflies, few studies revealed metabolite detection, were associated with host range mechanisms of evolutionary transition between specialists expansion in Nymphalini. Given the variation of ligand af- and generalists. Our analysis of nymphalid butterflies, in- finities among GRs, the increase in the number of GRs is cluding the generalist V. cardui,showed that rapid GR likely to reflect the repertoire expansion of detectable com- gene duplications have occurred during host range expan- pounds during host plant selection by ovipositing females. sion and that there were two taxon-specific GR subfamilies For herbivorous insects, especially butterflies, plant second- that may be related to secondary metabolite detection. ary metabolites are classified either as a stimulant or deter- Together, the results imply that the expansion of detect- rent, which respectively induces or prevents oviposition able secondary metabolites during host selection was an (Renwick and Chew 1994). Thus, the increase in GRs can adaptation to various host plants by generalist butterflies. be interpreted as the repertoire expansion of detectable These findings would establish the foundation for under- stimulants or deterrents, which is possibly an adaptation standing the evolution of plant–insect interactions, and to diverse host plants. how it has facilitated species diversification. If the increase in GRs represents an expansion of detectable deterrents, duplication in GRs would not be the direct cause of the host range expansion because adding new deterrents Supplementary Material into the repertoire would prevent colonization of novel plant Supplementary data are available at Genome Biology and taxa. However, once the host range is expanded by other Evolution online. mechanisms, insects may lay eggs on less suitable species for newly included host plant taxa, because they may not sense the compounds that represent toxicity of the new Acknowledgments taxa. In this situation, the increase in GRs would be favored in generalists because it can enhance the accuracy of host We thank Dr Jun Yokoyama, Dr Nobuaki Nagata, and Dr selection by expanding repertoires of detectable deterrents. Masaya Yago for providing information on butterfly ecology In contrast, GRs would not increase as long as the lineage and distribution; Yurie Hirosaki for assistance on cDNA library remains specialized on specific plant taxa. This concept is con- construction experiments; Michihiko Takahashi and Taisuke sistent with an explanation for herbivorous vertebrates tend- Hori for providing wild butterfly samples; Dr Rebecca P. ing to have larger numbers of bitter taste receptor genes: Duncan for advising methods to discriminate alleles and paral- recognition of various compounds would be adaptive for ac- ogs; and members of the Kawata lab at Tohoku University for curate discrimination of better plants for food (Li and Zhang general advice and discussions. RNA-sequencing experiments 2014). were supported by the Cooperative Research Grant of NODAI On the other hand, it is also possible that the increase in Genome Research Center at Tokyo University of Agriculture. GRs represents an expansion of detectable oviposition stimu- Computational analyses were partially performed on the lants. When butterflies include new compounds into their National Institute of Genetics (NIG) supercomputer at stimulant repertoire, they may start to recognize some plants the Research Organization of Information and Systems of that were previously avoided as oviposition sites. Thus, in this the NIG. 1360 Genome Biol. Evol. 10(6):1351–1362 doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Gustatory Receptor Genes and Host Range Expansion GBE Ling F, Dahanukar A, Weiss LA, Kwon JY, Carlson JR. 2014. The molecular Literature Cited and cellular basis of taste coding in the legs of drosophila. J Neurosci. Almeida FC, Sanchez-Gracia A, Campos JL, Rozas J. 2014. Family size 34(21):7148–7164. evolution in drosophila chemosensory gene families: a comparative Lynch M, Conery JS. 2000. The evolutionary fate and consequences of analysis with a critical appraisal of methods. Genome Biol Evol. duplicate genes. Science 290(5494):1151–1155. 6(7):1669–1682. Makino T, Kawata M. 2012. Habitat variability correlates with duplicate Benjamini Y, Hochberg Y. 1995. Controlling the false discovery rate: a content of drosophila genomes. Mol Biol Evol. 29(10):3169–3179. practical and powerful approach to multiple testing. J R Stat Soc B McBride CS. 2007. Rapid evolution of smell and taste receptor genes 57:289–300. during host specialization in Drosophila sechellia. Proc Natl Acad Sci Bouckaert R, et al. 2014. BEAST 2: a software platform for Bayesian evo- U S A. 104(12):4996–5001. lutionary analysis. PLoS Comput Biol. 10(4):e1003537–e1003536. McBride CS, Arguello JR, O’Meara BC. 2007. Five Drosophila genomes Briscoe AD, et al. 2013. Female behaviour drives expression and evolution reveal nonneutral evolution and the signature of host specialization in of gustatory receptors in butterflies. PLoS Genet. 9(7):e1003620. the chemoreceptor superfamily. Genetics 177(3):1395–1416. Cheng T, et al. 2017. Genomic adaptation to polyphagy and insecticides in Mitter C, Farrell B, Wiegmann B. 1988. The phylogenetic study of adaptive a major East Asian noctuid pest. Nat Ecol Evol. 1(11):1747–1756. zones: has phytophagy promoted insect diversification? Am Nat. Clyne PJ, Warr CG, Carlson JR. 2000. Candidate taste receptors in 132(1):107–128. Drosophila. Science 287(5459):1830–1834. Nei M, Niimura Y, Nozawa M. 2008. The evolution of animal chemo- Couto A, Alenius M, Dickson BJ. 2005. Molecular, anatomical, and func- sensory receptor gene repertoires: roles of chance and necessity. Nat tional organization of the Drosophila olfactory system. Curr Biol. Rev Genet. 9(12):951–963. 15(17):1535–1547. Niimura Y, Matsui A, Touhara K. 2014. Extreme expansion of the olfactory Dahanukar A, Lei YT, Kwon JY, Carlson JR. 2007. Two Gr genes underlie receptor gene repertoire in African elephants and evolutionary dynam- sugar reception in Drosophila. Neuron 56(3):503–516. ics of orthologous gene groups in 13 placental mammals. Genome Delventhal R, Carlson JR. 2016. Bitter taste receptors confer diverse func- Res. 24(9):1485–1496. tions to neurons. Elife 5:1–23. Nishida R, Ohsugi T, Kokubo S, Fukami H. 1987. Oviposition stimulants of Duncan RP, et al. 2014. Dynamic recruitment of amino acid transporters to a Citrus-feeding swallowtail butterfly, Papilio xuthus L. Experientia the insect/symbiont interface. Mol Ecol. 23(6):1608–1623. 43(3):342–344. Ehrlich PR, Raven PH. 1964. Butterflies and plants: a study in coevolution. Normand P, et al. 2006. Genome characteristics of facultatively symbiotic Evolution 18(4):586–608. Frankia sp. strains reflect host range and host plant biogeography. Gardiner A, Barker D, Butlin RK, Jordan WC, Ritchie MG. 2008. Drosophila Genome Res. 17(1):7–15. chemoreceptor gene evolution: selection, specialization and genome Nylin S, Slove J, Janz N. 2014. Host plant utilization, host range oscillations size. Mol Ecol. 17(7):1648–1657. and diversification in nymphalid butterflies: a phylogenetic investiga- Gouin A, et al. 2017. Two genomes of highly polyphagous lepidopteran tion. Evolution 68(1):105–124. pests (Spodoptera frugiperda, Noctuidae) with different host-plant Ozaki K, et al. 2011. A gustatory receptor involved in host plant ranges. Sci Rep. 7(1):1–12. recognition for oviposition of a swallowtail butterfly. Nat Commun. Grabherr MG, et al. 2011. Trinity: reconstructing a full-length transcrip- 2:542. tome without a genome from RNA-Seq data. Nat Biotechnol. Pearce SL, et al. 2017. Genomic innovations, transcriptional plasticity and 29(7):644–652. gene loss underlying the evolution and divergence of two highly po- Guo H, et al. 2017. Expression map of a complete set of gustatory receptor lyphagous and invasive Helicoverpa pest species. BMC Biol. genes in chemosensory organs of Bombyx mori. Insect Biochem Mol 15(1):1–30. Biol. 82:74–82. Pereyra PC, Bowers MD. 1988. Iridoid glycosides as oviposition stimulants Huang X, Renwick JAA, Sachdev-Gupta K. 1993. Oviposition stimulants for the buckeye butterfly, Junonia coenia (Nymphalidae). J Chem Ecol. and deterrents regulating differential acceptance of Iberis amara by 14(3):917–928. Pieris rapae and P. napi oleracea. J Chem Ecol. 19(8):1645–1663. Renwick JAA, Chew FS. 1994. Oviposition behavior in Lepidoptera. Annu Janz N, Nylin S. 1997. The role of female search behaviour in determining Rev Entomol. 39(1):377–400. host plant range in plant feeding insects: a test of the information Robertson HM, Warr CG, Carlson JR. 2003. Molecular evolution of the processing hypothesis. Proc R Soc B Biol Sci. 264(1382):701–707. insect chemoreceptor gene superfamily in Drosophila melanogaster. Janz N, Nylin S. 2008. The oscillation hypothesis of host-plant range and Proc Natl Acad Sci U S A. 100(Suppl 2):14537–14542. speciation. Specialization, Speciation, and Radiation: The Evolutionary Robinson GS, Ackery PR, Kitching IJ, Beccaloni GW, Hern andez LM. 2010. Biology of Herbivorous Insects. University of California Press. p. 203– HOSTS - A database of the world’s Lepidopteran Hostplants. London: Natural History Museum. Janz N, Nylin S, Wahlberg N. 2006. Diversity begets diversity: host expan- Sanchez-Gracia A, Vieira FG, Rozas J. 2009. Molecular evolution of the sions and the diversification of plant-feeding insects. BMC Evol Biol. major chemosensory gene families in insects. Heredity 6:4. 103(3):208–216. Katoh K, Standley DM. 2013. MAFFT multiple sequence alignment soft- Sato K, Tanaka K, Touhara K. 2011. Sugar-regulated cation channel ware version 7: improvements in performance and usability. Mol Biol formed by an insect gustatory receptor. Proc Natl Acad Sci U S A. Evol. 30(4):772–780. 108(28):11680–11685. Kwon JY, Dahanukar A, Weiss LA, Carlson JR. 2007. The molecular basis Sawyer SA. 1989. Statistical tests for detecting gene conversion. Mol Biol of CO2 reception in Drosophila. Proc Natl Acad Sci U S A. Evol. 6(5):526–538. 104(9):3574–3578. Simao FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. Letunic I, Bork P. 2016. Interactive tree of life (iTOL) v3: an online tool for 2015. BUSCO: assessing genome assembly and annotation the display and annotation of phylogenetic and other trees. Nucleic completeness with single-copy orthologs. Bioinformatics Acids Res. 44(W1):W242–W245. 31(19):3210–3212. Li D, Zhang J. 2014. Diet shapes the evolution of the vertebrate bitter taste Stamatakis A. 2014. RAxML version 8: A tool for phylogenetic analysis and receptor gene repertoire. Mol Biol Evol. 31(2):303–309. post-analysis of large phylogenies. Bioinformatics 30(9):1312–1313. Genome Biol. Evol. 10(6):1351–1362 doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 1361 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Suzuki et al. GBE Stolzer M, et al. 2012. Inferring duplications, losses, transfers and incom- Weiss LA, Dahanukar A, Kwon JY, Banerjee D, Carlson JR. 2011. The plete lineage sorting with nonbinary species trees. Bioinformatics molecular and cellular basis of bitter taste in Drosophila. Neuron 28(18):i409–i415. 69(2):258–272. Stothard P. 2000. The sequence manipulation suite: javaScript programs Wiens JJ, Lapoint RT, Whiteman NK. 2015. Herbivory increases diversifi- for analyzing and formatting protein and DNA sequences. cation across insect clades. Nat Commun. 6(1):8370. Biotechniques 28:1102–1104. Xu W, Papanicolaou A, Zhang H-J, Anderson A. 2016. Expansion of a Tamate SC, Kawata M, Makino T. 2014. Contribution of nonohnologous bitter taste receptor family in a polyphagous insect herbivore. Sci duplicated genes to high habitat variability in mammals. Mol Biol Evol. Rep. 6(1):23666. 31(7):1779–1786. Yang Z. 2007. PAML 4: phylogenetic analysis by maximum likelihood. Mol Wahlberg N. 2006. That Awkward Age for Butterflies: insights from the Biol Evol. 24(8):1586–1591. Age of the Butterfly Subfamily Nymphalinae (Lepidoptera: nymphali- Zhan S, Merlin C, Boore JL, Reppert SM. 2011. The monarch butterfly dae). Syst Biol. 55(5):703–714. genome yields insights into long-distance migration. Cell Wahlberg N, et al. 2009. Nymphalid butterflies diversify following near 147(5):1171–1185. demise at the Cretaceous/Tertiary boundary. Proc R Soc B Biol Sci. Zhang X, et al. 2007. Characterizing the expression of the human olfactory 276(1677):4295–4302. receptor gene family using a novel DNA microarray. Genome Biol. Wahlberg N, Rubinoff D. 2011. Vagility across Vanessa (Lepidoptera: nym- 8(5):r86. phalidae): Mobility in butterfly species does not inhibit the formation Zhang X, et al. 2004. High-throughput microarray detection of olfactory and persistence of isolated sister taxa. Syst Entomol. 36(2):362–370. receptor gene expression in the mouse. Proc Natl Acad Sci U S A. Wahlberg N, Wheat CW, Pen ~a C. 2013. Timing and patterns in the tax- 101(39):14168–14173. onomic diversification of Lepidoptera (butterflies and moths). PLoS One 8(11):e80875–e80878. Associate editor: Belinda Chang 1362 Genome Biol. Evol. 10(6):1351–1362 doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Genome Biology and Evolution Oxford University Press

Evolution of Gustatory Receptor Gene Family Provides Insights into Adaptation to Diverse Host Plants in Nymphalid Butterflies

Free
12 pages

Loading next page...
 
/lp/ou_press/evolution-of-gustatory-receptor-gene-family-provides-insights-into-o129YpOJkk
Publisher
Oxford University Press
Copyright
© The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution.
ISSN
1759-6653
eISSN
1759-6653
D.O.I.
10.1093/gbe/evy093
Publisher site
See Article on Publisher Site

Abstract

The host plant range of herbivorous insects is a major aspect of insect–plant interaction, but the genetic basis of host range expansion in insects is poorly understood. In butterflies, gustatory receptor genes (GRs) play important roles in host plant selection by ovipositing females. Since several studies have shown associations between the repertoire sizes of chemosensory gene families and the diversity of resource use, we hypothesized that the increase in the number of genes in the GR family is associated with host range expansion in butterflies. Here, we analyzed the evolutionary dynamics of GRs among related species, including the host generalist Vanessa cardui and three specialists. Although the increase of the GR repertoire itself was not observed, we found that the gene birth rate of GRs was the highest in the lineage leading to V. cardui compared with other specialist lineages. We also identified two taxon-specific subfamilies of GRs, characterized by frequent lineage-specific duplications and higher non-synonymous sub- stitution rates. Together, our results suggest that frequent gene duplications in GRs, which might be involved in the detection of plant secondary metabolites, were associated with host range expansion in the V. cardui lineage. These evolutionary patterns imply that the capability to perceive various compounds during host selection was favored during adaptation to diverse host plants. Key words: gustatory receptor, gene duplication, chemoreception, host range expansion, host plant selection, butterfly. Introduction ancestral species would result in a high speciation rate in Herbivorous insects comprise more than one-third of de- the following lineage. However, the evolutionary processes scribed species, and many authors claim that host plant asso- of host range expansion in herbivorous insects are poorly un- ciations have been major forces driving insect species derstood, because the genetic basis of host plant adaptation diversification (Ehrlich and Raven 1964; Mitter et al. 1988; has been largely unexplored in most insect taxa. Wiens et al. 2015). In particular, the host plant range, which Host plant selection by ovipositing females is one of the key refers to the variety of plant species consumed or used as determinants of host range (Janz and Nylin 1997). Because oviposition sites, is considered an important factor in adapta- the larvae of herbivorous insects generally have low dispersal tion of herbivorous insects. Several studies have suggested abilities, adult females have to recognize and lay eggs on that host range evolution has facilitated species diversification plants that are suitable for larval diet. In butterflies, adult (Janz et al. 2006; Janz and Nylin 2008; Nylin et al. 2014). females can discriminate their host and non-host plants by According to the “oscillation hypothesis” detailed in Janz touching the surface of the plant with their foreleg tarsi and Nylin (2008), expansion of the host plant range in an (Renwick and Chew 1994). This drumming behavior allows The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non- commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com Genome Biol. Evol. 10(6):1351–1362. doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 1351 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Suzuki et al. GBE butterflies to perceive soluble secondary metabolites on the repertoire size of GRs. V. cardui is one of the most polypha- plant, and adult females decide to oviposit or not based on gous butterfly species, with more than 10 plant orders the chemical blend (Nishida et al. 1987; Pereyra and Bowers recorded as its hosts (Robinson et al. 2010). Since it has 1988; Huang et al. 1993). From these observations, chemo- been estimated that host range expansion in V. cardui oc- sensory gene families have been seen as the biggest candidate curred within the genus Vanessa (Nylin et al. 2014), genomic genes regulating host selection behavior of butterflies. analysis among V. cardui and related host specialist species Although the precise genetic mechanisms underlying the provides an opportunity to address relationships between preference for specific hosts for oviposition are yet to be elu- host plant ranges and GRs in butterflies. Specifically, we in- cidated, Ozaki et al. (2011) showed that a gustatory receptor vestigated the numbers of gene gains and losses of GRs gene (GR) encoded a receptor for a host plant-specific com- within a monophyletic species group, structures and clusters pound (i.e., synephrine), and regulated the oviposition behav- of the GR gene phylogeny, and amino acid substitution rates ior of the Asian swallowtail butterfly, Papilio xuthus. across the GR gene family. Moreover, Briscoe et al. (2013) revealed that many GR genes had female-biased expression patterns in the postman butter- Materials and Methods fly, Heliconius melpomene, suggesting the importance of GRs Butterflies and Host Plants for host selection by females. Because drumming behavior is widely observed among butterflies, GRs likely play essential Four butterfly species from the tribe Nymphalini (V. cardui, host selection roles in other butterfly taxa, which also influ- Vanessa indica, Polygonia c-aureum,and Araschnia burejana) ences host range expansion. were used in this study. V. cardui was chosen as the generalist GRs comprise one of the biggest gene families in insect sample because it is the only butterfly species that is both gen- genomes. This gene family evolved under the birth-and-death eralist at the level of plant order (i.e., using more than three model, in which gene repertoires are shaped by multiple gene plant orders, according to Nylin et al. 2014) and available in gains (duplications) and losses (deletions or pseudogeniza- Japan. Other three specialist species are also commonly found tion), and the gene number varies among lineages (Nei in Japan, and each has a different phylogenetic distance from et al. 2008; Sanchez-Gracia et al. 2009). Several studies V. cardui. These specialists rely on a single plant order, the have reported possible associations between the repertoire Rosales, as their hosts. V. cardui was also recorded to utilize sizes of the GR family and the varieties of resource use in Rosales (family Urticaceae), but its primary host plants are from insects. In Drosophila, for example, two specialist species the order Asterales. Butterfly host use data were obtained from have fewer numbers of GRs than their non-specialist sister an online database HOSTS (Robinson et al. 2010). Nylin et al. species because of accelerated gene losses, which were ap- (2014) suggested the host range expansion in the ancestor of parently caused by relaxed selection after specializations on V. cardui occurred after the divergence of the V. cardui lineage chemically homogeneous environments (McBride 2007; and the V. indica lineage, which was around 20 million years McBride et al. 2007). Differences in the number of GRs ago (Wahlberg and Rubinoff 2011). have also been observed in the other direction, that is, host generalists: recent studies on genome sequences of extreme Sampling and RNA Extraction polyphagous herbivores in the family Noctuidae, such as the cotton bollworm moth Helicoverpa armigera and the tobacco Adult females were captured in the wild around Sendai, cutworm Spodoptera litura, revealed that they have remark- Japan, from 2015 to 2016. They were placed in cages with ably larger GR repertoires in the genomes compared with their host plants for oviposition. After oviposition, the eggs other lepidopteran host specialists (Xu et al. 2016; Cheng were collected and reared under constant environmental con- et al. 2017; Gouin et al. 2017; Pearce et al. 2017). These ditions of 25 C, 16L/8D until eclosion. Three female individ- observations raise the possibility that evolutionary transition uals of each species had all their legs dissected within 2 days from host specialist to generalist in butterflies is accompanied of eclosion. The legs of each individual were separately ho- by the increase in the number of GRs by gene duplications. mogenized as soon as possible after dissection, and the total However, the evolutionary relationships between the GRs and RNA was extracted from the homogenates using a Maxwell host expansion remain uncertain, because previous works did 16 LEV Plant RNA kit (Promega Corporation, WI, USA) follow- not present genomes for specialist species in Noctuidae, and ing the manufacturer’s protocol. The qualities of total RNA specialists used for comparisons were all phylogenetically dis- samples were measured using an Agilent 2100 Bioanalyzer tant (e.g., Bombyx mori). Also, the timings of host range ex- (Agilent Technologies, CA, USA). pansion within Noctuidae were not clearly understood. In this article, we analyzed the evolutionary dynamics of RNA-Sequencing and De Novo Assembly GRs in four closely related butterfly species, including the painted lady Vanessa cardui, to test our hypothesis that host Although detection and comparison of the gene repertoires range expansion is associated with the increase in the should be based on whole genome sequences, they were not 1352 Genome Biol. Evol. 10(6):1351–1362 doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Gustatory Receptor Genes and Host Range Expansion GBE available for our study species. Instead, we applied RNA- as the query for TBLASTN to search for contigs that were not sequencing and de novo assembly to identify candidate detected in the first round. genes. As insect GRs generally have very low expression levels The amino acid sequences of the GR genes of our study (Clyne et al. 2000; Ozaki et al. 2011), we focused our se- species (V. cardui, V. indica, P. c-aureum,and A. burejana) quencing efforts on female legs, in which the most diverse and those of five other insect species (H. melpomene, Briscoe GRs are expressed among body parts of the butterfly (Briscoe et al. 2013; D. plexippus, Zhan et al. 2011; P. xuthus, Ozaki et al. et al. 2013), rather than sequencing all tissues evenly at low- 2011; B. mori, Guo et al. 2017;and Drosophila melanogaster, sequencing depths. Our strategy for RNA-seq was to combine Robertson et al. 2003) were aligned using MAFFT 7.273 (Katoh 100 bp and 300 bp paired-end reads obtained using HiSeq and Standley 2013) with the L-INS-i algorithm. An initial gene 2500 and MiSeq sequencers (Illumina, CA, USA), respectively, phylogeny was constructed using the maximum likelihood to obtain the longest possible gene sequences. Briefly, paired- method in RAxML 8.2.8 (Stamatakis 2014), followed by boot- end cDNA libraries were constructed using 500 ng of total strap analyses with 100 replications. We utilized “-m RNA from each butterfly obtained using the TruSeq RNA PROTGAMMAAUTO” option to identify the most suitable pro- Sample Prep Kit v2 (Illumina) according to the manufacturer’s tein substitution model, and JTT model was assigned. protocol. Three libraries were obtained for the HiSeq system At this point, however, we observed many cases where one and one library for the MiSeq system from each species (i.e., GR gene of the nymphalid species seemed to be fragmented RNA samples from one individual of each species was used for into several contigs in the gene phylogeny. To avoid overesti- both HiSeq and MiSeq libraries (see supplementary fig. 1, mation of the number of GRs, we searched for contigs that Supplementary Material online for a schematic explanation). could be integrated into a single gene model. We first identified All RNA-seq runs were conducted at NODAI genome research sets of orthologous genes occupying close positions in the phy- center (Tokyo, Japan) between 2015 and 2016, except for a logeny. Nucleotide sequences of each gene set were then MiSeq library of A. burejana that was sequenced at JT aligned using MAFFT. Here, we set specific criteria to determine Biohistory Research Hall (Osaka, Japan) in 2016. whether the multiple contigs originated from a single gene: 1) If Low-quality raw RNA-seq reads were removed using the multiple contigs of one species aligned to their orthologous fastq_quality_filter of FASTX-toolkit 0.0.13 (http://hannonlab. genes without overlaps between contigs, they were combined cshl.edu/fastx_toolkit/; last accessed May 19, 2018), with a set- into a single gene model. The contig identities were later con- ting of -v -Q 33 -q 20 -p 30. Processed reads were assembled firmed by PCR and electrophoresis (Supplementary Material using Trinity 2.1.1 (Grabherr et al. 2011) with the –trimmo- online). For a few gene models, sequence gaps were deter- matic option on the default settings. All the reads from both mined by Sanger sequencing. 2) If multiple contigs of the spe- HiSeq and MiSeq data were collectively assembled (supple- cies overlapped in alignments, we calculated synonymous mentary fig. 1, Supplementary Material online). Qualities of substitution rates (dS) in the overlapping regions, following de novo assemblies were listed in supplementary table 1, the method used in Duncan et al. (2014).If dS was <0.19, Supplementary Material online. we collapsed the contigs into a single gene model (i.e., contigs were considered as alleles). Otherwise, we counted those con- tigs as paralogs. As dS is often considered as an index for diver- Gene Homology Search, Phylogenetic Analysis, and Gene gence times between homologs (Lynch and Conery 2000), Model Construction what we tried here was collapsing contigs with shorter diver- To detect candidate GR genes from the Trinity assemblies, gence times. The cut-off value for dS (0.19) was average pair- TBLASTN searches (e-value ¼ 1e-05) (https://blast.ncbi.nlm. wise dS between V. cardui and V. indica, calculated from 30 nih.gov/Blast.cgi; last accessed May 19, 2018) were performed randomly chosen BUSCO genes using codeml in PAML v4.8 usingthe aminoacidsequences of GRsof B. mori (Guo et al. (Yang 2007)(supplementary table 2, Supplementary Material 2017), Danaus plexippus (Zhan et al. 2011), and H. melpomene online). Because of this methodology, we had to repeat the (Briscoe et al. 2013) as input queries and our Trinity assemblies phylogeny construction process and gene alignment tests until as databases. Contig profiles selected in TBLASTN were further no contigs remained that could be combined with others. tested with BLASTX (e-value ¼ 1e-02) against the NCBI non- Hereafter, we considered each gene model as a unit of GR gene. redundant protein database (nr). Contigs that matched GRs of After obtaining the final set of GRs, a final gene phylogeny other insect species in BLASTX were considered to be candidate was constructed using the maximum likelihood method in GR genes of the focal species. For contigs having multiple iso- RAxML, with JTT substitution model assigned by forms, we chose the longest contig as the representative se- “PROTGAMMAAUTO” option. Bootstrap analyses were fur- quence. Amino acid sequences of candidate GR contigs were ther carried out with 500 replications. The phylogeny was predicted using TransDecoder (https://transdecoder.github.io; visualized using iTOL v3 (Letunic and Bork 2016). Based on last accessed May 19, 2018) and Sequence Manipulation Suite ligand information of several insect GRs from previous studies, (Stothard 2000). Subsequently, we repeated the same proce- we inferred the ligand affinities of each phylogenetic dure using amino acid sequences of the candidate GR contigs subfamily. Genome Biol. Evol. 10(6):1351–1362 doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 1353 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Suzuki et al. GBE As chemoreception plays a major role in host selection of genes exclusively expressed in other tissues or at other devel- butterflies, we also explored other chemosensory gene fam- opmental stages would have been ignored from the analysis. ilies expressed in female legs. Homology searches and con- In this case, the number of gene gains would be underesti- structions of gene models were also performed for olfactory mated, and the number of losses would be overestimated. receptor (OR), ionotropic receptor (IR), odorant-binding pro- Nevertheless, we believe that the overall pattern of gains and tein (OBP), and chemosensory protein (CSP) gene families in losses was reliable for several reasons. First, legs are one of the thesamemanneras GRs. primary gustatory-sensing organs for butterflies, thus we as- sumed that the numbers of GRs expressed in female legs could be used as substitutes for the numbers of the whole Comparisons of the Number of GRs GR repertoires. In fact, it has been shown that female legs of The difference in the number of GRs (i.e., the number of GR H. melpomene expressed the most diverse set of GRs among gene models) among species was determined. First, we its body parts (50 out of 73 total GRs, 68%) (Briscoe et al. counted the numbers of insect-specific BUSCO (single-copy 2013). The positive relationships between the number of orthologs conserved among the lineage) genes (v3, Insecta) expressed genes at a tissue and the total repertoire size (Simao et al. 2015) in each assembly to consider the variation have been reported in ORs. For example, D. melanogaster in the number of genes due to the de novo assembly proce- expressed 44 out of 62 total ORs (70%), and female H. mel- dure (supplementary table 3, Supplementary Material online). pomene expressed 67 out of 70 total ORs (96%) only in the We calculated the ratio of the number of GRs to that of antennae, the primary olfactory organ of insects (Couto et al. BUSCO in each assembly and tested the heterogeneity of 2005; Briscoe et al. 2013). Although phylogenetically distant the frequency among species using the pairwise Fisher test from insects, similar patterns were found for vertebrate ORs with the Benjamini–Hochberg correction (Benjamini and expressed at olfactory epithelium, the most important olfac- Hochberg 1995). Other chemosensory genes detected were tory organ for mammals: 295 out of 366 ORs (80%) and 817 also analyzed in the same procedure. out of 1,154 ORs (70%) were detected from human and mouse, respectively (Zhang et al. 2004, 2007). In contrast, Estimation of GR Gain/Loss Events and Birth/Death Rates the relationship may not be straightforward for GRs. In S. We estimated the numbers of gains and losses of GRs along litura, the most GR-rich tissue was larval maxilla, but it the species tree using Notung 2.8.1.7 (Stolzer et al. 2012), expressed only 84 out of 237 total GRs (35%) (Cheng et al. which reconciles a gene phylogeny onto a species tree. Input 2017). Adult legs possessed only 20 GRs (8%) (Cheng et al. species tree was manually drawn based on phylogenetic in- 2017). However, these percentages might have been under- formation of species reported in previous studies (Wahlberg estimated, because only 123 GRs were detected in their tran- 2006; Wahlberg et al. 2009, 2013; Wahlberg and Rubinoff scriptome analysis across all tissues of larvae and adults 2011). We used the “–phylogenomics” command imple- (Cheng et al. 2017). Taken together, these observations imply mented in Notung to estimate counts of gains (duplications) that the number of expressed chemosensory genes at an or- and losses (deletions or pseudogenization) of GRs occurring in gan is generally influenced by the size of total gene reper- each phylogenetic branch. We then tested the difference in toires, but the extent to which it reflects the total repertoires the ratio of the number of gains and losses among the most would depend on the importance of chemoreception at the recent branches for each species using the pairwise Fisher test organ. At least, we expect the numbers of expressed GRs in with the Benjamini–Hochberg correction (Benjamini and butterfly legs reflect high proportions of the whole repertoire Hochberg 1995). The same analysis was also performed for sizes, as leg chemoreception is commonly important among OR and IR gene family. butterflies. Second, the gene gain/loss estimation by Notung We then estimated the gene birth rate (ß) and the gene is based on the topology of the input gene phylogeny (Stolzer death rate (@), defined as the number of gene gains or losses et al. 2012). Since we collected complete GR repertoires from per million years per gene, for each phylogenetic branch after several species, our final GR gene phylogeny, in general, was the common ancestor of Nymphalini. Divergence times of highly supported by bootstrap analysis. Therefore, we assume lineages were inferred using BEAST 2.4.7. (Bouckaert et al. that the topology of the phylogeny would not dramatically 2014), based on the nucleotide sequences of five gene change even if missing GRs of Nymphalini were included, and regions (COI, EF-a, wingless, RpS5, and NADPH) of 34 species the overall pattern of gene gains and losses would not be obtained from Genbank (supplementary table 4, supplemen- affected. Further analysis using complete GR gene repertoires tary fig. 2, Supplementary Material online). We calculated the needs to be performed to verify this assumption. gene birth and death rate using the formula developed by Niimura et al. (2014) (supplementary table 5, Supplementary Evolutionary Rate Analysis Material online). However, we admit our analysis of gene gains and losses We took two approaches to infer d /d ratios across the GR N S based on leg transcriptomic data may be inaccurate, because family. First, we chose pairs of orthologous GRs between 1354 Genome Biol. Evol. 10(6):1351–1362 doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Gustatory Receptor Genes and Host Range Expansion GBE Table 1 Number of chemosensory genes (i.e., gene models) and BUSCO genes (v3, Insecta) detected from female leg transcriptomes. Species Vanessa cardui Vanessa indica Polygonia c-aureum Araschnia burejana Host Range Generalist Specialist Specialist Specialist GR 50 28 17 45 OR 25 15 17 24 IR 20 21 27 35 OBP 29 27 32 30 CSP 31 36 33 39 BUSCO 1618 1586 1609 1642 V. cardui and V. indica, and between V. cardui and A. bure- Table 2 jana, based on phylogenetic positions on the GR gene phy- P-values for the Pairwise Fisher’s Test Using Benjamini–Hochberg Correction, Testing Differences in the Ratio of the Number of GRs to logeny. If a GR gene had a one-to-many relationship with its that of BUSCOs among four species orthologs, we considered each different combination as one V. cardui V. indica P. c-aureum pair. In total, we obtained 24 orthologous pairs between V. V. indica 0.04275 — — cardui and V. indica (supplementary table 6, Supplementary P. c-aureum 0.00037 0.12034 — Material online), and 26 pairs between V. cardui and A. bur- A. burejana 0.6034 0.11378 0.00217 ejana (supplementary table 7, Supplementary Material on- line). After nucleotide alignments using MAFFT, we calculated pairwise d /d ratios with codeml in PAML 4.8 N S (Yang 2007). Second, we searched for orthologous gene gene conversion were detected, we re-estimated pairwise groups consisting three species (V. cardui, V. indica,and A. d /d ratios after removing the region of putative conversion burejana as the outgroup), and aligned those nucleotide N S from the alignments. sequences within groups. When one-to-many genes were included in the group, each different combination was con- sidered as an original sample. We obtained 14 orthologous Results groups (OGs) in total (supplementary table 8, Supplementary Detection and Comparisons of the Numbers of GRs Material online). d /d ratios for each branch in the unrooted N S three-species tree were calculated using codeml under branch Although a variation of total assembled bases was observed, model (Yang 2007). The estimation based on branch model is our de novo assemblies exhibited generally similar qualities likely to be more accurate than pairwise, but a problem in our among four species in terms of N50, average contig length, data set was that phylogenetic relationships of Nymphalini and the numbers of BUSCO (supplementary tables 1 and 3, GRs were so complicated it was difficult to find clear ortholog Supplementary Material online). The number of GRs found in groups among three species. Therefore, we integrated pair- V. cardui, V. indica, P. c-aureum,and A. burejana were 50, 27, wise dn/ds estimation to cover across the GR phylogeny. Most 17, and 45, respectively (tables 1 and 2). Most of the anno- GRs in our data set were partial sequences, and consequently, tated GRs were partial sequences. The frequencies of GRs in the analysis was performed only in aligned regions. We then the assembly (i.e., GR/BUSCO ratio) were significantly differ- tested the variation of mean d /d ratios using one-way N S ent between V. cardui and V. indica (pairwise Fisher test; P < ANOVA for phylogenetic subfamilies and using t-test be- 0.05), and between V. cardui and P. c-aureum (pairwise Fisher tween one-to-one and one-to-many pairs, respectively. For test; P < 0.001). However, the difference was not observed three-species OGs, we also tested the difference in selective between V. cardui and A. burejana (pairwise Fisher test; P ¼ pressures between the V. cardui lineage and the V. indica 0.603) (tables 1 and 2). For the other chemosensory gene lineage using paired Wilcoxon test. Any genes showing dS families, we did not observe any significant differences in < 0.01 were excluded from the analysis. the number of genes between V. cardui and other specialists Additionally, we searched for signatures of gene conver- (tables 1 and 2). sion among homologous genes of Nymphalini using GENECONV (Sawyer 1989), as gene conversion can modify the estimation of evolutionary rates. Protein sequences align- GR Gains and Losses along the Nymphalini Phylogeny ments of more than three genes (paralogs or homologs in the It was estimated that losses of GRs occurred more frequently close positions) were created by MAFFT 2.723 (Katoh and than gains in most branches in the Nymphalini lineage (fig. 2 Standley 2013) and used as input data. For genes in which and table 3). In contrast, the most recent branch leading to Genome Biol. Evol. 10(6):1351–1362 doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 1355 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Suzuki et al. GBE Table 3 The Estimated Numbers of Gains and Losses in the GR Family, Along With Gene Birth Rates (b) and Death Rates (@), Among the Most Recent Branches of Four Species. Species Vanessa cardui Vanessa indica Polygonia c-aureum Araschnia burejana Host Range Generalist Specialist Specialist Specialist No. of Gains 7 1 0 3 No. of Losses 9 25 40 24 Gain/Loss ratio 0.778 0.04 0 0.125 Birth rate (b) 0.0536 0.00101 0 0.00099 Death rate (d) 0.00689 0.02519 0.02853 0.00793 Two major taxon-specific GR subfamilies were observed in Table 4 P-values for the Pairwise Fisher’s Test Using Benjamini-Hochberg the gene phylogeny (fig. 2). One comprised nearly half of all Correction, Testing Differences in Gain/Loss Ratios of GRs. the genes, exhibiting frequent lineage-specific expansions. V. cardui V. indica P. c-aureum This subfamily only included GRs of Lepidopteran species (i.e., no Drosophila GRs), and was thus named the V. indica 0.0079 — — P. c-aureum 0.0003 0.4727 — “Lepidoptera-specific” (LS) subfamily. The other smaller sub- A. burejana 0.0487 0.6104 0.0916 family only contained GRs of butterfly species and was des- ignated the “butterfly-specific” (BS) subfamily (fig. 2). This subfamily originated from the fructose receptor clade, but, unlike the sugar and fructose subfamilies, it was characterized by several lineage-specific gene expansions. Putative species- specific gene duplications in the most recent branches of Nymphalini, including seven putatively duplicated genes in V. cardui exhibited a relatively high GR gain/loss ratio (seven the V. cardui lineage, were mapped on either LS or BS sub- gains/nine losses) (fig. 2 and table 3). Comparison of GR gain/ families (fig. 2). loss ratios among the most recent branches for each species revealed that the gain/loss ratio in V. cardui was significantly higher than those of the other species (pairwise Fisher test; P Evolutionary Rate Analysis < 0.01 for V. cardui–V. indica; P < 0.001 for V. cardui–P. c- aureum; P < 0.05 for V. cardui–A. burejana)(table 4). We We estimated pairwise d /d ratios between orthologous N S then estimated the gene birth and death rate at each phylo- genes, and those in lineages of V. cardui and V. indica based genetic branch. The branch leading to V. cardui showed the on the three-species tree, to examine patterns of selective highest birth rate (ß ¼ 0.00536), whereas the death rate was pressures acting on GRs among phylogenetic subfamilies similarto thatofthe branch leading to A. burejana (V. cardui: (i.e., sugar, fructose, CO , LS, and BS). For the pairwise anal- @ ¼ 0.00689, A. burejana: @ ¼ 0.00793) (table 3 and sup- ysis, significant differences in d /d were observed among N S plementary table 5, Supplementary Material online). subfamilies in both combinations of species (one-way However, expansion of the repertoire itself during host range ANOVA; P < 0.05 for V. cardui–V. indica; P < 0.01 for V. expansion was not observed (from 52 to 50 GRs, fig. 1). The cardui–A. burejana)(fig. 3). Particularly, LS and BS subfamilies same analysis for OR and IR families did not indicate any showed higher d /d ratios, whereas sugar, fructose and CO N S 2 lineage-specific accelerations of gene gains or losses (supple- receptor genes exhibited very low evolutionary rates (fig. 3). mentary fig. 3, Supplementary Material online). Moreover, GRs in one-to-many relationships had higher d /d N S ratios compared with one-to-one genes in both combinations (t-test; P < 0.01 for V. cardui and V. indica; P < 0.05 for V. Phylogenetic Analysis cardui and A. burejana)(fig. 3). The variation of d /d esti- N S Ligand information of several GRs could be estimated from mated for lineages of V. cardui and V. indica showed similar the gene phylogeny because they were included in the same patterns in terms of mean values (fig. 3). However, the differ- subfamilies as GRs whose ligands had been reported in past ences in d /d among subfamilies and between one-to-one N S studies, namely, CO (Kwon et al. 2007), sugar (Dahanukar and one-to-many genes were not significant for the V. cardui et al. 2007), and fructose receptors (Sato et al. 2011)(fig. 2). lineage, mainly because of low sample coverage and a few These subfamilies were mainly characterized by one-to-one outlier genes with d /d > 2(fig. 3). Although the mean N S orthologous relationships among distantly related species d /d was higher in the V. cardui lineage than in V. indica, N S (fig. 2). there was no significant difference in evolutionary rates 1356 Genome Biol. Evol. 10(6):1351–1362 doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Gustatory Receptor Genes and Host Range Expansion GBE FIG.1.—Estimation of gene gains and losses in the GR family. The results below the H. melpomene lineage are shown in the figure (see supplementary fig. 3, Supplementary Material online for the results on the whole phylogeny). Nymphalini lineage is boxed in dashed line. Numbers at each tip show the current repertoire size of GRs, and numbers at each node indicate an inferred GR repertoire size of a common ancestor. Estimated numbers of gene gains (þ) and losses () in GRs are shown on branches. The generalist species (V. cardui) is labeled in bold. The branch where host range expansion occurred (according to Nylin et al. 2014) was colored in red. Divergence times were estimated with BEAST 2.4.7. between orthologs of these species (supplementary table 8, Nymphalini (fig. 1). Moreover, the gene birth rate was highest Supplementary Material online). at the V. cardui branch, whereas the death rate was not We detected a possible gene conversion between clearly different from those in the other branches (table 3). VindGR23a and VindGR23b (supplementary table 9, These findings imply that gene duplications in the GR family Supplementary Material online). To examine the effect of occurred frequently in the course of host range expansion, gene conversion, we calculated pairwise d /d ratios for these which is consistent with our hypothesis. However, as Notung N S genes after removing regions of the putative gene conversion. simply considers clusters of conspecific genes in the phylog- As a result, estimated d /d ratios were slightly higher without eny (i.e., genes in one-to-many relationships) to be species- N S conversion than original values (supplementary table 9, specific gene duplication events, it was uncertain whether Supplementary Material online). those one-to-many genes represented actual duplications. Nevertheless, the results showed that one-to-many ortholo- gous pairs had higher d /d ratios than one-to-one pairs N S Discussion (fig. 3), which is consistent with the observation that dupli- The genetic basis underlying the transition from specialist to cated chemoreceptor genes have higher evolutionary rates generalist has been unexplored in herbivorous insects. In but- than non-duplicated homologs (Gardiner et al. 2008; terflies, GRs play important roles in host plant selection by Almeida et al. 2014). Therefore, we expect one-to-many ovipositing females (Ozaki et al. 2011). It has been suggested genes in our data set repersented recent duplication events that the repertoire sizes of this gene family are associated with in the GR family. the diversity of resource use (McBride 2007; McBride et al. In contrast, the expansion of the GR repertoire size itself 2007; Xu et al. 2016; Cheng et al. 2017; Gouin et al. 2017; was not observed along with the increased rate of gene du- Pearce et al. 2017). Thus, we hypothesized that evolutionary plication (fig. 1). Although this result was not consistent with events in the GR family, particularly the increase of the reper- our hypothesis, we note that it was largely influenced by our toire by gene duplications, were associated with adaptation to study design using partial GR repertoires, in which many GRs diverse host plants by butterflies. To test this hypothesis, we were likely to be missing. The algorithm of Notung predicts investigated characteristics of GRs among four closely related the number of genes at each node as the number of OGs species from the tribe Nymphalini, including the generalist V. retained in at least one species below the node (Stolzeretal. cardui and three specialists. 2012). For example, the estimated gene number at the com- mon ancestor of V. cardui and V. indica is the sum of the numbers of OGs found only in V. cardui, those found only Frequent Gene Duplications of GRs Were Associated With in V. indica, and those found in both species. Therefore, if Host Range Expansion partial gene repertoires are placed at the tips, there would be a strong bias that the gene numbers at nodes tend to be We observed that the gain/loss ratio was particularly higher at larger than those at tips, regardless of the actual changes in the branch leading to V. cardui than the other branches within Genome Biol. Evol. 10(6):1351–1362 doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 1357 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Suzuki et al. GBE FIG.2.—Phylogenetic relationships of GRs from nine insect species. The maximum likelihood tree was constructed with RAxML based on amino acid sequences of GRs. Bootstrap analysis was carried out with 500 replicates. Black dots indicate bootstrap support >80%. Subfamilies with putative ligand information are colored in yellow (sugar), orange (fructose), and gray (CO ). Taxon-specific subfamilies are colored in green (Lepidoptera-specific, LS) and light green (Butterfly-specific, BS). Putative species-specific gene duplications are labeled with colored bars. Vcar, V. cardui; Vind, V. indica;Pcau, P. c-aureum; Abur, A. burejana;Hmel, H. melpomene;Dple, D. plexippus;Pxut, P. xuthus;Bmor, B. mori;Dmel, D. melanogaster. the GR gene repertoires. Rather, we would like to focus on proposed as a variable explaining the variation of chemore- the fact that the increase in gene birth rate was observed ceptor gene repertoire sizes among Drosophila (Gardiner et al. despite the bias toward gene loss. Changes in the GR reper- 2008), which could also be affecting nymphalid butterflies. To toire sizes before and after host range expansion should be our knowledge, however, there is no evidence showing that further addressed in whole genome-based studies. genome size is significantly different among our study species, Another concern is that the variation of the GR repertoire thus we could not reach any conclusion about the relatively sizes between specialist and generalist was not entirely con- large GR repertoire for A. burejana. At least, the fact that sistent with previous studies (Xu et al. 2016; Cheng et al. lineages for V. cardui and A. burejana can be clearly separated 2017; Gouin et al. 2017; Pearce et al. 2017). Specifically, by the gene birth rates supported the importance of GRs the number of GRs detected from the generalist V. cardui during host range expansion. was significantly larger than V. indica and P. c-aureum, but Because we used only a single generalist species, it was not different from A. burejana. Again, this could be caused by difficult to identify whether the higher rate of gene duplica- the use of partial GR repertoires: Even if there is a significant tions in GRs was specifically associated with host plant range. difference in the repertoire sizes at the whole genome scale, However, several studies suggest that there could be associ- that difference could be reduced when we extracted only ations between gene duplications and adaptation to diverse genes expressed in female legs. Another explanation is that environments in various organisms. For example, species that A. burejana lineage, the most distant branch from V. cardui, survive under a variety of climatic conditions tend to have was influenced by other factors that maintained its GR reper- higher proportions of duplicated genes in the genomes toire relatively large. Besides host range, genome size was among Drosophila (Makino and Kawata 2012)and mammals 1358 Genome Biol. Evol. 10(6):1351–1362 doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Gustatory Receptor Genes and Host Range Expansion GBE FIG.3.—Variations of d /d ratios across the GR family, estimated in two ways. The upper section shows variations among phylogenetic subfamilies, N S analyzed with one-way ANOVA. The lower section shows comparisons between one-to-one and one-to-many genes, analyzed with t-test. One-to-many genes were not found in the V. indica lineage among ortholog sets for the three-species branch model estimation. (Tamate et al. 2014). In the soil bacteria genus Frankia, (fig. 2), and higher evolutionary rates compared with sugar metabolism-related genes were expanded by gene duplica- and CO clades (fig. 3). Based on these results, we assume tions in a strain infecting diverse plant orders (Normand et al. that genes in LS and BS are involved in the perception of plant 2007). Collectively, these observations raise the possibility that secondary metabolites by nymphalid butterflies. In gene duplications and gene family expansions are common Drosophila, the GR subfamily responding to bitter tastants, processes facilitating adaptation to diverse ecological niches such as plant-derived secondary metabolites, evolved faster across various organisms. Therefore, it would be reasonable than sugar and CO subfamilies (McBride et al. 2007; Weiss to assume that the gene duplication of GRs in the V. cardui et al. 2011; Ling et al. 2014; Delventhal and Carlson 2016; lineage are associated with adaptation to diverse host plants. supplementary fig. 4 and supplementary table 10, However, a more comprehensive study including multiple Supplementary Material online), which resembles the estima- generalist lineages will be necessary to investigate whether tion for LS and BS in our data. Moreover, more than half of the pattern we observed is peculiar to generalist lineages or lost genes during host specialization in Drosophila belonged specific to V. cardui. to the putative bitter receptor clade (McBride et al. 2007; supplementary table 10, Supplementary Material online), which is consistent with the frequent turnover of GR gene Taxon-Specific GR Subfamilies and Perception of Plant repertoires in LS/BS. As secondary metabolite repertoires on Secondary Metabolites plants are likely to be quite variable at both intra- and inter- The phylogenetic analysis identified two taxon-specific GR species scales, divergent evolution of GRs in LS/BS in terms of subfamilies, LS and BS (fig. 2). These subfamilies were char- sequences and gene repertoires might reflect their role as acterized by frequent lineage-specific gene gains and losses, receptors for secondary metabolites. The fact that GRs for including genes duplicated in the recent V. cardui lineage synephrine receptor in P. xuthus (Ozaki et al. 2011) and those Genome Biol. Evol. 10(6):1351–1362 doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 1359 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Suzuki et al. GBE expressed exclusively in female legs in H. melpomene (Briscoe case, duplication of GRs can trigger host range expansion. et al. 2013) were included in BS also imply that the subfamily Under this assumption, the number of GRs would be in- plays some roles in host selection. Although GRs in LS do not creased under selection in species or populations that favor have empirical knowledge to support their functions, three a more extensive host plant range, whereas the number GRs of H. armigera, which could be classified in LS in our gene would be stable as long as specialized hosts are favored and phylogeny, responded to extracts of cotton leaves, suggesting maintained. their functions as bitter taste receptors (Xu et al. 2016). We are unable to determine which scenario is more plau- Overall, these facts imply that GRs for secondary metabolite sible at this point. In fact, the two scenarios of adaptation may receptors of nymphalid butterflies are most likely to exist in have occurred simultaneously in the course of host range ex- the taxon-specific subfamilies. We acknowledge that we can- pansion. To understand the evolutionary relationships be- not ultimately verify their roles without functional analysis, tween the number of GRs and host range expansion in thus candidate bitter GRs should be deorphanized in the fol- detail, we have to answer more questions, such as the timings lowing studies. for duplication and fixation of GRs, ligand affinities and be- havioral influences of each receptor, and distribution of ligands on candidate host plant taxa in nature. Expansion of Detectable Secondary Metabolites May Be an Adaptation to a Wider Host Plant Range Conclusion Our findings together suggest that frequent gene duplica- Although a substantial number of studies have focused on tions in GRs, which might be responsible for secondary oviposition behavior of butterflies, few studies revealed metabolite detection, were associated with host range mechanisms of evolutionary transition between specialists expansion in Nymphalini. Given the variation of ligand af- and generalists. Our analysis of nymphalid butterflies, in- finities among GRs, the increase in the number of GRs is cluding the generalist V. cardui,showed that rapid GR likely to reflect the repertoire expansion of detectable com- gene duplications have occurred during host range expan- pounds during host plant selection by ovipositing females. sion and that there were two taxon-specific GR subfamilies For herbivorous insects, especially butterflies, plant second- that may be related to secondary metabolite detection. ary metabolites are classified either as a stimulant or deter- Together, the results imply that the expansion of detect- rent, which respectively induces or prevents oviposition able secondary metabolites during host selection was an (Renwick and Chew 1994). Thus, the increase in GRs can adaptation to various host plants by generalist butterflies. be interpreted as the repertoire expansion of detectable These findings would establish the foundation for under- stimulants or deterrents, which is possibly an adaptation standing the evolution of plant–insect interactions, and to diverse host plants. how it has facilitated species diversification. If the increase in GRs represents an expansion of detectable deterrents, duplication in GRs would not be the direct cause of the host range expansion because adding new deterrents Supplementary Material into the repertoire would prevent colonization of novel plant Supplementary data are available at Genome Biology and taxa. However, once the host range is expanded by other Evolution online. mechanisms, insects may lay eggs on less suitable species for newly included host plant taxa, because they may not sense the compounds that represent toxicity of the new Acknowledgments taxa. In this situation, the increase in GRs would be favored in generalists because it can enhance the accuracy of host We thank Dr Jun Yokoyama, Dr Nobuaki Nagata, and Dr selection by expanding repertoires of detectable deterrents. Masaya Yago for providing information on butterfly ecology In contrast, GRs would not increase as long as the lineage and distribution; Yurie Hirosaki for assistance on cDNA library remains specialized on specific plant taxa. This concept is con- construction experiments; Michihiko Takahashi and Taisuke sistent with an explanation for herbivorous vertebrates tend- Hori for providing wild butterfly samples; Dr Rebecca P. ing to have larger numbers of bitter taste receptor genes: Duncan for advising methods to discriminate alleles and paral- recognition of various compounds would be adaptive for ac- ogs; and members of the Kawata lab at Tohoku University for curate discrimination of better plants for food (Li and Zhang general advice and discussions. RNA-sequencing experiments 2014). were supported by the Cooperative Research Grant of NODAI On the other hand, it is also possible that the increase in Genome Research Center at Tokyo University of Agriculture. GRs represents an expansion of detectable oviposition stimu- Computational analyses were partially performed on the lants. When butterflies include new compounds into their National Institute of Genetics (NIG) supercomputer at stimulant repertoire, they may start to recognize some plants the Research Organization of Information and Systems of that were previously avoided as oviposition sites. Thus, in this the NIG. 1360 Genome Biol. Evol. 10(6):1351–1362 doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Gustatory Receptor Genes and Host Range Expansion GBE Ling F, Dahanukar A, Weiss LA, Kwon JY, Carlson JR. 2014. The molecular Literature Cited and cellular basis of taste coding in the legs of drosophila. J Neurosci. Almeida FC, Sanchez-Gracia A, Campos JL, Rozas J. 2014. Family size 34(21):7148–7164. evolution in drosophila chemosensory gene families: a comparative Lynch M, Conery JS. 2000. The evolutionary fate and consequences of analysis with a critical appraisal of methods. Genome Biol Evol. duplicate genes. Science 290(5494):1151–1155. 6(7):1669–1682. Makino T, Kawata M. 2012. Habitat variability correlates with duplicate Benjamini Y, Hochberg Y. 1995. Controlling the false discovery rate: a content of drosophila genomes. Mol Biol Evol. 29(10):3169–3179. practical and powerful approach to multiple testing. J R Stat Soc B McBride CS. 2007. Rapid evolution of smell and taste receptor genes 57:289–300. during host specialization in Drosophila sechellia. Proc Natl Acad Sci Bouckaert R, et al. 2014. BEAST 2: a software platform for Bayesian evo- U S A. 104(12):4996–5001. lutionary analysis. PLoS Comput Biol. 10(4):e1003537–e1003536. McBride CS, Arguello JR, O’Meara BC. 2007. Five Drosophila genomes Briscoe AD, et al. 2013. Female behaviour drives expression and evolution reveal nonneutral evolution and the signature of host specialization in of gustatory receptors in butterflies. PLoS Genet. 9(7):e1003620. the chemoreceptor superfamily. Genetics 177(3):1395–1416. Cheng T, et al. 2017. Genomic adaptation to polyphagy and insecticides in Mitter C, Farrell B, Wiegmann B. 1988. The phylogenetic study of adaptive a major East Asian noctuid pest. Nat Ecol Evol. 1(11):1747–1756. zones: has phytophagy promoted insect diversification? Am Nat. Clyne PJ, Warr CG, Carlson JR. 2000. Candidate taste receptors in 132(1):107–128. Drosophila. Science 287(5459):1830–1834. Nei M, Niimura Y, Nozawa M. 2008. The evolution of animal chemo- Couto A, Alenius M, Dickson BJ. 2005. Molecular, anatomical, and func- sensory receptor gene repertoires: roles of chance and necessity. Nat tional organization of the Drosophila olfactory system. Curr Biol. Rev Genet. 9(12):951–963. 15(17):1535–1547. Niimura Y, Matsui A, Touhara K. 2014. Extreme expansion of the olfactory Dahanukar A, Lei YT, Kwon JY, Carlson JR. 2007. Two Gr genes underlie receptor gene repertoire in African elephants and evolutionary dynam- sugar reception in Drosophila. Neuron 56(3):503–516. ics of orthologous gene groups in 13 placental mammals. Genome Delventhal R, Carlson JR. 2016. Bitter taste receptors confer diverse func- Res. 24(9):1485–1496. tions to neurons. Elife 5:1–23. Nishida R, Ohsugi T, Kokubo S, Fukami H. 1987. Oviposition stimulants of Duncan RP, et al. 2014. Dynamic recruitment of amino acid transporters to a Citrus-feeding swallowtail butterfly, Papilio xuthus L. Experientia the insect/symbiont interface. Mol Ecol. 23(6):1608–1623. 43(3):342–344. Ehrlich PR, Raven PH. 1964. Butterflies and plants: a study in coevolution. Normand P, et al. 2006. Genome characteristics of facultatively symbiotic Evolution 18(4):586–608. Frankia sp. strains reflect host range and host plant biogeography. Gardiner A, Barker D, Butlin RK, Jordan WC, Ritchie MG. 2008. Drosophila Genome Res. 17(1):7–15. chemoreceptor gene evolution: selection, specialization and genome Nylin S, Slove J, Janz N. 2014. Host plant utilization, host range oscillations size. Mol Ecol. 17(7):1648–1657. and diversification in nymphalid butterflies: a phylogenetic investiga- Gouin A, et al. 2017. Two genomes of highly polyphagous lepidopteran tion. Evolution 68(1):105–124. pests (Spodoptera frugiperda, Noctuidae) with different host-plant Ozaki K, et al. 2011. A gustatory receptor involved in host plant ranges. Sci Rep. 7(1):1–12. recognition for oviposition of a swallowtail butterfly. Nat Commun. Grabherr MG, et al. 2011. Trinity: reconstructing a full-length transcrip- 2:542. tome without a genome from RNA-Seq data. Nat Biotechnol. Pearce SL, et al. 2017. Genomic innovations, transcriptional plasticity and 29(7):644–652. gene loss underlying the evolution and divergence of two highly po- Guo H, et al. 2017. Expression map of a complete set of gustatory receptor lyphagous and invasive Helicoverpa pest species. BMC Biol. genes in chemosensory organs of Bombyx mori. Insect Biochem Mol 15(1):1–30. Biol. 82:74–82. Pereyra PC, Bowers MD. 1988. Iridoid glycosides as oviposition stimulants Huang X, Renwick JAA, Sachdev-Gupta K. 1993. Oviposition stimulants for the buckeye butterfly, Junonia coenia (Nymphalidae). J Chem Ecol. and deterrents regulating differential acceptance of Iberis amara by 14(3):917–928. Pieris rapae and P. napi oleracea. J Chem Ecol. 19(8):1645–1663. Renwick JAA, Chew FS. 1994. Oviposition behavior in Lepidoptera. Annu Janz N, Nylin S. 1997. The role of female search behaviour in determining Rev Entomol. 39(1):377–400. host plant range in plant feeding insects: a test of the information Robertson HM, Warr CG, Carlson JR. 2003. Molecular evolution of the processing hypothesis. Proc R Soc B Biol Sci. 264(1382):701–707. insect chemoreceptor gene superfamily in Drosophila melanogaster. Janz N, Nylin S. 2008. The oscillation hypothesis of host-plant range and Proc Natl Acad Sci U S A. 100(Suppl 2):14537–14542. speciation. Specialization, Speciation, and Radiation: The Evolutionary Robinson GS, Ackery PR, Kitching IJ, Beccaloni GW, Hern andez LM. 2010. Biology of Herbivorous Insects. University of California Press. p. 203– HOSTS - A database of the world’s Lepidopteran Hostplants. London: Natural History Museum. Janz N, Nylin S, Wahlberg N. 2006. Diversity begets diversity: host expan- Sanchez-Gracia A, Vieira FG, Rozas J. 2009. Molecular evolution of the sions and the diversification of plant-feeding insects. BMC Evol Biol. major chemosensory gene families in insects. Heredity 6:4. 103(3):208–216. Katoh K, Standley DM. 2013. MAFFT multiple sequence alignment soft- Sato K, Tanaka K, Touhara K. 2011. Sugar-regulated cation channel ware version 7: improvements in performance and usability. Mol Biol formed by an insect gustatory receptor. Proc Natl Acad Sci U S A. Evol. 30(4):772–780. 108(28):11680–11685. Kwon JY, Dahanukar A, Weiss LA, Carlson JR. 2007. The molecular basis Sawyer SA. 1989. Statistical tests for detecting gene conversion. Mol Biol of CO2 reception in Drosophila. Proc Natl Acad Sci U S A. Evol. 6(5):526–538. 104(9):3574–3578. Simao FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. Letunic I, Bork P. 2016. Interactive tree of life (iTOL) v3: an online tool for 2015. BUSCO: assessing genome assembly and annotation the display and annotation of phylogenetic and other trees. Nucleic completeness with single-copy orthologs. Bioinformatics Acids Res. 44(W1):W242–W245. 31(19):3210–3212. Li D, Zhang J. 2014. Diet shapes the evolution of the vertebrate bitter taste Stamatakis A. 2014. RAxML version 8: A tool for phylogenetic analysis and receptor gene repertoire. Mol Biol Evol. 31(2):303–309. post-analysis of large phylogenies. Bioinformatics 30(9):1312–1313. Genome Biol. Evol. 10(6):1351–1362 doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 1361 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Suzuki et al. GBE Stolzer M, et al. 2012. Inferring duplications, losses, transfers and incom- Weiss LA, Dahanukar A, Kwon JY, Banerjee D, Carlson JR. 2011. The plete lineage sorting with nonbinary species trees. Bioinformatics molecular and cellular basis of bitter taste in Drosophila. Neuron 28(18):i409–i415. 69(2):258–272. Stothard P. 2000. The sequence manipulation suite: javaScript programs Wiens JJ, Lapoint RT, Whiteman NK. 2015. Herbivory increases diversifi- for analyzing and formatting protein and DNA sequences. cation across insect clades. Nat Commun. 6(1):8370. Biotechniques 28:1102–1104. Xu W, Papanicolaou A, Zhang H-J, Anderson A. 2016. Expansion of a Tamate SC, Kawata M, Makino T. 2014. Contribution of nonohnologous bitter taste receptor family in a polyphagous insect herbivore. Sci duplicated genes to high habitat variability in mammals. Mol Biol Evol. Rep. 6(1):23666. 31(7):1779–1786. Yang Z. 2007. PAML 4: phylogenetic analysis by maximum likelihood. Mol Wahlberg N. 2006. That Awkward Age for Butterflies: insights from the Biol Evol. 24(8):1586–1591. Age of the Butterfly Subfamily Nymphalinae (Lepidoptera: nymphali- Zhan S, Merlin C, Boore JL, Reppert SM. 2011. The monarch butterfly dae). Syst Biol. 55(5):703–714. genome yields insights into long-distance migration. Cell Wahlberg N, et al. 2009. Nymphalid butterflies diversify following near 147(5):1171–1185. demise at the Cretaceous/Tertiary boundary. Proc R Soc B Biol Sci. Zhang X, et al. 2007. Characterizing the expression of the human olfactory 276(1677):4295–4302. receptor gene family using a novel DNA microarray. Genome Biol. Wahlberg N, Rubinoff D. 2011. Vagility across Vanessa (Lepidoptera: nym- 8(5):r86. phalidae): Mobility in butterfly species does not inhibit the formation Zhang X, et al. 2004. High-throughput microarray detection of olfactory and persistence of isolated sister taxa. Syst Entomol. 36(2):362–370. receptor gene expression in the mouse. Proc Natl Acad Sci U S A. Wahlberg N, Wheat CW, Pen ~a C. 2013. Timing and patterns in the tax- 101(39):14168–14173. onomic diversification of Lepidoptera (butterflies and moths). PLoS One 8(11):e80875–e80878. Associate editor: Belinda Chang 1362 Genome Biol. Evol. 10(6):1351–1362 doi:10.1093/gbe/evy093 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1351/4999382 by Ed 'DeepDyve' Gillespie user on 20 June 2018

Journal

Genome Biology and EvolutionOxford University Press

Published: May 18, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off