The German cockroach, Blattella germanica L. is an international urban pest that affects millions of residences on a global scale (Vargo, 2021). B. germanica impacts human health through the production of asthma and rhinitis-causing allergens, transmission of food-borne pathogens and psychological stress (Kopanic et al., 1994;Elgderi et al., 2006;Sohn and Kim, 2012). Up to 85% of inner city homes in the United States test positive for cockroach allergens and 60-93% of inner-city children with asthma are sensitized to cockroaches (Gore and Schal, 2007;Do et al., 2016;Pomés et al., 2017). Cockroaches also host pro-and eukaryotic microbes that contribute to house-dust microbiomes that intensify asthma (Roth and Willis, 1960;Pai et al., 2003;Carrasco et al., 2014;Crawford et al., 2015;Thorne et al., 2015;Wada-Katsumata et al., 2015;Lai, 2017;Turturice et al., 2017).
Insecticides are essential for efficiently overcoming health impacts of cockroaches (Schal and Hamilton, 1990;Pomés et al., 2017). However, insecticide resistance has been a formidable recurring barrier to effective cockroach control for decades (Scharf and Gondhalekar, 2021). As of 2016, the German cockroach was reported as having developed resistance worldwide to 42 distinct insecticide active ingredients in at least 219 documented cases (Zhu et al., 2016). Because cockroaches live in relatively closed populations (Crissman et al., 2010;Vargo et al., 2014;Vargo, 2021), resistance can build quickly, even with moderate insecticide selection pressure (Scharf et al., 1997a,b;Gondhalekar et al., 2013;Fardisi et al., 2019). Cockroach baits are widely used in management programs and have been an effective tool for controlling cockroaches and reducing pesticide loads in urban housing (e.g., Miller and Smith, 2020); however, resistance can readily develop even to bait insecticides (Gondhalekar et al., 2011(Gondhalekar et al., , 2013(Gondhalekar et al., , 2016;;Gondhalekar andScharf, 2012, 2013;Ko et al., 2016;Fardisi et al., 2019).
The goal of this research was to use a quantitative metatranscriptomics approach to better understand whole-body physiological responses in B. germanica to defined insecticide selection pressures. The insecticide indoxacarb was used as the selecting insecticide, which is an important bait active ingredient for cockroach control (Appel, 2003;Buczkowski et al., 2008;Gondhalekar et al., 2011). Through previous studies we documented early stages of resistance evolution to indoxacarb among field populations (Gondhalekar et al., 2011(Gondhalekar et al., , 2013) ) and verified hydrolysis and cytochrome P450-based oxidation as important steps in indoxacarb bioactivation and detoxification (Gondhalekar et al., 2016). Our specific objectives here were to: (1) identify whole-body mRNA expression profiles and candidate genes associated with indoxacarb resistance, and (2) investigate a subset of candidate resistance-associated genes in an independent and highly resistant field strain. Our findings reveal novel physiological insights into insecticide resistance evolution in this important global health pest; mainly that resistance evolves rapidly as a complex phenotype with multiple underlying mechanisms that include both xenobiotic detoxification and microbial clearance.
The Arbor Park (AP) strain was used for indoxacarb selection experiments within 12 months of its collection and is the source of the "Parental-F0, " "Selected-F6, " and "Control-F6" lines (Figure 1A). The AP strain was collected from an apartment in Gainesville, FL, United States after control failures with multiple insecticide products. The highly resistant field-collected "Oviedo-R" strain and the standard susceptible Johnson Wax (JWax-S) strain were included in post hoc validation experiments. The Oviedo-R strain was collected from a restaurant near Orlando, FL, United States, where indoxacarb-containing cockroach baits were used regularly for at least 3-4 years and was tested within 6 months of its collection. The JWax-S strain has been in culture for over 70 years and has never been exposed to synthetic organic insecticides. All of the above cockroach strains were reared in mixed life stages in 1,000 cm × 400 cm × 300 cm (L × W × H) plastic boxes with aerated lids, greased walls, cardboard harborage, and an ad libitum water source and rodent diet (#8604, Harlan Teklad, Madison, WI, United States). Rearing conditions consisted of a 12:12 (L:D) photocycle, 25-27 • C and ∼50% RH.
Technical grade indoxacarb (99.1% AI), blank bait matrix and formulated gel bait product containing indoxacarb (Advion R ) were provided by DuPont Inc. (Wilmington, DE, United States). All solvents and buffer components were purchased from Fisher Scientific (Pittsburgh, PA, United States) or Sigma (St. Louis, MO, United States). Other chemicals used in enzyme assays and other procedures are detailed in a previous report (Gondhalekar, 2011). These chemicals include the P450 substrate p-Nitroanisole, the esterase substrates p-nitrophenyl acetate and naphthyl acetate, the GST substrate chloro-dinitro benzene, protein extraction buffers, and native PAGE reagents and stains.
For selection experiments, >1,200 large nymphs (4th to 5th instar) were separated from the AP lab cultures and divided into six sub populations of ca. 200 nymphs. These life stages were used because they are among the most tolerant cockroach life stages (Koehler et al., 1993). A feeding delivery method was used for selections because it exactly represents field exposure to indoxacarb baits (Gondhalekar et al., 2011(Gondhalekar et al., , 2013;;Gondhalekar and Scharf, 2012). In brief, pellets of blank gel bait matrix (ca. 5-10 mg wet weight) were prepared manually and treated with a dose of indoxacarb in acetone that provided ca. 60 to 80% mortality (Supplementary Table 1). Large nymphs that were prestarved for 24-h were held individually with a single indoxacarb treated pellet in 1 oz. (30 mL) cups with vented lids. After 3 days, nymphs that had completely eaten bait pellets were transferred in groups of 100 into 17.8 cm × 17.8 cm × 6 cm disposable plastic Glad R boxes under conditions detailed above, where they remained for 7-10 days. Almost all nymphs exhibited intoxication symptoms at 3 days; however, ca. 25-35% completely recovered over the next 7-10 days. These surviving individuals were reared to adulthood and used as founders for the next generation (Scharf et al., 1997b). Selections continued in this manner for six generations (Parental or F0 to F5 generation). The above process happened with three replicate "selected" and "control" lines; control lines received lab diet only. Details of indoxacarb doses used and percent survival for each round of selection are given in Supplementary Table 1.
Bioassays performed included vial and feeding bioassays, each in two formats.
Vial bioassays were done in concentration-response and diagnostic-concentration formats. Concentration-response bioassays were done with adult males from the Parental (F0), Control F6 and Selected F6 generations following established protocols (Gondhalekar et al., 2011;Fardisi et al., 2019). Four to five concentrations providing 10-90% mortality were used for calculating lethal concentration (LC) estimates and associated parameters using probit analysis (see below). Diagnostic concentration bioassays were done by testing adult males of the JWax-S and Oviedo-R strains at previouslyestablished concentrations of 30 and 60 µg indoxacarb per vial (Gondhalekar et al., 2011(Gondhalekar et al., , 2013)). Three replicates of ten insects were conducted per concentration and isolate (n = 90). Control vials received acetone only.
Feeding bioassays were done in dose-response and no-choice formats. Dose-response bioassays were conducted with Parental (F0), Control F6 and Selected F6 generation adult males using a published protocol (Gondhalekar et al., 2011). At least four doses producing 10-90% mortality were tested against each subpopulation. Each replicate had ten insects and bioassays were repeated three times for each dose and sub-population. The dose-response data were analyzed by probit analysis (details below) and used for calculating lethal doses (LD) and associated parameters. Realized heritability (h 2 ) estimation was done on oral dose-response data to determine the proportion of phenotypic variance in resistance caused by additive genetic variation (Falconer, 1989;Tabashnik, 1992).
No-choice assays used formulated indoxacarb gel bait product (0.6% indoxacarb) with adult males of JWax-S strain and the field-selected highly resistant Mid-Florida strain using published protocols (Gondhalekar et al., 2011). The same protocols were followed for control treatments that were conducted using blank gel bait without indoxacarb. These tests were done on lab-reared individuals within 6 months of collecting the Oviedo-R strain. In these tests no alternative/competing food other than gel bait was present in the bioassay arenas. No-choice bait feeding assays were preferred for testing the Oviedo-R strain because within 6 months after field collection the population numbers were relatively low and the bait feeding bioassay required less insects as compared to the traditional dose or concentration-response tests. Disposable plastic GladWare boxes (17.8 cm × 17.8 cm × 6 cm; Clorox Co., Oakland, CA, United States) served as bioassay arenas and were provisioned with a water source, harborage and 0.5 g of formulated indoxacarb gel bait product in a plastic dish (Gondhalekar et al., 2011). Additional bait or blank matrix was provided if the insects consumed a majority of the initially provided bait. Mortality was recorded at 1, 3, 7, 14, and 21 days. Five replicates with 10 adult males per replicate (n = 50) were performed for each strain.
Enzyme assay methods were detailed previously (Gondhalekar, 2011). Esterase (hydrolysis) and P450 (O-demethylation) activity assays were done on F6-selected and control lines using the model substrates p-nitrophenyl acetate and p-nitroanisole, respectively. Esterase native PAGE was done using the model substrate betanaphthyl acetate. Esterase and P450 investigations were done on soluble and microsomal protein preparations, respectively, made from the same insect homogenates. Protein content was normalized using a commercial Bradford assay (Bio-Rad) with bovine serum albumin as a standard. All enzyme and protein assays were performed in triplicates representing three F6selected and three control lines.
Total RNA was isolated from three independent biological replicate samples of 10 whole adult male cockroaches from each of F6-selected and control lines. A two-step process was used that included the Promega SV Total RNA Isolation Kit (Madison, WI, United States) followed by the Bioline TRIsure kit (Taunton, MA, United States). Manufacturer protocols were followed for both kits with the exception that DNase treatment was excluded as the final step for the second kit. The use of two kits ensured the RNA samples were free of excess protein.
Total RNA yields ranged from 5.8-18.6 µg with A260/280 and 260/230 ratios in the range of 1.8-2.1. Sample quality was further assessed using an Agilent Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States) which verified above yields and provided acceptable RIN scores of 7.20-7.76. RNA samples were enhanced for messenger RNA (mRNA) using the Agilent Tru-Seq RNA prep kit before bar-coded sequencing libraries were made. Sequencing was done on the Illumina HiSeq 2000 platform by the Purdue University Genomics Core (West Lafayette, IN, United States). Sequencing reads were filtered using Phred quality scores and other parameters, and de novo transcriptome assembly performed from all six pooled replicate samples (3 selected and 3 control) using Trinity (Grabherr et al., 2011). Paired reads for individual replicate samples were then mapped to the de novo transcriptome using Bowtie2 (Langmead and Salzberg, 2012), which provided read counts used for differential expression analyses as detailed below. The de novo transcriptome assembly was used because a reference German cockroach genome was not yet available at the time sequencing was completed; however, new blast searches with significant contigs were performed in 2021 which confirmed origins in either the B. germanica genome or from other microbial sources.
The assembled contiguous sequences, i.e., "contigs" were analyzed by BLAST, Blast2GO and KAAS to assign identities and functional annotations. The contigs, as well as singleread "singletons, " were annotated using "Blast2GO" for cellular location (CL), biological process (BP), and molecular function (MF) (Conesa et al., 2005). BLAST searches were performed against the Genbank "nr" database available as of May 2012 at www.ncbi.nlm.nih.gov and re-verified in June 2021 (file provided). KAAS analysis was done to gain insights into possible pathways and gene networks involved in resistance (Moriya et al., 2007). KAAS is a rapid method that establishes orthologies to genes operating within conserved pathways using best hit information and Smith-Waterman scores.
Basic exploration of the data such as accessing data range, library sizes etc. was performed to ensure data quality. Three models were used for analysis of read-count data obtained for each contig and singleton: edgeR (v 2.9), DESeq (v 1.8.3), and voomlimma (v 1.2.0). An edgeR object was created by combining the counts matrix, library sizes, and experimental design (3 replicates each for selected and control lines, i.e., "samples") using the edgeR package. Normalization factors were calculated for the counts matrix, followed by estimation of common dispersion of counts. An Exact test for differences between the negative binomial distribution of counts for the selected and control replicates resulted in finding differential expression, which was then adjusted for multiple-hypothesis testing to generate a result file. A DESeq object analogous to the aforementioned edgeR object was created and used to generate normalization factors followed by dispersion estimates using DESeq package. The DESeq method tests for differences between the base means of the experimental conditions and differential expression (DE) results were reported in another result file (DE_analysis_DESeq.csv). A third method called 'voom' from the limma package was also used for DE analysis. The 'voom' function carries out log2 transformation of counts followed by mean-variance estimation and assigns weight to each transformed value. Linear model coefficients were then calculated using limma's design matrix and log2 transformed values. The linear model was fitted using an empirical Bayes method and differences between counts between selected and control replicates were calculated, which was then adjusted for multiple-hypothesis testing, and reported as result file. Venn diagrams were generated displaying the DE contigs with false discovery rate (FDR) P < 0.05 that were found to be common among all three analysis methods using the online tool Vennyfoot_2 .
Expression levels of 48 significantly up-and down-regulated genes identified from Illumina sequencing (Table 1) were investigated for validative purposes by qRT-PCR in the F6selected and control lines (three biological replicates each). PCR primers for target and reference genes were designed using Primer 3 2 (Untergasser et al., 2012) and are given in Supplementary Table 2. The efficiencies of qPCR primers used in this experiment were empirically determined and they were within the recommended range of 90-110%. Validative qRT-PCR was done on aliquots of the same RNA preparations used for Illumina sequencing above using an iCycler iQ real-time PCR detection system (Bio-Rad, Hercules, CA, United States) with Sybr Green product tagging (2x SensiMix Sybr and Fluorescein Kit; Quantace, Norwood, MA, United States). Each 20 µL qRT-PCR reaction in a 96-well format consisted of 10 µL SensiMix (Bioline, Taunton, MA, United States), 7 µL nanopure water, 1 µL each of forward and reverse primer (0.5 µM final concentration) and 1 µL cDNA. A published qRT-PCR temperature program was followed (Scharf et al., 2008). Three technical replicates were performed for each gene and cDNA preparation. The resulting critical threshold (Cq) data were analyzed by the 2 -CT CT method (Livak and Schmittgen, 2001).
Expression of a subset of 21 genes (indicated by asterisks * in Table 1 and Supplementary Table 2) that showed significant differential expression in transcriptome analyses was also quantified in the lab-susceptible JWax-S and resistant Oviedo-R strains. These genes included 13 host cockroach genes (8 P450s, 2 carboxylesterases, 1 chitinase, 1 transposable element, and 1 hypothetical protein) and 8 genes from eukaryotic/viral microbiota (4 virus, 2 gregarine, 1 coccidia, and 1 unknown). Total RNA was extracted from 1 to 2 weeks old adult males in three replicate groups of 10 for the above two strains. RNA extractions were done using the two-step process described under "Transcriptome Sequencing." cDNA synthesis was done from 500 ng total RNA using the iScript TM cDNA synthesis kit (Bio-Rad, Hercules, CA, United States). All qPCR procedures and expression analysis methods were similar to those mentioned above.
Probit analysis was used for LD and LC determinations using SAS software Version 9.2. Mortality data were corrected for control mortality (<10%) using Abbott's transformation prior to conducting Probit analysis. Resistance ratios (RRs) were calculated by dividing LC or LD estimates for the selected lines by corresponding values for susceptible lines. Significance of RRs was determined according to the procedure outlined by Robertson and Preisler (1992). Mean-separation analyses of insecticide toxicity (vial diagnostic and no choice bait feeding bioassays) were done by ANOVA and paired t-tests (P < 0.05) after arcsine transformation of raw mortality data. For enzyme activity assays, specific activity values for the F6-selected and susceptible sub-populations were compared by non-parametric Mann-Whitney U tests (P < 0.05). Methods used for differential gene expression analysis of read-count data from transcriptome sequencing are explained in a separate section above. Regression analyses were done using JMP Pro 15 software (SAS Institute, Cary, NC, United States) to compare (i) qRT-PCR results to Illumina read counts for the F6-selected and control lines (n = 48 2 https://bioinfo.ut.ee/primer3-0.4.0/ genes), and (ii) the F6-selected and control lines to each other and the Oviedo-R and JWax-S strains (n = 21 genes).
Six generations of selection of the parental AP strain with indoxacarb bait (Figure 1A) resulted in a selected strain having significant levels of resistance in both surface contact (23.8×; Figure 1B) and feeding bioassays (4.5×; Supplementary Table 2). The control strain left to inbreed without selection over the same timeframe acquired no resistance. The selection process led to heritable genetic changes with a realized heritability estimate (h2) of 0.28 (Supplementary Table 3). Biochemical assays on the F6-selected and control lines revealed >2× elevated esterase activity (P = 0.0004; Figure 1C) and ∼0.5× decreased P450 O-demethylation activity (P = 0.0004; Figure 1D).
From six total libraries representing three replicates each of the selected and control lines, 133 million paired-end sequence reads were obtained that contained >1.3 billion total nucleotide bases having an average read length of 98 base pairs (bp). The resulting overlapping sequences were assembled into 207,672 contiguous sequences, hereafter referred to as "contigs." Sequence reads are deposited in the NCBI GEO archive under accession number GSE188950.
Contig length ranged from 201 to 30,113 bp with an average length of 1,777 bp and a N50 size of 4,805 bp; i.e., half of the assembled metatranscriptome is covered by contigs ≥ 4,805 bp. Next, a subset of 48 assembled contigs was chosen for validation analysis by qRT-PCR using the same RNA preparations as were used for Illumina sequencing (Table 1). A regression plot of log2 transformed Illumina transcriptome read count (X) vs. log2 transformed qRT-PCR Cq values (Y) revealed a highly significant correlation (P < 0.0001, r 2 = 0.51) (Supplementary Figure 1). The latter result independently verifies the accuracy of the metatranscriptome results.
DE analysis by three different models yielded differing numbers of passing contigs (Figure 2). Different FDR p-values were considered (P < 0.05, 0.01, 0.001) but the greatest emphasis is placed on the P < 0.001 level. The most stringent analysis model was edgeR with 473 passing contigs (Figure 2A), followed by voom-limma with 1,089 (Figure 2B) and DESeq with 2,209 (Figure 2C). Another interesting feature of the datasets for all three models was the higher ratio of downregulated:upregulated contigs in the selected strain; i.e., 424:49 for edgeR, 960:129 for voom-limma, and 2,083:126 for DESeq. Overall, there were 236 significant DE contigs shared among all three analysis models at the FDR P < 0.01 level and log2FC of ±1 (Figure 2D). Only contigs passing in two or more models at the FDR < 0.001 log2FC ±1 level were considered for further analysis as detailed below.
BLAST2GO analyses for gene annotation revealed that 36.6% of contigs were annotated (76,053). For P < 0.001 passing contigs there were 779 GO terms for cellular location (CL; 99 upregulated, 680 downregulated), 1,514 for molecular function (MF; 209 up, 1,305 down) and 1,493 for biological process (BP; 265 up, 1,288 down) (Supplementary Figure 2A). CL terms potentially related to insecticide resistance include membrane, ribosome, mitochondria, microsome, endoplasmic reticulum and dendrite (Supplementary Figures 2B-D). Potentially relevant MF terms include electron carrier, hydrolase, monooxygenase, oxidoreductase, transferase and catalytic activity; and ATP, calcium, iron, heme, nucleotide and sugar binding. Lastly, relevant BP terms potentially linked to resistance include metabolic, transport, phosphorylation, glutathione conjugation and response to drug. KAAS analysis was used to gain further insights into responsive pathways and gene networks (Supplementary Figure 3). Upregulated pathways potentially linked to resistance
The blastx identities of many upregulated transcripts had logical links to resistance and detoxification. Upregulated cockroach genes included many gene families commonly associated with insecticide resistance including P450 oxidases and hydrolases (Figure 3) and glutathione-S-transferases (GSTs). Table 1 overviews the 21 contigs (shown by asterisks * ) used for post hoc validations and gives a general overview of upregulated transcripts that mainly match the B. germanica genome, along with downregulated transcripts that originate mainly from microbial sources. The identities of sequences from either cockroach or microbial origins were re-verified by Genbank blastX searches performed in September 2021 (Table 1).
Taxonomic compositions of significant DE contigs (P < 0.001), based on blastx database queries, indicate the majority of upregulated contigs are from the cockroach genome; whereas, the majority of downregulated contigs come from viruses and eukaryotic microbes (Figure 4). Top upregulated taxonomic matches at the domain level were overwhelmingly eukaryotic whereas downregulated contigs had >10-fold more numbers of matches to bacteria, viruses and archaea than insects. At the genera level the top taxonomic matches for upregulated contigs were nearly all insects (Blattella, Cryptotermes, Zootermopsis, Coptotermes, Timema, and Periplaneta). Downregulated transcripts at the genus level were dominated by eukaryotic microbes that are apparently commensal, pathogenic and/or parasitic (Gregarina, Cryptosporidium, Toxoplasma, Vitrella, Plasmodium, Neospora, etc.).
An independent resistant strain, Oviedo-R, was collected from the field after indoxacarb bait control failures and assayed alongside the standard susceptible JWax-S strain. The Oviedo-R strain was highly resistant with 0% mortality in vial bioassays on diagnostic concentrations of indoxacarb (Figure 5A). These diagnostic concentrations were approximately two and fourfold higher than the indoxacarb LC50 for the Control (F6) strain shown in Figure 1B. The JWax-S strain was highly susceptible with 100% mortality in the same assays. When tested in no-choice feeding bioassays with commercially-formulated indoxacarb bait, the Oviedo-R strain displayed exceptionally high resistance (ANOVA df = 1,48, F = 130.50, P < 0.001) (Figure 5B). The Oviedo-R strain entirely consumed 0.5 g indoxacarb bait per assay replicate by Days 7 and 14, at which time the bait was replenished. These results show high levels of physiological resistance to indoxacarb in the Oviedo-R strain with no involvement of a behavioral "aversion" component.
Finally, to test for common patterns of gene expression across strains, qRT-PCR analyses were performed on a subset of 21 significant up-and down-regulated contigs identified from the transcriptome analysis presented above (Table 1 and Supplementary Table 2). The four strains included in this analysis were Oviedo-R, JWax-S, indoxacarb Selected-F6 and Control-F6. From a series of regression analyses, there was significantly correlated transcript abundance between the Selected-F6 and Oviedo-R strains (Figure 6A), but no correlation when comparing Oviedo-R vs. Control-F6 (Figure 6B), Control-F6 vs. Selected-F6 (Figure 6C), and JWax-S vs. Selected F6 (Figure 6D). These results suggest common processes associated with indoxacarb resistance evolution in lab and field-selected cockroach populations.
This study reveals new insights into cockroach insecticide resistance evolution from multiple physiological perspectives. Through a combination of approaches including selection for FIGURE 5 | Indoxacarb bioassay results for the highly resistant Oviedo-R strain and the standard susceptible JWax-S strain. (A) Vial diagnostic concentration bioassays at two concentrations showing high-level mortality and the JWax-S strain and 0% mortality in the Oviedo-R strain. (B) No-choice feeding bioassays using commercial formulated indoxacarb bait showing rapid high-level mortality in the JWax-S strain and virtually no mortality in the Oviedo-R strain. * Asterisks indicate significant differences between strains at P < 0.0001.
resistance via feeding on insecticidal bait, coupled with different bioassay formats and metatranscriptome sequencing, we were able to observe that (1) the host cockroach mainly upregulates a range of detoxification resistance mechanisms, and (2) at the same time decreases its internal virus, parasite and/or pathogen levels. Further investigation of candidate gene expression in the highly resistant Oviedo-R strain suggests common phenomena that underlie resistance evolution. Taken together, these findings support the idea that high-level resistance evolution results from a dual process whereby the host tolerates the selecting insecticide through a number of potential mechanisms, and in parallel, host fitness is further increased as the body is cleared of parasites and pathogens. At present it remains unclear if the microbial impacts result from (i) direct effects of indoxacarb, (ii) indirect effects of antimicrobial preservatives included in the selecting bait matrix, or (iii) selection for general stress response mechanisms that confer both xenobiotic resistance and microbial immunity. With respect to the first two possibilities, it is unclear if these processes are associated generally with all insecticides, or specifically with ingested indoxacarb. Transcriptome sequencing revealed dozens of candidate upregulated genes from the host cockroach potentially involved in resistance, including detoxification enzymes (discussed below), transport mechanisms, host gut penetration barriers, and others. An unanticipated result was the downregulation thereby indicating disappearance of commensal, pathogenic and/or parasitic microbe transcript contigs after insecticide selection (also discussed below). Findings in both categories are corroborated by the GO and KAAS analyses, which provide additional independent confirmation for the sequence composition results relating to detoxification and parasite/pathogen disappearance. The qPCR correlation analyses between resistant and susceptible strains lends further strength to the above findings. Specifically, the regression analyses showed similar trends between independently-selected resistant strains for upregulation and downregulation of host detox mechanisms, and disappearance of a wide range of microbes. However, there appears to be wider variation with respect to disappearance of microbial contigs (right side of Figure 6A). This finding is logical given that different environments have different microbiomes associated with them that can, in turn, shape internal microbiomes of higher organisms living within them (Renelies-Hamilton et al., 2021;Tinker and Ottesen, 2021).
In terms of detoxification and metabolism, indoxacarb is a unique pro-insecticide that requires hydrolytic activation to its decarbomethoxyllated "DCJW" metabolite to become an active insecticide (Wing et al., 2000;Alves et al., 2008). A prior study investigating indoxacarb metabolism pathways in resistant and susceptible B. germanica strains revealed that hydrolysis and P450-oxidation were important to activation and resistance-associated detoxification, respectively (Gondhalekar et al., 2016). The current study reveals candidate genes involved in both oxidative and hydrolytic metabolism of indoxacarb; specifically, 20 differentially expressed (DE) P450 contigs and 14 DE esterase/hydrolase contigs. In theory, resistance could result from both upregulation of P450 oxidative enzymes and downregulation of hydrolases, both of which were detected in the present study.
With respect to detoxification genes and the B. germanica genome, a targeted analysis revealed 158 P450 genes total, which represents some gene expansions in relation to the available genomes of close insect relatives (Harrison et al., 2018a,b). The most expanded P450 families were Cyp4 (n = 59), followed by a subset of Cyp6 (n = 8). The present study identified resistanceassociated changes of expression of 12.6% of the B. germanica "cypome, " with the majority (14) being Cyp6 members and fewer being Cyp4s (4). It is also noteworthy that a transposable element potentially associated with Cyp4 expansion (Harrison et al., 2018b) had higher resistance-associated expression in both the indoxacarb-selected line and the field-selected Oviedo-R strain. Expansions previously identified for other detoxification genes in the B. germanica genome included 62 E4 esterases and a subgroup of 23 GST genes (Harrison et al., 2018b). The present study identified 14 DE esterases (4 up-and 10 down-regulated) and 11 GSTs (6 up-and 5 down-regulated). Significant changes in both P450 and esterase activity were noted with selection (Figure 1) but not GST activity (Gondhalekar, 2011).
Numerous DE P450 contigs identified here were close homologs to Cyp6J1 and 6K1 of B. germanica (Wen and Scott, 2001), but none were identical matches. Some of these Cyp6 homologs had 20-50× upregulation with indoxacarb selection. Based on these results it is likely that there is Cyp6 involvement in the detoxification pathway for indoxacarb (Gondhalekar et al., 2016) and more Cyp6 diversity in the B. germanica genome than initially suggested (Wen and Scott, 2001;Harrison et al., 2018b). A previous study found Cyp4G19 over expression in association with fipronil and imidacloprid resistance in B. germanica, but not indoxacarb resistance, which agrees with results of the present study finding no change of expression for Cyp4G P450s. Cyp4G19 was also linked to pyrethroid resistance, cuticular hydrocarbon production and cuticular penetrationbased resistance (Pridgeon et al., 2003;Guo et al., 2010;Pei et al., 2019;Chen et al., 2020;Hu et al., 2021), which are seemingly unrelated to the ingestion and gut uptake that occurred during the six generations of indoxacarb selection that were done in the present study.
The current study also identified DE homologs of Cyp9e2 and Cyp4C21 with possible roles in indoxacarb biotransformation. Both Cyp9e2 and Cyp4C21 were previously identified from B. germanica along with apparent homologous pseudogenes that may confound gene identification (Wen et al., 2001). Another upregulated P450 in the present study was Cyp15F1, which catalyzes the last step in juvenile hormone biosynthesis (Maestro et al., 2010;Tarver et al., 2012;Zhou et al., 2014). The Cyp15F1 homolog identified in the present study had a strong match to its corresponding genomic sequence reported in the B. germanica genome (Harrison et al., 2018a,b). It is not clear if Cyp15F1 participates in indoxacarb detoxification or another physiological process related to gut tissue remodeling after indoxacarb exposure, or in response to clearance of gut microbes (e.g., Sen et al., 2015). A previous study identified the P450 protein "P450MA" from a multi-resistant B. germanica strain (Scharf et al., 1998). P450MA was overexpressed after organophosphate insecticide selection and also was found to be over-expressed in organophosphate resistant B. germanica strains from different global origins (Scharf et al., 1997b(Scharf et al., , 1999)). While the molecular identity of P450MA remains unknown, it may be among the P450s identified in the present study, of which Cyp6 members were the most abundant.
Lastly, the identification of decreased P450 O-demethylation activity in the F6-selected strain does not appear to be an erroneous or trivial result. This is because there is good agreement between reduced O-demethylation activity in our F6selected indoxacarb-resistant strain and transcriptome results showing downregulation of several P450 contigs in the same genetic lineage. This latter finding also suggests that increased O-demethylation activity may not be a useful biomarker for indoxacarb resistance.
With regard to hydrolytic activity, it is considered more important for indoxacarb activation in B. germanica than detoxification (Gondhalekar et al., 2016). Esterases and associated hydrolases are overall not as well studied as the P450s discussed above, but several DE hydrolase and esterase homologs were identified through the current study. Several of these DE esterase contigs had significant matches to FE4 esterases from other insects with well-established roles in insecticide metabolism (Field and Foster, 2002;Srigiriraju et al., 2009). Esterase activity toward model substrates was also elevated in the F6-selected strain, which is consistent with the several upregulated esterase/hydrolase contigs that were identified through transcriptome sequencing. It is also possible that, despite being important for indoxacarb activation, increased hydrolytic activation could still enable greater detoxification by P450 and other detox enzymes and subsequent clearance from the body (Gondhalekar et al., 2016) by Glutathione-S-transferases, some of which were also upregulated.
Invertebrates were the first hosts of apicomplexan parasites like Plasmodium and Toxoplasma which later switched to vertebrate hosts during their evolution (Kopecná et al., 2006). Biotic associations of cockroaches with microbes have been known for over 100 years and include bacteria, archaea, viruses, protists, fungi, gregarines, and nematodes (Roth and Willis, 1960;Kakumanu et al., 2018). In the present study, transcript contigs from representatives of all of these groups were greatly reduced with indoxacarb selection. Omics approaches similar to those used here have previously identified gregarines and other eukaryotic microbes in insects and thus it should not be surprising that similar eukaryotic microbe transcripts were identified in the present study through the use of poly-A RNA tagging (McCarthy et al., 2011;Scharf, 2015;Scharf et al., 2017). Many of the microbial genera that were reduced by the selection process are potential pathogens to humans and companion animals (e.g., Cryptosporidium, Toxoplasma, Vitrella, Plasmodium, Neospora, etc.) and thus our findings on their disappearance may be highlighting a previously unknown/unacknowledged benefit of cockroach baits.
Gregarine protist contigs were the dominant downregulated contigs in the F6-indoxacarb selected line. One common parasitic gregarine, Gregarine blattarum, can infect multiple cockroach species including B. germanica, as well as other invertebrates (Yahaya et al., 2017). Gregarine infection has been shown to cause pathological effects in B. germanica as well as increase susceptibility to the fungal pathogen Metarhizium anisopliae and the growth regulator insecticide triflumuron (Lopes and Alves, 2005), which compels us to ask the question: is insecticide susceptibility caused in part by microbial pathogens or parasite stress? Cockroaches in culture are particularly susceptible to gregarine infections and such infections can be reduced by antimicrobial drugs (Clopton and Smith, 2002;Smith and Clopton, 2003). However, it does not appear that many of the sampled microbial groups are exclusively associated with laboratory rearing (Kakumanu et al., 2018). Gregarines are further known to accelerate larval development in cat fleas, Ctenocephalides felis (Alarcón et al., 2017) and have recently been identified in cucumber beetles, Acalymma vittatum, with no apparent impacts on fitness (Coco et al., 2020). To our knowledge no prior reports are available describing the effects of insecticides on gregarines, but a prior report does document impacts on related termite gut protists by the nicotinoid insecticide imidacloprid (Sen et al., 2015). While it is not clear if the disappearance of these internal microbes enhances host fitness to enable a rapid buildup of xenobiotic resistance, or vice-versa (i.e., if resistance happens first), our findings provide important new insights into the potential stepwise basis of resistance evolution and the complex physiological interactions involved.
This study provides new insights into cockroach insecticide resistance evolution from multiple physiological perspectives ranging from xenobiotic metabolism and excretion to pathogen and parasite resistance. The identities of dozens of candidate bioactivation and detoxification enzymes were revealed, namely cytochrome P450s and esterases/hydrolases, which agrees strongly with outcomes of preceding indoxacarb metabolomics work in B. germanica (Gondhalekar et al., 2016). Many of the genes studied here also had strong matches to the B. germanica genome (Harrison et al., 2018a,b) and thus our findings provide important annotations that have been lacking in many cases. Because we used a selection-based approach with temporally parallel non-selected controls that originated from the same genetic stock, and compared strains from both common and distinct genetic origins in a stepwise fashion, the insights provided have limited caveats. Thus, the correlated expression for a subset of candidate genes between independently-selected resistant strains suggests there are common, predictable patterns to resistance evolution across populations.
An unanticipated outcome was that numerous microbial transcripts were reduced with insecticide selection. In terms of the causative agents behind microbial disappearance, three possibilities exist: (i) direct antimicrobial effects by indoxacarb, (ii) indirect effects of antimicrobial preservatives in the bait matrix, or (iii) co-selection for dual detoxification and immune pathways. While the potentially direct effects of the selecting insecticide or bait matrix preservatives on eukaryotic microbes are logical, the causative factors underlying virus and bacterial declines are less clear and should be further investigated. Findings also revealed a potential added benefit of cockroach baits for curing cockroaches of potentially deleterious pathogens and associated allergens that can affect both humans and companion animals. A priori goals of this study did not include identification of effects on cockroach parasites and pathogens, and thus our experimental design was not optimized for gaining insights into host-microbial interactions. Future work thus needs to examine for similarities in genes, microbes and processes revealed here, with the ultimate goal of reducing impacts of insecticide resistance and concurrently creating healthier indoor urban environments.
Frontiers in Physiology | www.frontiersin.org
February 2022 | Volume 12 | Article 816675
http://bioinfogp.cnb.csic.es/tools/venny/index.html Frontiers in Physiology | www.frontiersin.org