Skip to main content
Download PDF

Transcriptome and population structure of glassy-winged sharpshooters (Homalodisca vitripennis) with varying insecticide resistance in southern California

BMC Genomics volume 23, Article number: 721 (2022) Cite this article

Abstract

Background

Homalodisca vitripennis Germar, the glassy-winged sharpshooter, is an invasive insect in California and a critical threat to agriculture through its transmission of the plant pathogen, Xylella fastidiosa. Quarantine, broad-spectrum insecticides, and biological control have been used for population management of H. vitripennis since its invasion and subsequent proliferation throughout California. Recently wide-spread neonicotinoid resistance has been detected in populations of H. vitripennis in the southern portions of California’s Central Valley. In order to better understand potential mechanisms of H. vitripennis neonicotinoid resistance, we performed RNA sequencing on wild-caught insecticide-resistant and relatively susceptible sharpshooters to profile their transcriptome and population structure.

Results

We identified 81 differentially expressed genes with higher expression in resistant individuals. The significant largest differentially expressed candidate gene linked to resistance status was a cytochrome P450 gene with similarity to CYP6A9. Furthermore, we observed an over-enrichment of GO terms representing functions supportive of roles in resistance mechanisms (cytochrome P450s, M13 peptidases, and cuticle structural proteins). Finally, we saw no evidence of broad-scale population structure, perhaps due to H. vitripennis' relatively recent introduction to California or due to the relatively small geographic scale investigated here.

Conclusions

In this work, we characterized the transcriptome of insecticide-resistant and susceptible H. vitripennis and identified candidate genes that may be involved in resistance mechanisms for this species. Future work should seek to build on the transcriptome profiling performed here to confirm the role of the identified genes, particularly the cytochrome P450, in resistance in H. vitripennis. We hope this work helps aid future population management strategies for this and other species with growing insecticide resistance.

Peer Review reports

Background

