RNAseq Analyses
Adapters were trimmed from the raw mRNAseq reads and low-quality reads
were removed using fastqc (Andrews 2010); an average of 10,869,280 reads
per sample were retained (Suppl. Table 1). Sequence reads are deposited
in the NCBI SRA under the BioProject Accessions: PRJNA753142,
PRJNA756317. We took advantage of the well-annotated, published genome
for T. castaneum to align trimmed reads to the most recentT. castaneum assembly (Tcas5.2) using BWA-MEM (avg. reads mapped
= 84.89%, s.d. = 0.023; Suppl. Table 1) . Gene count files were
generated using samtools . We used DESeq2 to perform differential
expression analyses.
First, to investigate the effects of Bt exposure and pesticide selection
regime on gene expression in the absence of pesticide exposure, we
subset only pesticide-free and pesticide-free + Bt libraries and
examined the effects of selection regime, Bt exposure, and their
interaction on gene expression (no pesticide model form: count matrix
~ Reg + BtTx + Reg:BtTx). Next, we analyzed the main
effects of pesticide exposure, pesticide selection regime, and Bt
exposure and the interactive effects between selection regime and
pesticide exposure, selection regime and Bt exposure, and pesticide and
Bt exposure for OP and Pyr separately (OP and Pyr model form: count
matrix ~ PTx + Reg + BtTx + Reg:PTx + Reg:BtTx +
PTx:BtTx). For each factor, the control selection regime, control
pesticide treatment, or control Bt treatment were set as the baseline
comparisons.
Differentially expressed transcripts were annotated using the EnsemblT. castaneum database and the R package ‘biomaRt’ . We manually
curated the annotated gene lists to categorize differentially expressed
transcripts into groups potentially impacted by pesticide resistance and
exposure and Bt treatment such as those involved in detoxification,
immunity, development, and cuticle associated transcripts (Suppl. Table
2). Gene ontology (GO) enrichment analyses for the subsets of
significantly differentially expressed gene sets were performed using
the online Gene Ontology Resource .
To analyze global patterns of gene expression, we used a weighted gene
co-expression network analysis with the R package ‘wgcna’ . This
approach allows us to group genes with correlated expression patterns
into modules and statistically associate gene module expression with
experimental factors. We constructed the initial network using all
replicate samples for all treatments (n = 126) with a variance
stabilizing transformation. A soft thresholding power of 6 was used to
construct the topological overlap matrix, and we then identified modules
with highly correlated module expression patterns (p < 0.05).