*3.5. RNAseq Analysis*

RNAseq analysis was carried out on a minimum of three biological replicates from the five resistant strains and the three laboratory susceptible colonies, Kisumu, N'Goussu and Moz. The correlation matrix shows high degrees of similarity between the two *An. arabiensis* populations, Gaoua-ara and Moz but no clear segregation according to species for the *An. gambiae* and *An. coluzzii* strains (Supplementary Figure S2). Hence, for all further analysis of differential expression between resistant and susceptible strains, Gaoua-ara was compared to Moz alone whereas the three *An. coluzzii* (Tiefora, VK72014, Banfora) and one *An. gambiae* s.s.(Bakaridjan) resistant strains were compared to the average values from the two susceptible strains of *An. coluzzii* (N'Gousso) and *An. gambiae* (Kisumu).

**Figure 6.** Frequency of point mutations associated with resistance. Data reported from samples genotyped in 2019. 995L, 1575N, 119G, 269A and 114I indicate the wildtype genotype (black bars); 995F, 995S, 1570Y, 119S, 296S and 114T indicate resistant genotype (green or purple bars). Heterozygote genotypes are shown with pink bars.

#### *3.6. Similarities between Strains*

The total number of genes differentially expressed across all the resistant compared to susceptible strains is shown in Supplementary Table S3. A total of 81 transcripts were up regulated in resistant versus susceptible strains with 73 down regulated. The upregulated transcripts show no enrichment but two P450s known to bind and/or metabolise pyrethroids (CYP6P3 and CYP6Z2) are amongs<sup>t</sup> the most highly upregulated genes and two glucoronosyl transfearses (UGT302H2 and UGT306A2) are also elevated in all strains.

Down regulated transcripts are strongly enriched for RNA processing (*p* = 1.25 × <sup>10</sup>−4) and do not contain genes previously associated with pyrethroid resistance.

GO term enrichment was explored for each individual resistant population. Whilst no GO terms were enriched across all five resistant populations, a number of similarities were seen across the four resistant *An. gambiae* and *An. coluzzii* colonies (Supplementary Figure S3). GO terms significant in up-regulated genes across each population include oxidoreductase activity, typically seen in resistant colonies [34,35] and related to cytochrome p450 activity, and terms related to neuronal signalling, potentially indicating changes in signalling and neurotransmitter activity are associated with resistance to these neurotoxic insecticides. Additionally, terms related to ATPase activity and GPCR signalling, both previously linked to insecticide resistance [36,37] are seen. There are similarities in GO enrichments in the down-regulated subset of genes, with terms related to transcription factor activity, translational regulation, regulation of dephosphorylation and phosphatase complexes, all repressed (Supplementary Figure S4).

#### *3.7. Differences between Strains*

The RNAseq data was then interrogated to identify both pathways and *a priori* candidate genes enriched in the up or down regulated genes in each resistant strain with *An. gambiae* and *An. coluzzii* compared to two susceptible controls. Analysis at the individual gene level revealed key differences between the strains. For example, 23 P450s are differentially expressed in one or more strains; as mentioned above, CYP6P3 and CYP6Z2 are up-regulated in all resistant strains but other known pyrethroid metabolisers including CYP6M2, CYP6P2, P4 and P5 and CYP9J5 and 9K1 [38,39] are also up-regulated in two or more strains (Figure 7). This analysis also identifies a number of additional P450s that are highly up-regulated in multiple strains but have not ye<sup>t</sup> been functionally characterised (e.g., CYP4H genes) which merit further study.

**Figure 7.** Heatmap showing cytochrome P450 genes that are significantly differentially expressed between the pyrethroid resistant and susceptible strains.

Recently, a number of genes with putative roles in sequestering pyrethroids were found to be over expressed in pyrethroid resistant populations from West Africa [40].

RNAseq data from the *An. gambiae* complex in Burkina Faso is supportive of a putative role of hexamerins in pyrethroid resistance in *An. arabiensis* and the VK7 strain of *An. coluzzii* (as shown previously [37]) (Figure 8). Suppression of the hexamerin AGAP001659 (highly upregulated in Gaoua-ara in this study) was previously associated with a reduction in pyrethroid resistance [37]. In addition, several alpha- cyrstallins were up-regulated in one or more of the pyrethroid resistance populations, with this gene family particularly enriched in the Banfora strain of *An. coluzzii* in agreemen<sup>t</sup> with earlier qPCR data [40]. Suppression of the alpha-cyrstallin AGAP007159, which is upregulated in multiple Burkina populations, has also been shown to result in a reduction in the resistance phenotype in VK7 2014.

**Figure 8.** Heatmap showing differential expression of genes in families putatively associated with insecticide sequestration between the pyrethroid resistant and susceptible strains.

Finally, we looked at expression of genes recently implicated in the cuticular hydrocarbon (CHC) synthetic pathway. This gene list was derived from transcripts encoding the six gene families (propionyl co A synthases, fatty acid synthases, elongases, desaturases, reductases and P450 decarbonylases) that are enriched in the sub epidermal oenocyte cells responsible for CHC production [41]. Several genes in this pathway were up-regulated in the Banfora, Bakaridjan and Gaoua-ara strains but, surprisingly, down-regulated in two of

the strains, Tiefora and VK7 2014 (Supplementary Figure S5). To date only two genes in this putative pathway have been functionally validated, CYP4G16 [42] and the fatty acid synthase FAS1899 [41]; both of these genes are upregulated in the pyrethroid resistant *An. arabiensis* strain (fold changes of 5.2 and 2- fold respectively) suggesting that cuticular resistance may be a particularly important resistance phenotype in this population.

#### *3.8. Evaluation of a Multiplex Gene Expression Assay for Metabolic Resistance*

RNAseq analysis provided a list of putative genes and pathways potentially contributing to the pyrethroid resistance phenotype in the different strains. However, simpler robust assays of gene expression are needed to further investigate the association between gene expression and resistance phenotype. To this end, the Taqman multiplex assay [32] was used to quantify relative expression of a subset of 8 insecticide detoxification genes in each of the resistant strains compared to their susceptible counterparts (to facilitate correlations with RNAseq data, expression levels from the *An. gambiae* and *An. coluzzii* resistant strains were compared to the average expression of the equivalent transcripts in the *An. gambiae* and *An. coluzzii* susceptible strains). The data generated in this study agreed well with previous Taqman multiplex P450 expression data for VK7, with the exception of CYP9K1 (where significant up-regulation was not detected in earlier generations). P450 levels in Banfora appear more variable between generations, consistent with recent findings that the resistance phenotype is less stable in this population than in other laboratory colonies [28]. Within the current study, there is generally good agreemen<sup>t</sup> between the qPCR (Supplementary Figure S6 and Supplementary Table S4) and RNAseq data, with the exception of CYP6P3 and CYP6Z1 (Table 2).


**Table 2.** Summary of correlation between results of detoxification multiplex qPCR and RNAseq data.