The glassy-winged sharpshooter (GWSS), Homalodisca vitripennis Germar (Hemiptera: Cicadellidae), is a xylem-feeding leafhopper, which is invasive to California, and has proliferated since its introduction in the 1990s [1, 2]. GWSS has a broad host range, with over 340 reported plant species according to the California Department of Food and Agriculture (CDFA) (https://www.cdfa.ca.gov/pdcp/Documents/HostListCommon.pdf). In California’s agricultural systems, citrus is the major feeding, over-wintering, and reproductive host for GWSS, and the distribution of the insect is closely associated with the major citrus-growing regions in the state [2]. Although citrus can sustain high densities of GWSS, the major economic impact of the insect is as a vector of several strains of the xylem-limited, plant pathogenic bacterium, Xylella fastidiosa. This pathogen is the causal agent of several important diseases of crops in California, including Pierce's disease (PD) of grapes, oleander leaf scorch, and almond leaf scorch, and outside of California is also responsible for Citrus variegated chlorosis [3]. In particular, contiguous plantings of citrus and grapes have resulted in a significant increase in the incidence of PD, given the annual movement of GWSS from citrus to adjacent vineyards during the spring when grapes come out of dormancy [4].

The introduction of GWSS to California initially led to outbreaks of the destructive PD in the Temecula valley [1, 5, 6]. To combat the threat to the viticulture industry, area-wide treatments with the systemic neonicotinoid insecticide, imidacloprid, were undertaken and were successful at reducing the population sizes in both the southern Central Valley and Southern California in general [5, 6]. As GWSS began to spread and proliferate in the California Central Valley during the early 2000s, area-wide treatments with neonicotinoids (acetamiprid and imidacloprid) and pyrethroids were also introduced in Tulare and Kern counties [7]. The area-wide programs were highly successful until about 2012, when the levels of control appeared to be compromised [7, 8]. Recent work has found that applications of neonicotinoids have led to high levels of resistance in some GWSS populations, and it is believed that this was one of the major contributing factors to the population resurgence in the region [9, 10].

Insecticide resistance generally occurs through several co-existing processes spanning behavioral (e.g. avoidance), and physiological mechanisms (e.g. cuticle modifications, detoxification by host or symbionts, and target site alterations) [11,12,13]. Investigations of neonicotinoid resistance have consistently found enhanced detoxification by constitutively overexpressed cytochrome P450 monooxygenases in many insects spanning members of the orders Coleoptera (e.g., Tribolium castaneum, Leptinotarsa decemlineata), Diptera (e.g., Bradysia odoriphaga) and Hemiptera (e.g., Bemisia tabaci, Myzus persicae, Laodelphax striatellus, Rhopalosiphum padi, and Nilaparvata lugens) [14,15,16,17,18,19,20,21,22,23]. Transcriptome-based profiling approaches have enabled a broader understanding of resistance mechanisms, with many insects, in addition to having up-regulated cytochrome P450s, displaying differential expression of genes with functions related to cuticle structure and assembly, and detoxification through esterases, glutathione-S-transferases, or ABC transporters [24,25,26,27,28,29,30,31].

Despite the high levels of neonicotinoid and especially imidacloprid resistance observed in several populations of California GWSS, little is known about the underlying mechanisms involved in conferring resistance in these populations. In this study, we sought to profile the transcriptome of wild-caught resistant and susceptible GWSS obtained from different southern California populations. Specifically, we sought to (i) identify differentially expressed genes (DEGs) between resistant and susceptible populations, (ii) assess DEGs for functional enrichment that might relate to resistance mechanisms, and (iii) assess population structure in coding regions between resistant and susceptible populations to order to identify single nucleotide variants (SNVs) that may correlate with resistance status. Based on previous studies of neonicotinoid resistance, we hypothesized that we would see upregulation of one or more cytochrome P450s, as well as other genes related to cuticle modifications and detoxification. Understanding expression patterns of gene candidates linked to insecticide resistance in GWSS may help inform long-term solutions for population management of this species.

Results

Transcriptome identifies gene candidates linked to insecticide-resistance status

Using DESeq2, we identified 607 DEGs between resistant and susceptible GWSS populations (Table S1). Of these, 81 had higher expression in resistant populations and 526 had higher expression in susceptible populations. Of the DEGs, 57.3% had a functional annotation match to at least one database. Insect cuticle proteins were the dominant function of DEGs with higher expression in susceptible populations (IPR000618; n = 28). In contrast, M13 peptidases (IPR000718; n = 11) and cytochrome P450s (IPR001128; n = 8) were the dominant functions found with higher expression in resistant populations. Ordination of overall transcriptome profiles did not depict a strong pattern related to resistance and while many genes were identified as differentially expressed, only a single gene had an obvious consistent pattern linked to resistance status, J6590_005969 (Fig. 1).

Fig. 1
figure 1

Gene expression differs between GWSS populations and with resistance status. A Principal Component Analysis (PCA) of variance stabilized transcriptomic count data representing the full dataset. Each point represents an individual transcriptome sample. Samples are colored by collection location, while shapes are used to display resistance status (circle = resistant, triangle = susceptible). Resistant and susceptible populations are further highlighted by ellipses representing the 95% confidence interval around the centroid of each group. B Heatmap showing the variance stabilized counts of the 25 most significant DEG between resistant and susceptible populations. Each column represents a GWSS RNAseq sample with letters (A-D) representing populations and numbers (1–4) representing replicate (see Table 1 for additional sample information)

Full size image

Based on previous insecticide literature [14,15,16,17,18,19,20,21,22,23], we had hypothesized that cytochrome P450 genes might be partially or fully responsible for insecticide resistance of GWSS in the California Central Valley. Of the 81 DEGs identified as having higher expression in the resistant populations, eight were annotated as cytochrome P450s. Of these, a single gene, J6590_005969, was also the most significantly differentially expressed gene between resistant and susceptible populations (Fig. 2). Phylogenetic approaches place J6590_005969 in a clade with other GWSS cytochrome P450s belonging to the CYP6A9 family. Interestingly, this locus is flanked on both sides by Helitron, Mutator and CACTA transposons, similar to other observations of repetitive elements flanking cytochrome P450s involved in detoxification (e.g., [32]). None of the elements surrounding J6590_005969 were predicted by the Extensive de novo TE Annotator (EDTA) to be structurally intact (i.e. containing TIR regions, transposases or protein domains).

Fig. 2
figure 2

A single cytochrome P450 (J6590_005969) is most closely linked to resistance status. A Depiction of the architecture of the genomic region housing J6590_005969, including the 10 kb up and downstream of the candidate gene. Log-transformed coverage from resistant individuals is mapped to the region in gray. The gene (exons) is displayed in yellow. Repetitive elements as identified and annotated by EDTA are also shown in green (CACTA TIR transposon), blue (helitron), orange (Mutator TIR transposon) and pink (simple repetitive region). B Phylogeny of cytochrome P450s in the GWSS genome. J6590_005969 is highlighted in the tree in orange. A full phylogeny of cytochrome P450s in GWSS and their respective transcriptomic expression across samples can be found in Figure S1. C Mean of variance stabilized counts of reads mapping to J6590_005969 in resistant and susceptible populations. Bars represent standard error

Full size image

Functional enrichment of GO terms supports a landscape of resistance in GWSS

GO term enrichment analysis was performed on the DEGs to identify significantly over-enriched functional terms (Table 2, Bonferroni-adjusted p < 0.05). For DEGs with higher expression in insecticide-resistant GWSS, we found that GO terms for seven molecular functions (MF) and 15 biological processes (BP) were over-enriched, but no terms for cellular compartments (CC) were enriched. For DEGs with higher expression in susceptible GWSS, we found that GO terms for three MFs and three CCs were over-enriched, but found no enrichment in terms for BPs. Generally, over-enriched MF terms in resistant GWSS were represented by gene clusters with predicted functions such as cytochrome P450s (GO:0005506, GO:0016705, GO:0020037, GO:0046906), the M13 protease neprilysin (GO:0004222, GO:0008237) and vitellogenins (GO:0005319) (Fig. 3). In contrast, gene clusters represented by over-enriched MF terms in susceptible GWSS were dominated by genes with functions related to cuticle structural proteins (GO:0042302, GO:0005198) and peritrophins (GO:0008061). The functions of genes represented in over-enriched BP terms in insecticide-resistant GWSS include vitellogenins (GO:0006869,GO:0010876) and a cluster of three gene copies of the body color gene yellow (GO:0018958, GO:0046189, GO:0042440, GO:0046148, GO:0019953, GO:0032504, GO:0044703, GO:0048609, GO:0051704, GO:0000003, GO:0022414, GO:1901617, GO:1901615). The functions of over-enriched CC terms in susceptible insects include peritrophins (GO:0005576), and NADH-ubiquinone oxidoreductases and MICOS complex subunits (GO:0019866, GO:0005743).

Fig. 3
figure 3

Gene concept network showing links between genes with shared over-enriched molecular function GO terms. A Over-enriched GO terms from the molecular function (MF) category representing genes with higher expression in resistant GWSS populations. B Over-enriched MF GO terms representing genes with higher expression in susceptible GWSS populations. Each cluster is annotated with the dominant gene function as determined by the gene annotation and/or from best matches to UniProt (black boxes)

Full size image

No evidence of broad-scale population structure from coding region variants

Observed heterozygosity was significantly lower than expected (Bartlett’s test, p < 0.001) and FIS was 0.25 indicating high levels of non-random mating (e.g. inbreeding) in sampled GWSS. Overall Fst was 0.02, indicating low differentiation between populations (i.e., A-D; Table 1) and possibly high levels of migration or gene flow between populations, which is consistent with the relatively small geographic range studied here. Within-population FST were 0.028 and 0.033 for resistant and susceptible populations, respectively, and were 0.070, 0.081, 0.087, and 0.033 for population A (Tulare susceptible), B (Temecula susceptible), C (General Beale resistant) and D (Tulare resistant), respectively. Pairwise (between-population) FST was 0.011 between resistant and susceptible populations. Pairwise FST was similarly low for between A, B, C and D populations, ranging from 0.012 to 0.036 (Fig. 4).

Table 1 Sample information
Full size table
Fig. 4
figure 4

Low population differentiation across GWSS individuals from different collection locations and with different resistant statuses. A Pairwise FST values are shown here as heatmap with each value shown and FST values colored on a gradient with lowest values in blue and higher values in yellow. B PCA of population variation based on the SNV data. Each point represents an individual variant profile. Samples are colored by collection location while shapes are used to display resistance status (circle = resistant, triangle = susceptible). Resistant and susceptible populations are further highlighted by ellipses representing the 95% confidence interval around the centroid of each group

Full size image

Further, population structure results from fastSTRUCTURE (marginal likelihood maximized at K = 1, Figure S2) and Landscape and Ecological Association Studies (LEA) (lowest cross-entropy at K = 1; Tracy Widom test p > 0.05 for all PCA eigenvalues, Figure S3) were both supportive of weak differentiation between populations. Despite a lack of overall population structure, using PCA ordinations, we observed some minimal separation by the general location individuals were collected from (Fig. 4). Thus, when performing analysis of molecular variance (AMOVA) tests, we used resistance status alone (therefore assuming no population structure) and in a nested structure to account for separation by collection location or population. AMOVA tests on resistance status alone were significant (p < 0.01), but only explained a tiny portion, 1.78%, of observed variation. Most of the variation, 98.22%, was found within individuals, further supporting a panmictic (randomly mating) population. When using a nested structure, AMOVA tests of resistance status were not significant (p > 0.05), indicating no detectable signal of resistance when accounting for local population structure. In these cases, collection location or population were significant (p < 0.01), accounting for 4.05% of variation, indicating that while there is no detectable broad-scale population structure, there may still be some local adaptation occurring. We used OutFlank to identify Fst outliers as likely variants potentially under selection due to resistance status, but no significant outliers were detected (q > 0.1).

We used SNPEff annotations of the functional effect of variants to identify SNVs that might contribute to the increased expression of J6590_005969. Unfortunately, the low levels of expression of J6590_005969 in susceptible GWSS prevented our ability to call SNVs in this gene for susceptible samples. However, given that the reference genome represents a susceptible genotype [33], we still investigated the variants called in the resistant GWSS individuals. We found no SNVs that were homozygous or heterozygous for an alternate allele that was predicted to affect J6590_005969 function.

Discussion

Transcriptome results identified a cytochrome P450 candidate strongly linked to resistance

As has been seen in previous neonicotinoid resistance studies [14,15,16,17,18,19,20,21,22,23], we observed eight cytochrome P450s with higher expression in insecticide-resistant GWSS, with a robust pattern of expression in a single overexpressed cytochrome P450, J6590_005969, similar to CYP6A9, linked to resistance status. This gene candidate was flanked by repetitive elements, whose presence was previously hypothesized to be characteristic of xenobiotic-metabolizing P450s [32]. The mechanisms behind overexpression of many cytochrome P450s vary and part of our inability to identify SNVs responsible for J6590_005969’s higher expression may be because our analyses focused on coding regions. In other insects, overexpression of cytochrome P450s has been found to be regulated by trans and/or cis regulatory elements [34,35,36]. In N. lugens, for example, SNVs in the promoter region of CYP6AY1 were observed to enhance promoter activity and hypothesized to be acting as cis-acting factors that enhance expression in resistant individuals [37]. In Drosophila melanogaster, the characterization of several CYP6A genes with higher expression in resistant individuals suggested that a defective repressor might be involved in regulation of DDT resistance [38]. Future work should focus on confirming J6590_005969’s role in resistance and characterizing the upstream promoter region of this gene to investigate the molecular mechanism driving its overexpression in resistant individuals. New advances demonstrating the efficiency of using CRISPR mutagenesis in GWSS may dramatically accelerate opportunities to examine the functional capabilities of J6590_005969 in vivo [39, 40].

Functional enrichment suggests a complex multi-gene response to insecticides

The transcriptome of GWSS has been previously described as being dominated by genes annotated with GO terms related to molecular binding, catalytic (e.g., hydrolase and oxidoreductases), and transporter activity [41]. These categories broadly match with the GO terms observed to be over-enriched in insecticide-resistant populations here (Table 2). Specifically, in addition to cytochrome P450s, we characterized an overall transcriptome pattern in resistant individuals supportive of additional genes being directly involved in or indirectly affected by insecticide resistance mechanisms including genes related to detoxification, immune response and cuticle modifications. Genes indirectly affected by resistance mechanisms may be the result of trade-offs and represent a potential fitness cost of insecticide resistance [42,43,44,45,46]. Although we were able to observe an over-enrichment of functions that may be related to resistance, 42.7% of DEGs had no functional annotation. Additional novel genes or functions may be involved in resistance that are not discussed here. Studies focused on the molecular characterization of unannotated genes in this and other insect species are needed to help close these annotation gaps.

Table 2 Over-enriched GO terms for genes that are differentially expressed between resistant and susceptible GWSS
Full size table

We observed upregulation of neprilysin (M13 peptidases) and vitellogenin-like genes in resistant GWSS. Neprilysin and neprilysin-like proteins are zinc metalloendopeptidases and are type II integral membrane proteins that turn off signaling events at the cell surface [47]. Previous studies in insects have seen altered expression of M13 peptidases during metamorphosis and immune responses [48,49,50,51]. Vitellogenins are important for insect reproduction. They play a role in immune responses and protect against oxidative stress induced by insecticides in bees [52,53,54]. Additionally, the expression of vitellogenins in the white-backed planthopper has been shown to be altered by insecticide application [55]. Enhanced expression of neprilysin and vitellogenin-like genes in GWSS therefore may be related to increased immune response, or alternatively to a reproductive fitness cost, in resistant individuals.

Compared to susceptible GWSS, resistant individuals had a lower expression of cuticle structural and peritrophin-like genes. The cuticle is the first barrier of protection for insects against insecticides. Expression changes in genes predicted to be involved in cuticle structure have been observed in a variety of insect species [16, 27, 28, 56]. While some studies of hemipteran pests have observed an upregulation of cuticle genes, others in M. persicae and Aphis gossypii, have reported their downregulation such as was seen here for GWSS [24, 27]. While future work is needed to confirm this, one possibility is that these expression changes, whether up or down, are leading to cuticle modifications that could contribute to resistance (e.g. thicker cuticles). Alternatively, the down regulation of these genes may indicate a fitness cost, with resistant individuals prioritizing detoxification mechanisms over barrier protection. Peritrophins are proteins with chitin-binding domains that are an integral part of the peritrophic membrane, which lines the insect gut. The peritrophic membrane is thought to aid in digestion and protection from toxins [57]. Previous work has shown altered expression of peritrophin-like proteins linked to cycloxaprid application [58] and silencing peritrophins can lead to higher imidacloprid susceptibility in termites [59]. While it's unclear here why we observe reduced expression of pertriphonins in resistant individuals, it is possible that these genes have roles in resistance mechanisms through altered gut structure or detoxification. It is also possible that these represent another trade-off related to gut structure and digestion in resistant individuals.

A lack of broad-scale population structure is consistent with previous work

GWSS has been proliferating in California since its initial invasion in the 1990s, where it was introduced via humans likely through nursery shipments [1, 2]. Given this recent introduction, and the relatively small geographic scale examined here, it is perhaps not surprising to see limited broad-scale population structure, especially when considering coding regions. Previous work by Stenger et al. used a GWSS reovirus, due to its faster rate of evolution, to date the introduction of GWSS to California [2]. Their results pointed at an introduction followed by a bottleneck expansion, which is in agreement with the levels of inbreeding observed here. Our results are also consistent with a study by Smith (2015), who used the mitochondrial cytochrome oxidase I (COI) gene to investigate population structure of GWSS across the United States [60]. Their results showed two distinct groups of haplotypes of GWSS in the United States, a group of populations from east of the Mississippi River including Louisiana, Mississippi, Alabama and Florida, and a group composed of populations west of the Mississippi River from Texas and California. However, they found that neither group of haplotypes had sufficient genetic structuring to further differentiate populations within the two groups. Furthermore, the data supported the hypothesis that GWSS populations in California most likely originated from a source in Texas, arising from random distribution by humans, and not from GWSS from east of the Mississippi River. The subsequent distribution of GWSS throughout California could promote gene-flow within populations that would contribute to the low differentiation and lack of population structure seen here.

Conclusion

We provide the first, to our knowledge, characterization of the transcriptome of neonicotinoid resistant and susceptible H. vitripennis. We identified a suite of candidate genes linked to resistance status including a highly expressed cytochrome P450 (J6590_005969), and observed additional expression patterns supportive of multi-gene roles in resistance mechanisms or fitness trade-offs. While we found no evidence of broad-scale population structure, this may be due to the recent introduction of GWSS to California, relatively small geographic range investigated here, or continued gene flow due to accidental distribution by humans. Follow up work is needed to investigate the specific functional roles and molecular mechanisms responsible for the DEGs identified here, particularly the upregulated cytochrome P450 (J6590_005969), and to confirm whether these DEGs can affect GWSS resistance in vivo. Additionally, future studies should also consider whether obligate and facultative microbial symbionts of GWSS are involved in conferring insecticide resistance (e.g. [61]). We believe that this work serves as a foundation for future studies of insecticide resistance in GWSS and other Hemipteran insects.

Methods

Sample collection and sequencing

Sharpshooters were previously collected in August 2019 from California citrus groves in Porterville (Tulare-Organic), Temecula (Temecula-Organic), Bakersfield (General Beale-Organic) and Terra Bella (Tulare-Conventional) (Fig. 5) as part of a multi-year monitoring program and were confirmed to have varying levels of neonicotinoid resistance [9]. Baseline susceptible levels for imidacloprid were determined in 2003 using populations that had never been exposed to imidacloprid. Using a range of imidacloprid doses, insects were treated by topically applying insecticide to their abdomens. Mortality was assessed at 48 h, and LD50s (the lethal dose that results in 50% mortality) derived from the dose–response curves using probit analysis. Resistance was defined when the LD50 of a field population was statistically significant from the LD50 of a susceptible population. However, it was not possible to determine an LD50 for the General Beale-Organic or the Tulare-Conventional populations due to the high levels of resistance. Therefore, a discriminating dose of 500 ng/insect was chosen to distinguish susceptible and resistant insects [9]. In 2020, insects from the General Beale-Organic population were tested at the discriminating dose, resulting in only 16% mortality (FJB unpublished). The latter result continued a consistent pattern of high resistance in GWSS collected from that region of California since the first tests were conducted in 2017. In contrast, tests with insects from the Temecula-Organic population confirmed full susceptibility to imidacloprid in 2019 [9], and in 2020 when the most recent data were determined (FJB unpublished).

Fig. 5
figure 5

Map depicting collection locations. Map diagram showing California with points depicting different populations. Points are labeled by population (A-D; Table 1), colored by collection location, and shapes are used to display resistance status (circle = resistant, triangle = susceptible)

Full size image

Four sharpshooters from each of these four locations were chosen from the 2019 collections to produce a total of 16 sharpshooter transcriptomes (Table 1) [9, 33]. For the resistant locations (Tulare-Conventional and General Beale-Organic), GWSS were treated with imidacloprid to confirm resistance (as described above) and healthy survivors (resistant GWSS) were sequenced. While for the susceptible locations (Tulare-Organic and Temecula-Organic), a proportion of collected GWSS were tested with imidacloprid to confirm susceptibility levels and untreated GWSS from the same collections were sequenced. For each sharpshooter, RNA was extracted from adult prothoracic leg tissue using the Monarch Total RNA Mini Kit (New England Biolabs, Ipswich, MA). Paired-end RNA-Seq libraries were constructed with NEBNext Ultra II Directional RNA prep (New England Biolabs, Ipswich, MA) and sequenced on NovaSeq 6000 to produce an average of 87 M paired reads per library (minimum library 51 M, max library 124 M reads). The RNASeq data is available on NCBI Genbank under BioProject PRJNA717315. Computational scripts associated with analysis in this manuscript are available on GitHub and archived in Zenodo [62].

Transcriptome expression and functional enrichment

The H. vitripennis UCR_GWSS_1.0 assembly and the annotation described previously and available on NCBI at JAGXCG010000000 were used here as a reference [33]. EDTA v.1.9.4 was used to annotate repetitive elements for this assembly using the following parameters: –anno 1 –evaluate 1 –sensitive 1 –step all [63]. Transcriptome reads were aligned against the H. vitripennis UCR_GWSS_1.0 reference genome using STAR v.2.7.9a to generate bam files and then these files were summarized at the gene level using featureCounts v.1.6.2 [64, 65]. These count tables were uploaded into R v.4.1.2 and analyzed using the DESeq2 package v.1.34.0 to examine the log2 fold change (i.e. differential expression) of genes between insecticide-resistant and susceptible populations [66, 67]. We subsequently focused on genes with Bonferroni corrected p-values < 0.05 that had a log2 fold change >|2|. Enrichment analysis of Gene Ontology (GO) terms was performed on differentially expressed genes to identify GO terms that were significantly over-enriched (p < 0.05). This analysis was performed for each of the three GO classes (i.e., biological processes [BP], molecular functions [MF], and cellular components [CC]). Differential expression and GO enrichment results were visualized in R using clusterProfiler v.4.2.2 [68]. When graphing expression trends, count data was normalized using the variance stabilizing transformation. The R package gggenomes v.0.9.5.9000 was used to visualize genes of interest [69].

Predicted cytochrome P450 genes were aligned using MUSCLE v.3.8.425 [70]. The resulting alignment was then trimmed using the -automated1 option in trimAl v.1.4.1 [71]. A maximum likelihood tree was then built from the trimmed alignment using IQTREE2 v.2.1.3 with 1000 bootstraps [72]. This phylogeny was imported into R for visualization with ggtree v.3.2.1 [73].

Variant calling

Transcriptome bam files were processed using ‘AddOrReplaceReadGroups’, ‘MarkDuplicates’ and ‘​​SplitNCigarReads’ to assign reads to new sample groups, flag duplicate reads and split reads containing N’s using Picard tools (http://broadinstitute.github.io/picard). The variants (SNVs and indels) were genotyped relative to the H. vitripennis UCR_GWSS_1.0 reference genome using the HaplotypeCaller step in GATK v4.0 [74]. Predicted variants were filtered using GATK’s SelectVariants call with the following parameters: for SNVs, -window-size = 10, -QualByDept < 2.0, -MapQual < 40.0, -QScore < 100, -MapQualityRankSum <  − 12.5, -StrandOddsRatio > 4.0, -FisherStrandBias > 60.0, -ReadPosRankSum <  − 8.0; for indels, -window-size = 10, -QualByDepth < 2.0, -MapQualityRankSum <  − 12.5, -StrandOddsRatio > 10.0, -FisherStrandBias > 200.0, -ReadPosRank <  − 20.0, -InbreedingCoeff <  − 0.8. Variants were subsequently annotated with snpEff [75]. We then used VCFtools v.0.1.16 to investigate missingness across samples using –missing-indv, and further filtered the variant table using –missingness 0.75 and –mac 3 [76]. Resulting variant tables contained 300,365 polymorphic SNVs. Using VCFtools, we converted final VCF files into plink format for some analyses [77].

Population analysis

Variant tables were imported into R in VCF and plink formats for analysis. We used hierfstat v.0.5–10 to calculate basic population statistics, such as observed and expected heterozygosity, using the basic.stats function [78]. We also used the hierfstat package to calculate the fixation index statistics including FST and FIS both within and between populations. We then used fastSTRUCTURE to assess broad-scale population structure using K values from 1 to 30 [79]. The ‘chooseK.py’ function was used to assess which K provided the best marginal likelihood value. We also investigated broad-scale structure using a complementary method, LEA, with K values from 1 to 30 across ten iterations [80]. We then performed AMOVA tests on population stratification using poppr v.2.9.3 and ade4 v.1.7–18 [81, 82]. We ran AMOVAs on resistance status (resistant, susceptible) alone and then in nested hierarchies with collection county (Tulare, Temecula, General Beale), population (A-D), or both. Finally, we used OutFLANK v.0.2 to search for and identify FST outliers that might be linked to resistance status [83].

Availability of data and materials

The transcriptome data is publicly available on DDBJ/ENA/GenBank under BioProject PRJNA717315. Computational scripts associated with the analysis of this data are available on GitHub and archived in Zenodo [62].

References

  1. Blua MJ, Phillips PA, Redak RA. A new sharpshooter threatens both crops and ornamentals. Calif Agric. 1999;53:22–5.

    Google Scholar 

  2. Stenger DC, Sisterson MS, French R. Population genetics of Homalodisca vitripennis reovirus validates timing and limited introduction to California of its invasive insect host, the glassy-winged sharpshooter. Virology. 2010;407:53–9.

    Article  PubMed  CAS  Google Scholar 

  3. Hopkins DL. Xylella Fastidiosa: Xylem-Limited Bacterial Pathogen of Plants. Annu Rev Phytopathol. 1989;27:271–90.

    Article  Google Scholar 

  4. Redak RA, Purcell AH, Lopes JRS, Blua MJ, Mizell RF 3rd, Andersen PC. The biology of xylem fluid-feeding insect vectors of Xylella fastidiosa and their relation to disease epidemiology. Annu Rev Entomol. 2004;49:243–70.

    Article  PubMed  CAS  Google Scholar 

  5. Hix R, Toscano N, Redak R, Blua M. Area-wide management of the glassy-winged sharpshooter in the Temecula Valley. In: Proceedings of the Pierce’s Disease Research Symposium. Sacramento: California Dept of Food and Agriculture; 2002. 157–8.

  6. Toscano NC, Hix R, Gispert C. Riverside county glassy-winged sharpshooter area-wide management program in the Coachella and Temecula valleys. In: Proceedings of the Pierce’s Disease Research Symposium. Sacramento: California Department of Food and Agriculture; 2005. 384–7.

  7. Stone-Smith B, Stewart-Leslie J, Kunkel G, Hardy S, Appleby B, Borges D. Area-wide pest management of the glassy-winged sharpshooter in Tulare County. In: Proceedings of the Pierce’s Disease Research Symposium. Sacramento: California Department of Food and Agriculture; 2005:376–9.

  8. Haviland D, Stone-Smith B. Monitoring glassy-winged sharpshooter and Pierce’s disease in Kern County, California. In: Proceedings of the Pierce’s Disease Research Symposium. Sacramento: California Dept of Food and Agriculture; 2016. 75–80.

  9. Byrne FJ, Redak RA. Insecticide resistance in California populations of the glassy-winged sharpshooter Homalodisca vitripennis. Pest Manag Sci. 2021;77:2315–23.

    Article  PubMed  CAS  Google Scholar 

  10. Andreason SA, Prabhaker N, Castle SJ, Ganjisaffar F, Haviland DR, Stone-Smith B, et al. Reduced susceptibility of Homalodisca vitripennis (Hemiptera: Cicadellidae) to commonly applied insecticides. J Econ Entomol. 2018;111:2340–8.

    Article  PubMed  CAS  Google Scholar 

  11. Bass C, Denholm I, Williamson MS, Nauen R. The global status of insect resistance to neonicotinoid insecticides. Pestic Biochem Physiol. 2015;121:78–87.

    Article  PubMed  CAS  Google Scholar 

  12. Dawkar VV, Chikate YR, Lomate PR, Dholakia BB, Gupta VS, Giri AP. Molecular insights into resistance mechanisms of lepidopteran insect pests against toxicants. J Proteome Res. 2013;12:4727–37.

    Article  PubMed  CAS  Google Scholar 

  13. Dang K, Doggett SL, Veera Singham G, Lee CY. Insecticide resistance and resistance mechanisms in bed bugs, Cimex spp. (Hemiptera: Cimicidae). Parasit Vectors. 2017;10:318.

    Article  PubMed  PubMed Central  Google Scholar 

  14. Karunker I, Benting J, Lueke B, Ponge T, Nauen R, Roditakis E, et al. Over-expression of cytochrome P450 CYP6CM1 is associated with high resistance to imidacloprid in the B and Q biotypes of Bemisia tabaci (Hemiptera: Aleyrodidae). Insect Biochem Mol Biol. 2008;38:634–44.

    Article  PubMed  CAS  Google Scholar 

  15. Bass C, Zimmer CT, Riveron JM, Wilding CS, Wondji CS, Kaussmann M, et al. Gene amplification and microsatellite polymorphism underlie a recent insect host shift. Proc Natl Acad Sci. 2013;110:19460–5.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  16. Puinean AM, Foster SP, Oliphant L, Denholm I, Field LM, Millar NS, et al. Amplification of a cytochrome P450 gene is associated with resistance to neonicotinoid insecticides in the aphid Myzus persicae. PLoS Genet. 2010;6: e1000999.

    Article  PubMed  PubMed Central  Google Scholar 

  17. Bass C, Carvalho RA, Oliphant L, Puinean AM, Field LM, Nauen R, et al. Overexpression of a cytochrome P450 monooxygenase, CYP6ER1, is associated with resistance to imidacloprid in the brown planthopper Nilaparvata lugens. Insect Mol Biol. 2011;20:763–73.

    Article  PubMed  CAS  Google Scholar 

  18. Ding Z, Wen Y, Yang B, Zhang Y, Liu S, Liu Z, et al. Biochemical mechanisms of imidacloprid resistance in Nilaparvata lugens: Over-expression of cytochrome P450 CYP6AY1. Insect Biochem Mol Biol. 2013;43:1021–7.

    Article  PubMed  CAS  Google Scholar 

  19. Liang X, Xiao D, He Y, Yao J, Zhu G, Zhu KY. Insecticide-mediated up-regulation of cytochrome P450 genes in the red flour beetle (Tribolium castaneum). Int J Mol Sci. 2015;16:2078–98.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  20. Elzaki MEA, Miah MA, Wu M, Zhang H, Pu J, Jiang L, et al. Imidacloprid is degraded by CYP353D1v2, a cytochrome P450 overexpressed in a resistant strain of Laodelphax striatellus. Pest Manag Sci. 2017;73:1358–63.

    Article  PubMed  CAS  Google Scholar 

  21. Zhu F, Moural TW, Nelson DR, Palli SR. A specialist herbivore pest adaptation to xenobiotics through up-regulation of multiple Cytochrome P450s. Sci Rep. 2016;6:20421.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  22. Chen C, Wang C, Liu Y, Shi X, Gao X. Transcriptome analysis and identification of P450 genes relevant to imidacloprid detoxification in Bradysia odoriphaga. Sci Rep. 2018;8:2564.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Wang K, Zhang M, Huang Y, Yang Z, Su S, Chen M. Characterisation of imidacloprid resistance in the bird cherry-oat aphid, Rhopalosiphum padi, a serious pest on wheat crops. Pest Manag Sci. 2018;74:1457–65.

    Article  PubMed  CAS  Google Scholar 

  24. Meng J, Chen X, Zhang C. Transcriptome-based identification and characterization of genes responding to imidacloprid in Myzus persicae. Sci Rep. 2019;9:13285.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Zhou C, Yang H, Wang Z, Long GY, Jin DC. Comparative transcriptome analysis of Sogatella furcifera (Horváth) exposed to different insecticides. Sci Rep. 2018;8:8773.

    Article  PubMed  PubMed Central  Google Scholar 

  26. Yang N, Xie W, Jones CM, Bass C, Jiao X, Yang X, et al. Transcriptome profiling of the whitefly Bemisia tabaci reveals stage-specific gene expression signatures for thiamethoxam resistance. Insect Mol Biol. 2013;22:485–96.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  27. Pan Y, Peng T, Gao X, Zhang L, Yang C, Xi J, et al. Transcriptomic comparison of thiamethoxam-resistance adaptation in resistant and susceptible strains of Aphis gossypii Glover. Comp Biochem Physiol Part D Genomics Proteomics. 2015;13:10–5.

    Article  PubMed  CAS  Google Scholar 

  28. Gao Y, Kim K, Kwon DH, Jeong IH, Clark JM, Lee SH. Transcriptome-based identification and characterization of genes commonly responding to five different insecticides in the diamondback moth Plutella xylostella. Pestic Biochem Physiol. 2018;144:1–9.

    Article  PubMed  CAS  Google Scholar 

  29. Antony B, Johny J, Abdelazim MM, Jakše J, Al-Saleh MA, Pain A. Global transcriptome profiling and functional analysis reveal that tissue-specific constitutive overexpression of cytochrome P450s confers tolerance to imidacloprid in palm weevils in date palm fields. BMC Genomics. 2019;20:440.

    Article  PubMed  PubMed Central  Google Scholar 

  30. Zhou CS, Lv HH, Guo XH, Cao Q, Zhang RX, Ma DY. Transcriptional analysis of Bemisia tabaci MEAM1 cryptic species under the selection pressure of neonicotinoids imidacloprid, acetamiprid and thiamethoxam. BMC Genomics. 2022;23:15.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Zhang B-Z, Ma K-S, Liu J-J, Lu L-Y, Chen X-L, Zhang S-P, et al. Differential expression of genes in greenbug (Schizaphis graminum Rondani) treated by imidacloprid and RNA interference. Pest Manag Sci. 2019;75:1726–33.

    Article  PubMed  CAS  Google Scholar 

  32. Chen S, Li X. Transposable elements are enriched within or in close proximity to xenobiotic-metabolizing cytochrome P450 genes. BMC Evol Biol. 2007;7:46.

    Article  PubMed  PubMed Central  Google Scholar 

  33. Ettinger CL, Byrne FJ, Collin MA, Carter-House D, Walling LL, Atkinson PW, et al. Improved draft reference genome for the Glassy-winged Sharpshooter (Homalodisca vitripennis), a vector for Pierce’s disease. G3. 2021;11:jkab255.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  34. Maitra S, Dombrowski SM, Basu M, Raustol O, Waters LC, Ganguly R. Factors on the third chromosome affect the level of cyp6a2 and cyp6a8 expression in Drosophila melanogaster. Gene. 2000;248:147–56.

    Article  PubMed  CAS  Google Scholar 

  35. Liu N, Scott JG. Genetic analysis of factors controlling high-level expression of cytochrome P450, CYP6D1, cytochrome b5, P450 reductase, and monooxygenase activities in LPR house flies Musca domestica. Biochemical Genetics. 1996;34:133–48.

    Article  PubMed  CAS  Google Scholar 

  36. Cariño FA, Koener JF, Plapp FW, Feyereisen R. Constitutive overexpression of the cytochrome P450 gene CYP6A1 in a house fly strain with metabolic resistance to insecticides. Insect Biochem Mol Biol. 1994;24:411–8.

    Article  PubMed  Google Scholar 

  37. Pang R, Li Y, Dong Y, Liang Z, Zhang Y, Zhang W. Identification of promoter polymorphisms in the cytochrome P450 CYP6AY1 linked with insecticide resistance in the brown planthopper Nilaparvata lugens. Insect Mol Biol. 2014;23:768–78.

    Article  PubMed  CAS  Google Scholar 

  38. Maitra S, Dombrowski SM, Waters LC, Ganguly R. Three second chromosome-linked clustered Cyp6 genes show differential constitutive and barbital-induced expression in DDT-resistant and susceptible strains of Drosophila melanogaster. Gene. 1996;180:165–71.

    Article  PubMed  CAS  Google Scholar 

  39. de Souza Pacheco I, Doss A-LA, Vindiola BG, Brown DJ, Ettinger CL, Stajich JE, et al. Efficient CRISPR/Cas9-mediated genome modification of the glassy-winged sharpshooter Homalodisca vitripennis (Germar). Sci Rep. 2022;12:1–16.

    Google Scholar 

  40. Pacheco ID, Walling LL, Atkinson PW. Gene editing and genetic control of hemipteran pests: progress challenges and perspectives. Front Bioeng Biotechnol. 2022;10: 900785.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Nandety RS, Kamita SG, Hammock BD, Falk BW. Sequencing and de novo assembly of the transcriptome of the glassy-winged sharpshooter (Homalodisca vitripennis). PLoS ONE. 2013;8: e81681.

    Article  PubMed  PubMed Central  Google Scholar 

  42. Silva AX, Bacigalupe LD, Luna-Rudloff M, Figueroa CC. Insecticide resistance mechanisms in the green peach aphid Myzus persicae (Hemiptera: Aphididae) II: Costs and benefits. PLoS ONE. 2012;7: e36810.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  43. Homem RA, Buttery B, Richardson E, Tan Y, Field LM, Williamson MS, et al. Evolutionary trade-offs of insecticide resistance - The fitness costs associated with target-site mutations in the nAChR of Drosophila melanogaster. Mol Ecol. 2020;29:2661–75.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  44. Diniz DFA, de Melo-Santos MAV, de Mendonça Santos EM, Beserra EB, Helvecio E, de Carvalho-Leandro D, et al. Fitness cost in field and laboratory Aedes aegypti populations associated with resistance to the insecticide temephos. Parasites & Vectors. 2015;8:1–15.

    Article  Google Scholar 

  45. Gordon JR, Potter MF, Haynes KF. Insecticide resistance in the bed bug comes with a cost. Sci Rep. 2015;5:10807.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  46. Osoro JK, Machani MG, Ochomo E, Wanjala C, Omukunda E, Munga S, et al. Insecticide resistance exerts significant fitness costs in immature stages of Anopheles gambiae in western Kenya. Malar J. 2021;20:259.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  47. Turner AJ, Isaac RE, Coates D. The neprilysin (NEP) family of zinc metalloendopeptidases: genomics and function. BioEssays. 2001;23:261–9.

    Article  PubMed  CAS  Google Scholar 

  48. Zhao X, Mita K, Shimada T, Okano K, Quan GX, Kanke E, et al. Isolation and expression of an ecdysteroid-inducible neutral endopeptidase 24.11-like gene in wing discs of Bombyx mori. Insect Biochem Mol Biol. 2001;31:1213–9.

    Article  PubMed  CAS  Google Scholar 

  49. Aguilar R, Jedlicka AE, Mintz M, Mahairaki V, Scott AL, Dimopoulos G. Global gene expression analysis of Anopheles gambiae responses to microbial challenge. Insect Biochem Mol Biol. 2005;35:709–19.

    Article  PubMed  CAS  Google Scholar 

  50. Wilson CL, Shirras AD, Isaac RE. Extracellular peptidases of imaginal discs of Drosophila melanogaster. Peptides. 2002;23:2007–14.

    Article  PubMed  CAS  Google Scholar 

  51. Zhu Y, Johnson TJ, Myers AA, Kanost MR. Identification by subtractive suppression hybridization of bacteria-induced genes expressed in Manduca sexta fat body. Insect Biochem Mol Biol. 2003;33:541–59.

    Article  PubMed  CAS  Google Scholar 

  52. Banerjee BD, Seth V, Ahmed RS. Pesticide-induced oxidative stress: perspective and trends. Rev Environ Health. 2001;16:1–40.

    Article  PubMed  CAS  Google Scholar 

  53. Amdam GV, Simões ZLP, Hagen A, Norberg K, Schrøder K, Mikkelsen Ø, et al. Hormonal control of the yolk precursor vitellogenin regulates immune function and longevity in honeybees. Exp Gerontol. 2004;39:767–73.

    Article  PubMed  CAS  Google Scholar 

  54. Seehuus S-C, Norberg K, Gimsa U, Krekling T, Amdam GV. Reproductive protein protects functionally sterile honey bee workers from oxidative stress. Proc Natl Acad Sci U S A. 2006;103:962–7.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  55. Zhou C, Yang X-B, Yang H, Long G-Y, Jin D-C. Effects of Sublethal Concentrations of Insecticides on the Fecundity of Sogatella furcifera (Hemiptera: Delphacidae) via the Regulation of Vitellogenin and Its Receptor. J Insect Sci. 2020;20:14.

    Article  PubMed  PubMed Central  Google Scholar 

  56. Cao C, Sun L, Wen R, Shang Q, Ma L, Wang Z. Characterization of the transcriptome of the Asian gypsy moth Lymantria dispar identifies numerous transcripts associated with insecticide resistance. Pestic Biochem Physiol. 2015;119:54–61.

    Article  PubMed  CAS  Google Scholar 

  57. Rodríguez-de la Noval C, Rodríguez-Cabrera L, Izquierdo L, Espinosa LA, Hernandez D, Ponce M, et al. Functional expression of a peritrophin Alike SfPER protein is required for larval development in Spodoptera frugiperda (Lepidoptera: Noctuidae). Sci Rep. 2019;9:2630.

    Article  PubMed  PubMed Central  Google Scholar 

  58. Qi H, Liu T, Lu Q, Yang Q. Molecular insights into the insensitivity of Lepidopteran pests to Cycloxaprid. J Agric Food Chem. 2020;68:982–8.

    Article  PubMed  CAS  Google Scholar 

  59. Sandoval-Mojica AF, Scharf ME. Silencing gut genes associated with the peritrophic matrix of Reticulitermes flavipes (Blattodea: Rhinotermitidae) increases susceptibility to termiticides. Insect Mol Biol. 2016;25:734–44.

    Article  PubMed  CAS  Google Scholar 

  60. Smith PT. Mitochondrial DNA variation among populations of the glassy-winged sharpshooter. Homalodisca coagulata J Insect Sci. 2005;5:41.

    PubMed  CAS  Google Scholar 

  61. Ettinger CL, Byrne FJ, Redak RA, Stajich JE. Metagenome-assembled genomes of bacterial symbionts associated with insecticide-resistant and -susceptible individuals of the glassy-winged sharpshooter (Homalodisca vitripennis). Microbiol Resour Announc. 2022;11: e0050622.

    Article  PubMed  Google Scholar 

  62. Ettinger CL, Stajich JE. stajichlab/GWSS_Resistance. 2022. https://doi.org/10.5281/zenodo.6959134.

    Book  Google Scholar 

  63. Ou S, Su W, Liao Y, Chougule K, Agda JRA, Hellinga AJ, et al. Benchmarking transposable element annotation methods for creation of a streamlined, comprehensive pipeline. Genome Biol. 2019;20:275.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  64. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21.

    Article  PubMed  CAS  Google Scholar 

  65. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30:923–30.

    PubMed  CAS  Google Scholar 

  66. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.

    Article  PubMed  PubMed Central  Google Scholar 

  67. R Core Team. R: a language and environment for statistical computing. 2021.

    Google Scholar 

  68. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16:284–7.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  69. Hackl T, Hackl T. gggenomes: a grammar of graphics for comparative genomics. 2021.

    Google Scholar 

  70. Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004;32:1792–7.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  71. Capella-Gutiérrez S, Silla-Martínez JM, Gabaldón T. trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics. 2009;25:1972–3.

    Article  PubMed  PubMed Central  Google Scholar 

  72. Minh BQ, Schmidt HA, Chernomor O, Schrempf D, Woodhams MD, von Haeseler A, et al. IQ-TREE 2: new models and efficient methods for phylogenetic inference in the genomic era. Mol Biol Evol. 2020;37:1530–4.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  73. Yu G, Smith DK, Zhu H, Guan Y, Lam TT. ggtree : an r package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods Ecol Evol. 2017;8:28–36.

    Article  Google Scholar 

  74. Poplin R, Ruano-Rubio V, DePristo MA, Fennell TJ, Carneiro MO, Van der Auwera GA, et al. Scaling accurate genetic variant discovery to tens of thousands of samples. bioRxiv. 2018. https://doi.org/10.1101/201178.

  75. Cingolani P, Platts A, Wang LL, Coon M, Nguyen T, Wang L, et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly. 2012;6:80–92.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  76. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, et al. The variant call format and VCFtools. Bioinformatics. 2011;27:2156–8.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  77. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81:559–75.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  78. Goudet J. hierfstat, a package for r to compute and test hierarchical F-statistics. Mol Ecol Notes. 2005;5:184–6.

    Article  Google Scholar 

  79. Raj A, Stephens M, Pritchard JK. fastSTRUCTURE: variational inference of population structure in large SNP data sets. Genetics. 2014;197:573–89.

    Article  PubMed  PubMed Central  Google Scholar 

  80. Frichot E, François O. LEA: An R package for landscape and ecological association studies. Methods Ecol Evol. 2015;6:925–9.

    Article  Google Scholar 

  81. Kamvar ZN, Tabima JF, Grünwald NJ. Poppr: an R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction. PeerJ. 2014;2: e281.

    Article  PubMed  PubMed Central  Google Scholar 

  82. Dray S, Dufour A-B. The ade4 Package: Implementing the Duality Diagram for Ecologists. Journal of Statistical Software. 2007;22:1–22.

    Article  Google Scholar 

  83. Whitlock MC, Lotterhos KE. Reliable detection of loci responsible for local adaptation: inference of a null model through trimming the distribution of FST. Am Nat. 2015;186:S24–36.

    Article  PubMed  Google Scholar 

Download references

Acknowledgements

We would like to thank Beatriz G. Vindiola for useful feedback throughout this project.

Funding

J.E.S. is a CIFAR Fellow in the program Fungal Kingdom: Threats and Opportunities and partially supported by USDA Agriculture Experimental Station at the University of California, Riverside and NIFA Hatch projects CA-R-PPA-5062-H. This work was supported by the Pierce's Disease Control program (sponsor award #: 14–0379-000-SA-2) to F.J.B. and R.R., a California Department of Food and Agriculture (CDFA) agreement # 007011–003 to R.R. and F.J.B., and a CDFA agreement # 20–0267 to P.W.A., L.L.W., R.R, and J.E.S. The funders had no role in the design of the study, collection, analysis, data interpretation or writing of the manuscript.

Author information

Authors and Affiliations

  1. Department of Microbiology and Plant Pathology, University of California, Riverside, Riverside, CA, USA

    Cassandra L. Ettinger & Jason E. Stajich

  2. Department of Entomology, University of California, Riverside, Riverside, CA, USA

    Frank J. Byrne, Inaiara de Souza Pacheco, Dylan J. Brown, Peter W. Atkinson & Richard A. Redak

  3. Department of Botany and Plant Sciences, University of California, Riverside, Riverside, CA, USA

    Linda L. Walling

  4. Institute for Integrative Genome Biology, University of California, Riverside, Riverside, CA, USA

    Linda L. Walling, Peter W. Atkinson & Jason E. Stajich

Authors
  1. Cassandra L. Ettinger
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Frank J. Byrne
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Inaiara de Souza Pacheco
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Dylan J. Brown
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Linda L. Walling
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Peter W. Atkinson
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Richard A. Redak
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Jason E. Stajich
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

C.E. wrote the manuscript text, prepared figures, and analyzed data. J.E. supervised, performed preliminary analysis, and deposited data to public servers. P.A. and I.D.S.P. performed preliminary molecular work. F.B. collected samples and performed resistance assays. D.B. assisted with resistance assays. R.R., F.B. and J.E. conceptualized the work. L.W., R.R., J.E., P.A., and F.B. obtained funding. All authors discussed the results, edited, and reviewed the manuscript. The author(s) read and approved the final manuscript.

Corresponding authors

Correspondence to Cassandra L. Ettinger or Jason E. Stajich.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table s1. 

Differentially expressed genes between insecticide-resistant and susceptible GWSS. DESeq2 was used to identify differentially expressed genes (DEGs) between insecticide-resistant and susceptible GWSS populations. We provide here a table of all significant 607 DEGs including their gene ID, whether the gene was upregulated in resistant or susceptible populations, the log2 fold change, Bonferonni corrected p-value, InterPro annotation, and associated GO terms. DEGs in the table are ordered by p-value.

Additional file 2:

Figure s1. Phylogeny of all cytochrome P450s and their relative expression levels. The eight cytochrome P450s that were differentially expressed between insecticide-resistant and susceptible glassy winged sharpshooters, and which had higher expression in resistant individuals, are highlighted in the tree in orange. A heatmap displays the variance stabilized counts for each cytochrome P450 across all sharpshooters sampled.

Additional file 3: Figure s2. 

STRUCTURE results also indicate no broad-scale population structure. STRUCTURE plots for K=1 to K=4 populations. Marginal likelihood is maximized at K = 1.

Additional file 4: Figure s3. 

LEA results support that overall population structure is indicative of weak differentiation. (A) LEA cross-entropy across K=1 to K=30. There is a break in the cross-entropy values at K=16 which is equal to the number of individuals in this study. Cross-entropy was lowest at K = 1. (B) LEA sparse nonnegative matrix factorization predicted ancestry proportions for K=1 to K=4.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Ettinger, C.L., Byrne, F.J., de Souza Pacheco, I. et al. Transcriptome and population structure of glassy-winged sharpshooters (Homalodisca vitripennis) with varying insecticide resistance in southern California. BMC Genomics 23, 721 (2022). https://doi.org/10.1186/s12864-022-08939-1

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-022-08939-1

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us