Transcriptomic analysis of Anopheles gambiae from Benin reveals overexpression of salivary and cuticular proteins associated with cross-resistance to pyrethroids and organophosphates

doi:10.21203/rs.3.rs-3408547/v1

Download PDF

Research Article

Transcriptomic analysis of Anopheles gambiae from Benin reveals overexpression of salivary and cuticular proteins associated with cross-resistance to pyrethroids and organophosphates

https://doi.org/10.21203/rs.3.rs-3408547/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 06 Apr, 2024

Read the published version in BMC Genomics →

You are reading this latest preprint version

Background

Insecticide resistance (IR) is one of the major threats to malaria vector control programs in endemic countries. However, the mechanisms underlying insecticide resistance are poorly understood. Thus, investigating gene expression patterns related to IR can offer important insights into the molecular basis of IR in mosquitoes. In this study, RNA-Seq was used to characterize gene expression in Anopheles gambiae surviving exposure to pyrethroids (deltamethrin, alphacypermethrin) and an organophosphate (pirimiphos-methyl).

Results

Larvae of An. gambiae s.s. collected from Bassila and Djougou in Benin were reared to adulthood and phenotyped for IR using a modified CDC intensity bottle bioassay. The results showed that mosquitoes from Djougou were more resistant to pyrethroids (5X deltamethrin: 51.7% mortality; 2X alphacypermethrin: 47.4%) than Bassila (1X deltamethrin: 70.7%; 1X alphacypermethrin: 77.7%), while the latter were more resistant to pirimiphos-methyl (1.5X: 48.3% in Bassila and 1X: 21.5% in Djougou). RNA-seq was then conducted on resistant mosquitoes, non-exposed mosquitoes from the same locations, and laboratory-susceptible An. gambiae s.s. Kisumu strain. The results showed overexpression of detoxification genes, including cytochrome P450s (CYP12F2, CYP12F3, CYP4H15, CYP4H17, CYP6Z3, CYP9K1, CYP4G16, and CYP4D17), carboxylesterase genes (COEJHE5E, COE22933) and glutathione S-transferases (GSTE2 and GSTMS3) in all three resistant mosquito groups analyzed. Genes encoding cuticular proteins (CPR130, CPR10, CPR15, CPR16, CPR127, CPAP3-C, CPAP3-B, and CPR76) were also overexpressed in all the resistant groups, indicating their potential role in cross resistance in An. gambiae. Salivary gland protein genes related to ‘salivary cysteine-rich peptide’ and ‘salivary secreted mucin 3’ were also over-expressed and shared across all resistant groups.

Conclusion

Our results suggest that in addition to metabolic enzymes, cuticular and salivary gland proteins could play an important role in cross-resistance to multiple classes of insecticides in Benin. These genes warrant further investigation to validate their functional role in An. gambiae resistance to insecticides.

Insecticide resistance

An. gambiae

vector surveillance

RNA-seq

differential gene expression

Insecticide-based vector control approaches such as insecticide treated nets (ITNs) and indoor residual spraying (IRS) are core methods to break human-vector contact, thus reducing transmission of malaria [1]. The insecticides traditionally used in vector control fall within 4 main classes: pyrethroids, carbamates, organochlorines, and organophosphates, with pyrroles and neonicotinoids being recently repurposed from agricultural pesticides [2]. Because of this reliance on insecticide-based control methods, the emergence and spread of insecticide resistance in mosquito vectors is becoming a serious threat.

Insecticide resistance is caused by multiple mechanisms including target site insensitivity, such as knockdown resistance (kdr), and metabolic resistance brought about by the increased production of enzymes capable of breaking down insecticides [3]. Knockdown resistance (kdr) arises from mutations in the voltage-gated sodium channel gene which is a target of pyrethroid and organochlorine insecticides, while the acetylcholinesterase (ace-1) gene is a target of organophosphate and carbamate insecticides. Mutations such as L1014S and L1014F in the kdr gene [4] and the G280S mutation in the ace-1 gene [5] cause structural changes that decrease the ability of certain insecticides to bind with their target sites. Additionally, the overexpression of genes encoding detoxification enzymes (e.g., cytochrome P450s, glutathione S-transferases and carboxylesterases) are known to contribute to metabolic resistance in mosquitoes [6, 7] and can be markers of IR in mosquito populations [8]. Because these changes resulting in IR have a genetic basis, robust sets of molecular markers could provide sensitive and timely diagnoses of insecticide resistance in mosquito populations, thereby enabling the implementation of effective insecticide resistance management (IRM) strategies.

Next-generation sequencing approaches such as RNA-seq and whole genome sequencing (WGS) enable the understanding of the genetic variations that result in changes to mosquito biology and behavior in response to environmental factors, including insecticide pressure. Transcriptomic analyses allow for the investigation of gene expression and polymorphic variations associated with specific phenotypes [9, 10]. Furthermore, it allows the identification of novel candidate genes related to a specific insecticide or multi-insecticide resistance [9, 11, 12]. A recent study of the malaria vector Anopheles arabiensis found multiple highly overexpressed genes related to cuticular-associated proteins and salivary gland proteins associated with pyrethroid and organophosphate resistance, suggesting roles of these lesser-understood gene groups in cross-resistance [12].

The objective of this study was to use transcriptomic data to identify candidate genes associated with insecticide resistance in An. gambiae s.s collected from Djougou and Bassila, two sites in Benin with different levels of resistance to pyrethroid or organophosphate insecticides.

Resistance profile of An. gambiae from Bassila and Djougou

The mortality of the mosquitoes in each insecticide bioassay is presented in Fig. 1 and Additional file 1. There was a significant difference in the mortality rate between mosquitoes from Bassila and mosquitoes from Djougou (P < 0.05). Mosquitoes from Djougou were more resistant to alphacypermethrin and deltamethrin than mosquitoes from Bassila. At 1X exposure to alphacypermethrin and deltamethrin, the average percent mortality in Djougou was zero, while it was 77.5% (sd: 2.28) for alphacypermethrin and 65.5% (sd: 3.03) for deltamethrin in Bassila. At 2X exposure, the average percent mortality in Djougou was 46.7% (sd: 1.06) and 26.68% (sd: 2.5) for alphacypermethrin and deltamethrin, respectively, and in Bassila, it was 82.8% (sd: 0.1) and 80% (sd: 5.54), respectively. At 5X exposure, the average percent mortality in Djougou was 77.4% (sd: 1.65) and 50.8% (sd: 2.98) for alphacypermethrin and deltamethrin, respectively, and 100% and 86.8% (sd: 0.72) in Bassila (Additional file 1).

In contrast to the pyrethroids, mosquitoes from Bassila were more resistant than mosquitoes from Djougou when tested against pirimiphos-methyl. The mortality at 1X exposure was 22.0% (sd: 0.31) in Djougou and 2.1% (sd: 0.58) in Bassila. At 1.5X exposure, the mortality was 75.2% (sd: 1.66) in Djougou and 26.3% (sd: 3.81) in Bassila (Additional file 1).

RNA-Seq data quality control and mapping

Raw reads generated ranged from 45–112 million for mosquitoes from Bassila, 51–108 million for mosquitoes from Djougou and 77–102 million for the susceptible Kisumu strain mosquitoes. After filtering, more than 98% of the reads were retained in all experiments and mapped to the whole reference genome of Anopheles gambiae PEST (VectorBase release 48). The percentage of reads mapped to the reference genome ranged between 61% and 70.6% for Bassila, 51.4% and 76% for Djougou and 69.2% and 75% for the susceptible Kisumu strain (Additional file 2), and 70% and 79% of the alignments (read pairs) were successfully assigned to the exonic features of the gene set AgamP4 (Additional file 3).

Differential gene expression analysis

EdgeR was used to perform differential gene expression (DGE) analysis between the different mosquito groups (R-S), (C-S) and (R-C). A fold-change (FC) > 2 and a false discovery rate (FDR) < 0.01 were used to identify differentially expressed genes (Fig. 2). The DGE results are summarized in Table 1 and Additional file 4.

Table 1

Differential gene expression summary for alphacypermethrin, deltamethrin and pirimiphos-methyl
Insecticide	Site	Condition	Number of genes tested		DE genes (\|FC\|>2 & adjP < 0.05)		DE genes (\|FC\|>2 & adjP < 0.01)
Insecticide	Site	Condition	Number of genes tested		UP	Down	UP	Down
Alphacypermethrin	Djougou	DA vs DU (R-C)		9743	31	30	13	12
Alphacypermethrin	Djougou	DA vs KIS (R-S)		10041	949	860	850	755
Deltamethrin	Bassila	BD vs BU (R-C)		9920	41	65	28	47
	Bassila	BD vs KIS (R-S)		10146	966	819	943	805
	Djougou	DD vs DU (R-C)		9937	68	24	45	17
	Djougou	DD vs KIS (R-S)		10154	887	573	822	555
Pirimiphos-methyl	Bassila	BP vs BU (R-C)		9740	12	50	8	39
	Bassila	BP vs KIS (R-S)		10080	968	776	926	758
	Djougou	DP vs DU (R-C)		9571	94	114	55	81
	Djougou	DP vs KIS (R-S)		9979	953	764	932	750
Unexposed	Bassila	BU vs KIS (C-S)		9999	886	696	863	686
Unexposed	Djougou	DU vs KIS (C-S)		10014	757	609	698	576

Differential gene expression associated with alphacypermethrin.

Genes associated with alphacypermethrin resistance were derived from three different comparisons using mosquitoes from Djougou. A total of 1274 genes (698 upregulated and 576 downregulated), 1605 genes (850 upregulated and 755 downregulated) and 25 genes (13 upregulated and 12 downregulated) were significantly differentially expressed in the C-S, R-S and R-C comparisons, respectively (Table 1). A total of 4 key upregulated genes were shared between the three comparison sets, with only one characterized protein belonging to the cuticular protein RR-1 family 75 (Additional file 5-A). Comparing R-S and C-S comparisons, genes overexpressed in both sets are shown in (Additional file 6). Interestingly, a carboxylesterase (COE22933), some cuticular proteins (CPR76, CPR75, CPR16, TWDL1, CPR81, and CPAP3-A1b), some members of the cytochrome P450 family (CYP4D17, and CYP12F2) and some salivary gland proteins (salivary secreted peptide (AGAP013060), kDa salivary (AGAP004316)) showed higher fold changes in the R-S comparison compared to the C-S comparison. Along with these genes, DE-cadherin-like isoform X1 (AGAP029696), COX2, TEP1, ESP, methyltransferase 22 isoform X1 (AGAP011879), microfibril-associated glyco 4-like (AGAP010531) and BTB/POZ domain-containing protein KCTD16 (AGAP000805) also showed higher fold changes in the R-S set of genes (Additional file 6). Importantly, a comparison between the R-S and R-C groups available in (Additional file 7) showed that the cuticular protein (CPR75) and the carboxylesterases (COEunkn) were differentially expressed and shared among both comparisons.

Focusing on R-S comparisons, carboxylesterases (COE12O, COEunkn, and COEBE3C), cuticular proteins (CPAP3-A1c, CPAP3-A1a, and CPAP3-D), cytochrome P450s (CYP12F3, CYP306A1, and CYP4G16), salivary gland proteins (salivary secreted mucin 3 (AGAP009473), secreted salivary gland (AGAP001989)), and an odorant-binding protein (AGAP012867) were overexpressed and were not significantly expressed in the C-S set of genes (Additional file 6, Fig. 3). Additionally, other genes (FC ~ 11.2-380.5) overexpressed in the R-S set of genes were related to neuronal acetylcholine receptor subunit non-alpha-3-like (AGAP010169), autophagy 12-like (AGAP012847), DE-cadherin-like isoform X1 (AGAP029698), antigen 5 like allergen Cul n 1-like (AGAP006417) or NADH dehydrogenase subunit 6 (AGAP028386), and CEC1 (Additional file 6).

Differential gene expression associated with deltamethrin

In Bassila, 1549 genes (863 upregulated and 686 downregulated), 1748 genes (943 upregulated and 805 downregulated) and 75 genes (28 upregulated and 47 downregulated) were significantly differentially expressed in C-S, R-S and R-C comparisons, respectively (Table 1). No DEGs were shared between the three groups (Additional file 5-B). In Djougou, the number of DEGs significantly differentially expressed in the C-S, R-S and R-C comparisons were 1274 genes (698 upregulated and 576 downregulated), 1377 genes (822 upregulated and 555 downregulated) and 62 genes (45 upregulated and 17 downregulated), respectively (Table 1). Unlike in Bassila, five DEGs were upregulated and shared between the three groups (Additional file 5-C). Among them, two had retrievable annotations and were related to SERAC1 isoform X1 (AGAP011044) and serine protease (CLIPB5).

A total of 1240 DEGs (661 upregulated and 579 downregulated) and 887 DEGs (506 upregulated and 381 downregulated) were shared between the R-S and C-S groups from Bassila and Djougou, respectively. The DEGs in Bassila included a carboxylesterase (COE120), five cuticular proteins (CPCFC1, CPR125, CPR140, and CPLCG5), two cytochrome P450s (CYP9M1, and CYP6P2), two glutathione S-transferases (GSTE2, and GSTE4) and six salivary gland proteins (D7r2, SG7, SG3, SG1c, D7L1, and D7L2) with markedly higher expression in the R-S group compared to the C-S group. In Djougou, 7 cuticular proteins (CPAP3-A1b, CPR76, TWDL1, CPAP3-B, CPR81, CPR75, and CPAP3-E), 3 cytochrome P450s (CYP6AA2, CYP12F2, and CYP6Z3) and a salivary gland protein (SG2) exhibited higher expression in the R-S group than in the C-S group (Additional file 6, Fig. 3). Other DEGs with higher expression in the R-S group than in the C-S group were related to TEP9, GNBPB4, indirect flight muscle (AGAP011514), and an apyrase-like (AGAP005457) in Bassila, and St1, CLIPB5, TEP4, ESP, CLIPC7, and glycine-rich 5-like (AGAP001508), SERAC1 isoform X1 (AGAP011044), and myrosinase 1-like (AGAP006426) in Djougou (Additional file 6).

Comparing the R-S and R-C groups, CYP6Z3 in Bassila, and CPR9, CPR144 and CYP4H24 in Djougou were shared among both comparisons. In addition, the angiopoietin-like salivary protein (AGAP007041) in Bassila and CYP314A1 in Djougou were overexpressed in only the R-C group, together with other genes such as fibrinogen A (AGAP011228), CLIPB12, PPO6, and PGRPLB in Djougou (Additional file 7).

Focusing on DEGs within the R-S group, a total of 951 (546 upregulated and 405 downregulated) genes were commonly differentially expressed in deltamethrin survivors from both the Bassila and Djougou populations relative to the susceptible strain (Fig. 4-A). The upregulated genes included one carboxylesterase (COEJHE5E), seven cuticular proteins (CPLCP3, CPR130, CPR30, CPR59, CPAP3-B, CPR76, and TWDL1), five cytochrome P450s (CYP4H15, CYP12F2, CYP6Z3, CYP9K1, and CYP9M1), three glutathione S-transferases (GSTE2, GSTD7, and GSTE4), a salivary gland protein (SG2), and an odorant binding protein (OBP47). Interestingly, the detoxification genes CYP12F2 and CYP6Z3 exhibited higher transcription activity in deltamethrin survivors from Djougou (FC = 33.6 & 24.3 respectively) than in those from Bassila (FC = 15.1 & 6.7 respectively), suggesting their potential association with the differing levels of resistance between these two field populations (Additional file 8, Fig. 3).

Differential gene expression associated with pirimiphos-methyl

Differentially expressed genes in Bassila included 1549 genes (863 upregulated and 686 downregulated), 1684 genes (926 upregulated and 758 downregulated) and 47 genes (8 upregulated and 39 downregulated) in C-S, R-S and R-C respectively (Table 1). However, in Djougou, 1274 genes (698 upregulated and 576 downregulated), 1682 genes (932 upregulated and 750 downregulated) and 136 genes (55 upregulated and 81 downregulated) were differentially expressed in C-S, R-S and R-C, respectively (Table 1).

In Bassila, a total of 1220 DEGs (665 upregulated and 555 downregulated) were shared between the R-S and C-S groups (Additional file 5-D). Some of them exhibited higher fold change expression in the R-S group than in the C-S group, and these included four cuticular proteins (CPR16, CPR30, CPLCG5, and CPCFC1), six cytochrome P450s (CYP9M1, CYP6Z2, CYP6P3, CYP6Z3, CYP307A1, and CYP12F2), two glutathione S-transferases (GSTS1, and GSTD7) and two salivary gland proteins (SG7, and SG9) (Additional file 6). In Djougou, 982 DEGs (545 upregulated and 437 downregulated) were shared between the R-S and C-S comparisons (Additional file 5-E). From these DEGs, those that exhibited higher fold changes included one carboxylesterase (COE22933), seven cuticular proteins (CPR30, CPR76, CPR16, CPLCG1, CPR81, CPLCX2, and CPAP3-E), seven cytochrome P450s (CYP6M2, CYP6P3, CYP6Z3, CYP4D22, CYP9K1, CYP6AA2, and CYP6Z2), two glutathione S-transferases (GSTE2, and GSTD7), two salivary gland proteins (SG7, and SG9) and one odorant binding protein (OBP25) (Fig. 3). Other DEGs with retrievable annotations were two trypsin-related proteases (TRYP7, and AGAP012842), a glycine-rich cell wall structural-like (AGAP008892) and others (Additional file 6).

Comparing the R-S and R-C groups, no metabolic genes were shared among the two comparisons in Bassila, while CYP6Z3 and GSTD11 were shared in Djougou. Additionally, two carboxylesterases (AGAP011507, COEAE2E), PPO6, PPO9, PGRPLB, and fibrinogen A (AGAP011228) were overexpressed in only the R-C group in Djougou (Additional file 7).

Focusing on the DEGs in the R-S group, a total of 1191 DEGs (665 upregulated and 526 downregulated) were commonly differentially expressed in the pirimiphos-methyl survivors from both Bassila and Djougou when compared to the susceptible strain (Fig. 4-B). These shared DEGs included two carboxylesterases (COE22933, and COEJHE5E), 11 cuticular proteins (CPLCG5, CPR130, CPR30, CPCFC1, CPR76, CPR75, CPAP3-B, CPR127, CPR15, CPR16, and CPR144), seven cytochrome P450s (CYP6Z3, CYP12F2, CYP6M2, CYP9K1, CYP6P4, CYP6P3, and CYP6Z2), three glutathione S-transferases (GSTD7, GSTE2, and GSTE4), four salivary gland proteins (SG7, SG2, SG9, and salivary cysteine-rich peptide (AGAP011460)), and two odorant binding proteins (OBP47, OBP26) (Additional file 8, Fig. 3).

Genes associated with multiple insecticide resistance

A total of 500 DEGs (265 upregulated and 235 downregulated) were shared by mosquitoes that survived exposure to either alphacypermethrin, deltamethrin or pirimiphos-methyl (Fig. 5-A, Additional file 9). Among those genes, 10 were cuticular proteins (CPAP3-B, CPAP3-E, CPR10, CPR130, CPR15, CPR16, CPR30, CPR75, CPR76, and CPR81), 8 were cytochrome P450s (CYP12F2, CYP307A1, CYP4C26, CYP4C27, CYP6Z3, CYP9J3, CYP9K1, and CYP9L3), 3 were glutathione S- transferases (GSTD7, GSTE2, and GSTMS3) and 3 were undefined salivary gland protein genes Fig. 5-C, Table 2).

Table 2

Significantly differentially expressed genes of interest in R–S and C–S (FDR < 0.01 and FC > 2)
Gene categories	Gene type	Gene ID	Gene Symbol	FC (DA_KIS)	FC (BD_KIS)	FC (DD_KIS)	FC (BP_KIS)	FC (DP_KIS)	FC (BU_KIS)	FC (DU_KIS)
Detoxification genes	Carboxylesterase genes (COE)	AGAP010390	COE12O	NS	4,37	NS	NS	NS	4,21	NS
		AGAP003568	COE22933	2,77	NS	NS	NS	2,17	NS	2,17
		AGAP005372	COEBE3C	NS	2,09	NS	NS	NS	2,34	NS
		AGAP005837	COEJHE5E	NS	2,27	2,47	2,27	3,34	2,84	3,34
		AGAP013509		0,38	NS	NS	NS	NS	NS	0,37
	Cytochrome P450s genes (CYPs)	AGAP002429	CYP314A1	NS	NS	NS	NS	0,19	NS	0,19
		AGAP000192	CYP4C26	0,09	0,26	0,18	0,12	0,09	0,30	0,09
		AGAP009246	CYP4C27	2,92	3,52	2,39	2,08	2,12	3,73	2,12
		AGAP009240	CYP4C35	NS	3,09	NS	NS	NS	2,58	NS
		AGAP009241	CYP4C36	0,37	NS	NS	NS	0,40	NS	0,40
		AGAP008358	CYP4H17	2,45	NS	NS	NS	NS	NS	3,09
		AGAP002419	CYP4D22	NS	2,14	NS	2,22	3,96	2,00	3,96
		AGAP001864	CYP4H15	4,84	5,54	4,79	3,35	NS	5,35	5,39
		AGAP013128	CYP6AA2	NS	NS	2,37	2,55	2,68	2,48	2,68
		AGAP010961	CYP6AK1	0,28	0,38	0,29	NS	NS	0,45	0,41
		AGAP008020	CYP12F2	38,07	15,11	33,64	31,97	29,82	27,46	29,82
		AGAP012957	CYP4D17	3,76	NS	NS	NS	NS	NS	2,03
		AGAP005992	CYP302A1	0,27	NS	0,25	NS	0,38	NS	0,38
		AGAP002869	CYP6P2	NS	2,64	NS	NS	NS	2,07	NS
		AGAP010077	CYP303A1	NS	NS	NS	NS	14,07	NS	14,07
		AGAP008212	CYP6M2	NS	NS	NS	6,44	5,47	5,81	5,47
		AGAP001039	CYP307A1	0,35	0,15	0,45	0,20	0,27	0,17	0,27
		AGAP002865	CYP6P3	2,08	2,66	NS	4,58	4,79	3,10	4,79
		AGAP000818	CYP9K1	2,34	3,85	6,26	6,06	10,41	7,22	10,41
		AGAP002867	CYP6P4	NS	6,77	NS	4,59	6,55	6,07	6,55
		AGAP009374	CYP9M1	NS	3,14	NS	3,92	NS	2,46	NS
		AGAP012295	CYP9L1	0,43	NS	NS	NS	NS	NS	0,49
		AGAP012293	CYP9L3	0,37	0,43	0,21	0,40	0,29	0,45	0,29
		AGAP012296	CYP9J5	0,40	NS	NS	NS	0,48	NS	0,48
		AGAP012291	CYP9J3	0,35	0,41	0,47	0,30	0,26	0,39	0,26
		AGAP008217	CYP6Z3	23,07	NS	24,26	36,86	NS	26,13	23,99
		AGAP008218	CYP6Z2	3,09	NS	2,90	4,21	5,43	2,82	5,43
	Glutathione S-transferase genes (GSTs)	AGAP009193	GSTE4	NS	3,29	NS	2,13	NS	2,81	NS
		AGAP004164	GSTD1-4	NS	0,45	NS	0,47	NS	0,50	NS
		AGAP009194	GSTE2	2,69	6,20	3,24	3,95	3,58	4,90	3,58
		AGAP009946	GSTMS3	2,53	2,94	3,27	4,13	4,71	3,37	4,71
		AGAP010404	GSTS1	NS	2,85	NS	3,26	NS	2,56	NS
		AGAP004163	GSTD7	3,18	3,56	2,80	5,62	3,95	4,54	3,95
Cuticullar proteins	CPs	AGAP009790	CPAP3-B	2,71	5,92	3,57	4,39	3,63	4,74	3,63
		AGAP000047	CPR130	4,88	19,51	3,98	13,03	6,58	15,36	6,58
		AGAP007980	CPCFC1	NS	20,16	NS	11,36	NS	10,78	NS
		AGAP008444	CPLCG1	0,25	2,06	NS	NS	0,49	2,08	0,49
		AGAP008447	CPLCG4	NS	32,11	NS	25,61	NS	31,88	NS
		AGAP009405	CPAP3-E	2,37	2,54	3,03	3,13	3,84	3,51	3,84
		AGAP003334	CPLCX2	0,29	NS	0,19	NS	0,36	NS	0,36
		AGAP009759	CPLCP12	NS	9,51	NS	7,10	NS	8,66	NS
		AGAP000987	CPAP3-A1b	2,66	NS	3,47	NS	NS	NS	2,20
		AGAP001664	CPR1	NS	NS	NS	NS	9,26	NS	9,26
		AGAP002994	CPR10	2,53	2,44	2,03	2,38	2,91	2,45	2,91
		AGAP000820	CPR125	NS	3,83	NS	NS	NS	2,44	NS
		AGAP008449	CPLCG5	NS	41,81	NS	32,75	NS	30,06	NS
		AGAP009758	CPLCP11	NS	3,91	NS	NS	NS	3,66	NS
		AGAP000344	CPR127	2,99	4,68	NS	4,15	2,53	4,21	2,53
		AGAP003308	CPAP3-C	NS	3,95	NS	NS	NS	3,60	NS
		AGAP006868	CPR140	NS	7,06	NS	2,86	NS	4,55	NS
		AGAP005456	CPR15	3,62	3,67	2,97	4,03	3,85	4,60	3,85
		AGAP005459	CPR16	5,53	3,04	2,57	4,57	4,28	3,56	4,28
		AGAP006829	CPR59	NS	6,40	NS	NS	NS	5,64	NS
		AGAP006283	CPR70	NS	NS	NS	0,41	NS	0,39	NS
		AGAP009871	CPR75	NS	2,79	2,45	4,72	2,73	4,63	2,73
		AGAP009874	CPR76	9,99	5,24	5,27	7,87	7,67	8,71	7,67
		AGAP011937		NS	5,15	NS	2,37	NS	4,40	NS
		AGAP006009	CPR30	4,07	5,72	4,13	10,52	12,52	9,04	12,52
		AGAP006369	CPR144	NS	NS	NS	4,53	NS	4,83	NS
		AGAP011936		NS	2,32	NS	NS	NS	2,24	NS
		AGAP008512		NS	2,24	NS	2,06	NS	2,04	NS
		AGAP001748		NS	3,62	NS	3,03	NS	3,57	NS
		AGAP000352	TWDL1	3,52	NS	3,43	NS	NS	NS	2,64
		AGAP009879	CPR81	2,86	3,35	2,66	3,29	3,53	3,51	3,53
		AGAP002726	CPR9	NS	2,82	NS	NS	NS	2,33	NS
Salivary gland proteins	SGPs	AGAP008216	SG7	NS	4,48	NS	5,53	NS	2,40	NS
		AGAP008279	D7L2	NS	2,94	NS	NS	NS	2,25	NS
		AGAP008282	D7r2	NS	4,78	NS	NS	NS	2,43	NS
		AGAP008278	D7L1	NS	4,42	NS	NS	NS	3,35	NS
		AGAP008280	D7r5	0,31	NS	NS	NS	0,46	NS	0,46
		AGAP000612	SG1	0,41	NS	0,39	NS	NS	NS	0,42
		AGAP003841	SG10	0,07	NS	NS	0,29	0,20	0,48	0,20
		AGAP000607	SG1c	NS	2,78	NS	NS	NS	2,08	NS
		AGAP003713		NS	NS	NS	NS	0,43	NS	0,43
		AGAP006507	SG3	NS	4,59	NS	NS	NS	3,25	NS
		AGAP013423	SG9	NS	NS	NS	2,51	2,62	2,07	2,62
		AGAP006899		NS	2,02	NS	NS	NS	2,08	NS
		AGAP006442		0,21	NS	0,46	NS	0,47	NS	0,47
		AGAP006506	SG2	NS	4,29	3,59	3,49	2,95	3,69	2,95
		AGAP007851		2,28	NS	2,32	4,15	4,28	3,51	4,28
		AGAP013060		7,62	3,92	4,31	4,20	5,47	4,72	5,47
		AGAP000152		NS	NS	0,27	NS	0,40	NS	0,40
		AGAP008782		NS	NS	NS	2,81	NS	2,79	NS
		AGAP001989		NS	3,79	NS	2,07	NS	2,82	NS
		AGAP005822		NS	2,82	NS	NS	NS	2,29	NS
		AGAP011026		NS	2,76	NS	3,03	NS	2,37	NS
		AGAP012926		21,96	21,98	34,61	18,38	34,50	30,36	34,50
		AGAP006495		0,22	NS	0,22	NS	0,35	NS	0,35
		AGAP004316		5,75	10,10	9,75	5,32	NS	11,86	4,92
		AGAP006432		NS	NS	NS	NS	0,45	NS	0,45

Non-synonymous target site mutations

Target site mutations in the Ace-1 and in the kdr genes were identified through the analysis of single nucleotide polymorphisms (SNPs). The analyzed samples were composed of pools of seven mosquitoes, and a single population included three replicates (pools). The frequency at which polymorphisms appeared in a population was derived from the depth coverage at each position with the contribution of each mosquito in the pools. We found the G280S mutant allele frequency in the Ace-1 gene to be 50% for Bassila compared to 66% for Djougou (Fig. 6). For the kdr gene, we observed the L995F mutation in 33% of samples from Bassila and 100% of samples from Djougou. Other mutations observed in that gene are shown in Fig. 6.

Gene ontology annotation and enrichment analysis

Gene ontology enrichment (GOE) analysis was conducted on differentially expressed genes (up and downregulated) for all R-S comparisons. Gene ontologies are classified into three classes: biological process (BP), molecular function (MF) and cellular component (CC) (Additional file 10).

Among BP GO terms, terms related to proton transport (GO:0015986), aerobic respiration (GO:0009060) and mitochondrial electron transport (GO:0006123, GO:0006120) were shared in the BD vs KIS, BP vs KIS, DD vs KIS and DP vs KIS comparisons. Furthermore, GO terms associated with glycolytic processes (GO:0006096) were enriched in BP vs KIS and DP vs KIS; while the carbohydrate metabolic processes (GO:0005975) were enriched in BD vs KIS and DD vs KIS (Additional file 11). Considering the CC GO terms, terms associated with mitochondrial respiratory chain complex I (GO:0005747) III (GO:0005750) and IV(GO:0005751) were enriched in BD vs KIS, BP vs KIS, DD vs KIS and DP vs KIS with mitochondrial proton-transporting ATP synthase complex, and coupling factor F(o) (GO:0000276) was enriched in only the BP vs KIS and DP vs KIS comparisons (Additional file 11). Concerning terms related to MF, macromolecular complex binding related terms (GO:0044877) and sulfur cluster binding (GO:0051539) were enriched in all R-S comparisons except DA vs KIS. Hydrogen ion transmembrane transporter activity (GO:0015078) was enriched in both BP and DP vs KIS, highlighting its strong correlation with exposure to pirimiphos-methyl (Additional file 11).

RNA-Seq data validation using quantitative PCR

The expression patterns of four genes (SG7, CYP9K1, CYP6P3 and COEJHE5E) were validated in relation to two housekeeping genes (40S ribosomal protein S7; RPS7 and Actin5c) (Fig. 7). Most of the qPCR results supported the directionality of the expression level changes observed after RNA sequencing. Nevertheless, for the COEJHE5E gene, while downexpression was observed in DP vs KIS after qPCR, an overexpression was observed in the RNA-seq dataset.

In this study, we analyzed the gene expression profiles of two An. gambiae populations following exposure to key insecticides used for malaria vector control. In addition to detecting differential expression of genes encoding detoxification enzymes, genes coding for salivary gland and cuticular proteins were also overexpressed in mosquitoes from both study sites.

The two mosquito populations displayed variable levels of resistance to alphacypermethrin, deltamethrin and pirimiphos-methyl. A higher level of resistance to both alphacypermethrin and deltamethrin was observed in Djougou than in Bassila despite both being located in Donga district in northern Benin. Such heterogeneity in the levels of insecticide resistance despite relative geographic proximity has also been described elsewhere [13, 14]. The high pyrethroid resistance detected in the Djougou population could be related to the intensive use of pyrethroids in vector control by the NMCP which has intensified ITN distribution every three years since 2008 in Djougou [15]. In addition, there is intensive use of carbamates, organophosphates and pyrethroid pesticides for agricultural activity in Donga district [15, 16]⁠. In contrast, in Bassila, apart from being an agricultural site, no NMCP interventions have been distributed to date. The lower alphacypermethrin and deltamethrin resistance in mosquitoes from this locality might thus be related to reduced selective pressure. For this reason, the increased resistance to pirimiphos-methyl is surprising as it was higher than that of Djougou.

Cytochrome P450s are usually present in abundance in insects [17] and their primary function is to metabolize pheromones and xenobiotics [17–20]. In the Djougou mosquitoes that survived alphacypermethrin exposure, two members of this family, CYP4H17 and CYP4D17, were overexpressed. Both genes have been found to be overexpressed in mosquitoes resistant to permethrin [21, 22] with an overexpression of CYP4D17 in mosquitoes exposed to 0.05% deltamethrin [23]. Their over-expression only in mosquitoes resistant to alphacypermethrin here is quite interesting given that they were not overexpressed in deltamethrin resistant from the same population. Further, CYP4G16 and CYP12F3 were overexpressed in Djougou mosquitoes that survived both alphacypermethrin and deltamethrin exposure. Overexpression of CYP4G16 has been demonstrated to be important in cuticular hydrocarbon content enrichment and thus reducing the uptake of pyrethroid insecticides [23, 24]. On the other hand, over-expression of CYP12F3 has not yet been associated with insecticide resistance. Functional validation of these four genes associated with the high intensity of resistance recorded to these class two pyrethroids would confirm the role played in conferring resistance.

Several CYP450s were found to be commonly shared among all R-S comparisons including CYP9K1, CYP6Z3 and CYP12F2. All three were shown to be consistently overexpressed in both pyrethroid and pirimiphos-methyl resistant mosquitoes. CYP9K1 has previously been shown to be involved in the metabolism of deltamethrin and pyriproxyfen in An. gambiae [23, 25] and therefore its overexpression in the pyrethroid-resistant groups is not surprising⁠. Additionally, a similar study in An. arabiensis showed its association with resistance to both organophosphates and pyrethroids [12],⁠ consistent with our findings here of its association with resistance to both chemical classes. CYP6Z3 has previously been associated with pyrethroid metabolism [24] but here we observed its overexpression in survivors exposed to both pyrethroids and pirimiphos-methyl indicating a potential role in cross-resistance as has been reported before [26–29]. CYP12F2 has been previously associated with permethrin resistance in An. arabiensis [30] [31]⁠ and was associated with both pyrethroid and organophosphate resistance in our study, suggesting a role in the detoxification of both insecticide classes. Moreover, CYP6Z2, CYP9M1 and CYP4H15 were found to be overexpressed in 4 out of 5 R-S comparison groups, with the exception of BD vs KIS. The Bassila population was not as intensely resistant to pyrethroids as Djougou, which suggests that these three genes could potentially contribute to intensified pyrethroid resistance as well as pirimiphos-methyl resistance.

Other metabolic genes including three glutathione-S-transferase genes (GSTD7, GSTE2, GSTE4) and two carboxylesterase genes (COEJHE5E, COE22933) were overexpressed in all R-S comparisons. This is consistent with the understanding that GSTs contribute to resistance against multiple insecticide classes [32]⁠. They participate in the detoxification of xenobiotics and metabolize secondary products from other metabolic activities (CYP450s, COEs) [32, 33]⁠. COEJHE5E has been recently found to be overexpressed in deltamethrin and permethrin resistant An. coluzzii population [23], while genomic signals of resistance to deltamethrin were found around the COE22933 locus in An. gambiae and An. coluzzii [34]. A functional validation of these two genes might enhance understanding of their role in pyrethroid resistance.

In addition to metabolic genes, salivary gland and cuticular proteins were also differentially expressed. In our study, a wide variety of genes encoding cuticular proteins (CPs) were found to be overexpressed in R-S comparisons. Among them, CPR10, CPAP3-E, CPAP3-B, CPRd30, CPR130, CPR15, CPR16, CPR76, CPR75, and CPR81 were overexpressed in all R-S comparisons. Their mechanism of function is not well understood but alterations to the mosquito cuticle by thickening or hardening could inhibit the penetration of insecticides and other toxicants [35]. Additionally, recent studies highlighted the massive production of chitin to enhance metabolic gene functions, such as CYPs, which could lead to cuticle hardening in Aedes aegypti [36]. If this is the case, this mechanism would lead to cross-resistance of insecticides regardless of class. CPAP3s belong to a group of proteins that maintain the structural integrity of the cuticle [37]⁠. CPAP3s with an obstructor-E function have been found to play a role in multiple insecticide resistance in An. gambiae [38]. In this study, CPAP3-B was overexpressed in all R-S comparisons making this gene a potential candidate marker for cross-resistance. Other cuticular proteins that belong to the CPR family have been associated with multiple-resistance phenotypes [35]⁠. In this study, CPR10, CPR30, CPR130, CPR15, CPR16, CPR75, CPR81 and CPR76 were overexpressed, although their roles in insecticide resistance have not yet been investigated. Additional cuticular genes were observed to be more specifically overexpressed in response to exposure to certain classes of insecticides; CPR9, CPR59, CPLCP3, TWDL1, CPAP3-A1b and CPAP3-A1c in response to pyrethroids and CPR144 in response to pirimiphos-methyl exposure. The involvement of CPLCP3 in cuticle barrier formation has been previously described and its overexpression has been associated with mosquito survival following deltamethrin exposure [39]. The overexpression of CPAP3-A1b and CPAP3-A1c has also been associated with resistance to permethrin, deltamethrin or lambdacyhalothrin [40–42].

Recently, salivary gland proteins have been reported to play a potential role in insecticide resistance [12]. In this study, three salivary gland genes stood out among the overexpressed genes in resistant mosquitoes. A threonine serine-rich mucin (SG9: AGAP013423, FC ranging from 2.07 to 2.62) and a salivary cysteine-rich peptide (AGAP011460, FC ranging from 2.03 to 3.38) were overexpressed in mosquitoes resistant to pirimiphos-methyl while salivary secreted mucin 3 (AGAP009473, FC ranging from 2.04 to 2.21) was overexpressed in Djougou mosquitoes resistant to alphacypermethrin and deltamethrin. The salivary protein SG9 is a mucin protein. The role of mucin protein in mosquitoes remains uncharacterized to date, unlike in vertebrates (including humans). They are found in the peritrophic membrane and have been reported to only be induced after blood ingestion in female mosquitoes [43]⁠. Here, we found a slight overexpression (FC ranging from 2.51 to 2.62) of this gene in mosquitoes resistant to pirimiphos-methyl. It has been demonstrated that salivary gland mucin-like protein in Drosophila could perform an immune defense reaction [44],⁠ and for this reason, we hypothesize that SG9 could perform a similar function in protecting and defending the respiratory wall against the penetration of organophosphate molecules. In humans, for example, exposure to low-level OPs triggers increased mucin secretion in asthmatic patients [45]⁠. This might also be the case in mosquitoes, highlighting the need to further investigate mucin-like protein functions.

On the other hand, AGAP011460 and AGAP009473 are two uncharacterized genes that were both highly overexpressed in mosquitoes resistant to pirimiphos-methyl and both pyrethroids. Further characterization and functional annotation are needed to understand their involvement in insecticide resistance.

In addition to differential gene expression, target site mutations are also important genetic markers of phenotypic insecticide resistance. Mutations such as L995F or L995S in the kdr gene [4] and G280S mutations in the Ace-1 gene [5] have been reported to confer resistance to insecticides. The data on target site mutations presented here are relatively low resolution, as they were extrapolated from RNA-Seq data on pooled mosquitoes. Nevertheless, these data provide some insight into the presence of these mutations in the two populations. The mutation L995F was found in Djougou with an allele frequency of one (100%) and in Bassila with an allele frequency of 0.33 (33%). This is likely the result of evolution due to the selective pressure caused by the intensive use of pyrethroid insecticides in Djougou [15, 16]⁠. A frequency of 0.66 (66%) of the G280S mutation which can confer resistance to organophosphates was detected in Djougou. although this population was less resistant to pirimiphos-methyl than Bassila where the mutation frequency was 0.5 (50%). Previous research has shown that the G280S mutation may not be a strong predictor of resistance to all organophosphates [46], which supports this study’s findings. The overexpression of other genes described above could be important contributors to the resistant phenotypes observed.

The analysis of RNA-Seq data allowed us to describe differentially expressed genes and target site mutations that are associated with pyrethroid and organophosphate resistance in An. gambiae from two locations in Benin: Bassila and Djougou. Multiple genes, including members of the cytochrome P450 family (CYP4H17, CYP4D17 and CYP12F2), salivary gland proteins (SG9, AGAP011460 and AGAP009473) and cuticular proteins (CPR30, CPR130, CPR15, CPR16, CPR76, CPAP3-A1b and CPAP3-A1c) were found overexpressed in resistant mosquitoes after exposure to alphacypermethrin, deltamethrin, and pirimiphos-methyl. The DEGs described here are potential molecular markers of insecticide resistance that could be incorporated into Benin’s NMCP insecticide resistance surveillance and management strategy once validated.

Study sites and Samples

Anopheles gambiae larvae were collected in August 2019 from Djougou (9° 42' 29.1312'' N and 1° 39' 58.8672'' E) and in October 2019 from Bassila (9° 0' 23.0148'' N and 1° 39' 50.1264'' E) (Fig. 8) from-puddles, swamp areas and drains near irrigated croplands using the dipping method. Mosquito larvae were brought back to the insectary of the ‘Environnement, Gestion des Données et Formation Universitaire (EGDFU)’ unit of the Tropical Infectious Diseases Research Center (TIDRC) of Benin Republic and reared to the adult stage under insectary conditions (insecticide-free environment, 27 ± 2°C ambient temperature, 70 ± 8% relative humidity and 12 h:12 h light:dark photoperiod). Larvae were fed ad libitum with TetraMin Baby® fish food and emerged adults were fed on 10% honey solution.

Anopheles gambiae from the insecticide susceptible reference Kisumu laboratory strain were reared in the insectary at the Centers for Disease Control and Prevention (CDC), Atlanta, Georgia, USA. Mosquitoes were maintained at a constant 27 ± 2°C and 70 ± 10% humidity on a 14 h:10 h hour light:dark cycle and adults were provided 10% sucrose ad libitum.

Insecticide resistance intensity bioassays

Insecticide resistance intensity assays using modified CDC bottle bioassays were conducted on the reared adult mosquitoes using 1X, 2X, and 5X the diagnostic doses of deltamethrin and alpha-cypermethrin and 0.5X, 1X and 1.5X the diagnostic doses of pirimiphos-methyl. Stock solutions (10X) of the diagnostic doses (1X for deltamethrin: 12.5 µg/bottle, alphacypermethrin: 12.5 µg/bottle and pirimiphos-methyl: 20 µg/bottle) were prepared by diluting technical grade insecticide in 50 mL of absolute ethanol. Lower doses were obtained by serial dilution.

For each concentration, four bottles were coated with 1 ml of insecticide solution and one control bottle was coated with absolute ethanol. Alphacypermethrin and deltamethrin coated bottles were covered to keep them protected from light and to allow evaporation of the solvent overnight while the bottles coated with pirimiphos-methyl were left covered and protected from light for 8 hours.

Approximately 10 to 25 mosquitoes aged 4–5 days and fed with 10% honey solution were released in each bottle using an aspirator. The exposure time was set to 30 min, after which the mosquitoes were removed from the bottles and sorted into “alive” and “knocked-down” groups. Mosquitoes alive after 30 minutes of exposure to insecticide were considered resistant [47, 48]. Mosquitoes from the control bottle were labeled as unexposed.

Mosquito species identification

The legs of the resistant and unexposed mosquitoes were removed, and the rest of the body was immediately stored in RNA later. DNA extraction from the legs was performed following the protocol of Myriam and Cecile (2003) using the cetyltrimethyl ammonium bromide (CTAB) technique. Mosquito species identif ication was carried out to distinguish members of the An. gambiae s.l species complex using species-specific PCR primers for An. gambiae s.s, and An. arabiensis. The PCR mix had a total volume of 20 µl composed of 5 µl of DNA, 8.4 µl of H₂O, 2.5 µl of reaction buffer 1x, 0.1 µl of Taq DNA polymerase, 1 µl of MgCl₂, and 1 µl of each primer (Universal forward primer: 5'-GTGTGCCCCTTCCTCGATGT-3'; An. gambiae s.s reverse primer:5'-CTGGTTTGGTCGGCACGTTT-3' and An. arabiensis reverse primer: 5'-AAGTGTCCTTCTCCATCCTA). The PCR cycling conditions were: 94°C for 5 min, followed by 25 amplification cycles (94°C for 30 s, 72°C for 30 s) and a final elongation step at 72°C for 5 min. The PCR products were visualized using 1.5% agarose gels stained with 5 µl of BET. The DNA bands on the gel (390 bp for An. gambiae s.s. and 315 bp for An. arabiensis) were compared to a 1 kb reference ladder.

RNA extraction, RNA-Seq library preparation and sequencing

Four-to-five-day-old adult nonblood-fed female mosquitoes from the susceptible Kisumu strain were killed by freezing and stored at -80°C until RNA extraction. Mosquitoes from the field populations unexposed to insecticides and those that survived to the highest dose of each insecticide were stored in RNAlater and shipped to the Entomology Branch laboratory at the CDC, Atlanta, USA, for RNA extraction, library preparation and sequencing.

RNA extraction was conducted using the Arcturus PicoPure RNA isolation kit (Life Technologies, USA) according to the manufacturer’s instructions from three biological replicates with pools of 6 mosquitoes each from the following groups: mosquitoes from Djougou phenotyped as resistant to pirimiphos-methyl (DP), alphacypermethrin (DA) or deltamethrin (DD); mosquitoes from Djougou unexposed to insecticides (DU); mosquitoes from Bassila phenotyped as resistant to pirimiphos-methyl (BP) or deltamethrin (BD); mosquitoes from Bassila unexposed to insecticides (BU) and mosquitoes from the susceptible Kisumu strain (KIS). The Agilent 4200 TapeStation was used to measure the RNA concentration and integrity. Ribosomal RNA was depleted using the Ribo-Zero™ Magnetic Core Kit and Ribo-Zero™ rRNA Removal kit (Illumina, USA). Library preparation was carried out using the ScriptSeq v2 RNA-Seq Library Preparation Kit (Epicenter, Illumina) according to the manufacturer’s instructions. Libraries were purified using Agencourt AMPure XP beads (Beckman Coulter, USA). The quantity and size distributions of the libraries were assessed using the Agilent DNA ScreenTape assay. Equimolar amounts of each library were pooled and sequenced (2x125 bp paired-end) on an Illumina HiSeq 2500 sequencer, using v2 chemistry. Sequencing was performed at the Biotechnology Core Facility at CDC, Atlanta, USA.

RNA-Seq data quality control and mapping

Quality control was performed on the raw dmultiplexed paired end sequencing reads obtained from the sequencing center using FastqC [49]. Reads from lane 1 and lane 2 were concatenated and trimmed to remove polyG tails, polyX at the 3' ends, bases that did not meet the minimum quality score of 20 and paired reads where one or both were shorter than 50 bp using Fastp [50]. The trimmed reads (R1/R2) were aligned to the An. gambiae PEST reference genome (GenBank assembly identifier = GCA_000005575.2), directly downloaded from Vectorbase (release48) using 'subjunc', which is part of the subread aligner package, version 2.0.1 [51] with default settings. The alignments were filtered to remove low quality mapping reads (< 10) and sorted using SAMtools (v.1.10) [52].

FeatureCounts from the subread package version 2.0.1 [51, 53] was used to count tags that overlapped the coding sequence (CDS) features by at least 1 bp in the sense orientation of the gene set AgambiaePest (structural annotation version = AgamP4.13). Tag counts represent the number of sequence reads that originated from a particular gene. The higher the number of counts, the more reads associated with that gene, and the assumption that there is a higher level of expression of that gene in the sample. The tag count result was used as input for differential expression analysis using edgeR, version 4.1.1 [54] .

Differential expression analysis was performed using the package EdgeR to decide whether, for a given gene, an observed difference in tag counts was significantly different from what would be expected due to natural random variation. Thus, to remove the effect of noise and genes with very low expression, only genes where at least a tag count of 50 or more was obtained across all libraries were considered.

CalcNormFactors, a function of the edgeR package that uses the TMM (trimmed mean M-values) method, was used to normalize the number of tags between samples, by finding a set of scaling factors for the library sizes that minimize log-fold changes between samples for most genes. DEGs between the different comparisons were selected after multiple testing using the decideTests function, of the limma package[55]. An absolute fold-change > 2 and FDR (false discovery Rate) ≤ 0.01 were used as statistical cutoffs to tag a gene as a DEG.

Comparisons were made between i) resistant field mosquitoes and laboratory susceptible mosquitoes (R-S): (DA vs KIS; DD vs KIS; DP vs KIS; BD vs KIS; BP vs KIS); ii) unexposed field mosquitoes and laboratory susceptible mosquitoes (C-S): (DU vs KIS; BU vs KIS) and iii) resistant field mosquitoes and unexposed field mosquitoes (R-C): (DA vs DU; DD vs DU; DP vs DU; BD vs BU; BP vs BU). R-S and C-S comparisons were made following the assumptions that constitutive resistance genes would be differentially expressed between both bioassay survivors and the unexposed mosquitoes when compared to the susceptible strain. The R-C comparison was made to detect gene expression that may be induced due to insecticide exposure. In addition, insecticide-specific comparisons were made (DD vs BD and DP vs BP) to detect DEGs associated with the same insecticide across the different mosquito populations.

Detection of non-synonymous target site mutations

Target site mutations were detected from the variant calling analysis. Sorted bam files were used to count the read coverage at each genomic position and provide genotype likelihoods using the ‘mpileup’ methods of BCFtools [52, 56]⁠. Single nucleotide variants (SNVs) were then detected using bcftool ‘call’ and the resulting vcf file was filtered with vcfutils.pl [52, 56]⁠.

In each population, the frequencies of known point mutations in the Ace-1 gene (G280S) and in the kdr gene (A1125V, A1746S, A1934V, D466H, E1597G, F1920S, I1527T, I1868T, I1940T, K1603T, L995F, L995S, M490I, N1575Y, P1874L, P1874S, T791M, V1254I, V1853I, V402L, and V1853I) were recorded and visualized.

Gene ontology enrichment analysis

Gene ontology enrichment analysis was used to classify the differentially expressed genes based on the specific biological functions using Goatools [57]⁠. This approach allows enrichment analysis to indicate which molecular functions, biological processes or cellular components were overrepresented (enriched) in a DEG list compared to an annotated list of the whole genome of An. gambiae obtained from blast2go [58] (Additional file 12)⁠. The, GO term enrichment analysis of up- and downregulated genes was carried out using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) [59, 60]. The P-values used to assess significantly enriched GO terms associated with the DEG list were calculated based on Fisher's exact test and corrected by the Benjamini-Hochberg multiple testing correction method. Thus, an FDR adjusted P-value < 0.05 was used to report significantly enriched GO terms from the list of DEGs (Additional file 10). Visualization of the enriched GO terms was performed using the package ‘ggplot2’ following the protocol of Bonnot and others [61]⁠.

RNAseq data validation using quantitative PCR

To validate the RNA-seq data, a set of six genes (SG7, CYP9K1, CYP6P3, COEJHE5E, RPS7, Actinc5) was used. RNA from three replicates of samples resistant to alphacypermethrin, deltamethrin or pirimiphos-methyl from Djougou (same batch of mosquitoes used for RNA sequencing) was used to synthesize cDNA using the HighCapacity cDNA reverse transcription kit (Applied Biosystems) with oligo-dT20 (NEB), according to the manufacturer’s instructions. The primers used are listed in Additional file 13. Standard curves of Ct values for each gene were generated using a serial dilution of cDNA, allowing assessment of PCR efficiency. qPCR amplification was carried out on an Agilent Technologies Stratagene Mx3005P using PowerUp SYBR Green Master Mix (Applied Biosystems). cDNA from each sample was used as a template in a three-step program: 50°C for 2 minutes denaturation at 95°C for 10 minutes, followed by 40 cycles of 15 seconds at 95°C, 1 minute at 60°C and a final step of 15 seconds at 95°C, 1 minute at 60°C, and 15 seconds at 95°C. The relative expression level and fold change (FC) of each target gene from resistant field samples relative to the susceptible lab strain were calculated using the 2 − ΔΔCT method [62] incorporated in Python script (https://github.com/dany-gaga/from_Ct_to_logFoldChange). Housekeeping genes encoding ribosomal protein S17 (RPS17; AGAP010592) and Actin5C (AGAP000651) were used for normalization.

Insecticide resistance

ITNs

Insecticide-treated nets

IRS

indoor residual spraying

DGE

differential gene expression

Fold-change

CPs

Cuticular proteins

CYPs

Cytochrome P450 genes

GSTs

Glutathione-S-transferases

COE

Carboxylesterases

SNPs

Single nucleotide polymorphisms

GOE

Gene ontology

Biological process

Molecular function

Cellular component

NMCP

National Malaria Control Program

LLINs

Long-lasting insecticidal nets

SGPs

Salivary gland proteins

RPS17

Ribosomal protein S1

CDS

coding sequence

Gene ontology

DAVID

Database for Annotation, Visualization, and Integrated Discovery

CTAB

Cetyl trimethyl ammonium bromide

Ethics approval and consent to participate

Ethics approval was granted by the Kenya Medical Research Institute (KEMRI), scientific Ethics Review Unit steering committee (SERU 3275). Verbal consent was obtained from the head of each household prior to mosquito collections.

Consent for publication

All authors have approved the publication of this manuscript.

Availability of data and materials

The dataset supporting the conclusions of this article is available at Sequence Read Archive (SRA) under the accession number PRJNA982704.

Competing interests

The author declares that they have no competing interests.

Funding

This work was supported by the Bill and Melinda Gates Foundation, Grand Challenges Grant No. INV024969 and Investment Grant OPP1210769.

Authors' contributions

EO, LSD, LI, ND and AL: project design and funding acquisition; HS, DD and LI: pipeline design, formal analysis, writing of the original draft and editing; DO, SO, NM, EO, LSD, ND and AL: manuscript review and editing.

Acknowledgments

We would like to thank all of the ‘Unité Environnement, Gestion des Données et Formation Universitaire’ (UEGDFU), Pierre Marie Sovegnon for his assistance during larval collection and ideas brainstorming, Adande Bernard Medjigbodo for his help during the bioassays and Claudia Corredor for her help during RNA extractions. Furthermore, we express our gratitude to Yassimine Djibril, Flavia Hounwanou, Horace Agossadou, Nicole Dzuris for facilitating the administrative tasks related to this study. The views expressed in this manuscript are those of the authors and do not necessarily reflect the official policy or position of the Centers for Disease Control and Prevention.

Raghavendra K, Barik TK, Reddy BPN, Sharma P, Dash AP. Malaria vector control: From past to future. Parasitol Res. 2011;108:757–79.
WHO. World Malaria Report 2019. Geneva. 2019.
Liu N. Insecticide Resistance in Mosquitoes: Impact, Mechanisms, and Research Directions. http://dx.doi.org/101146/annurev-ento-010814-020828. 2015;60:537–59.
Chandre F, Darrier F, Manga L, Akogbeto M, Faye O, Mouchet J, et al. Status of pyrethroid resistance in Anopheles gambiae sensu lato. Bull World Health Organ. 1999;77:230–4.
Weill M, Luffalla G, Mogensen K, Chandre F, Berthomieu A, Berticat C, et al. Insecticide resistance in mosquito vectors. Nature. 2003;423:136–7.
Riveron JM, Tchouakui M, Mugenzi L, D.Menze B, Chiang M-C, Wondji CS. Insecticide Resistance in Malaria Vectors: An Update at a Global Scale. Towar Malar Elimin - A Leap Forw. 2018. https://doi.org/10.5772/INTECHOPEN.78375.
Tchouakui M, Miranda JR, Mugenzi LMJ, Djonabaye D, Wondji MJ, Tchoupo M, et al. Cytochrome P450 metabolic resistance (CYP6P9a) to pyrethroids imposes a fitness cost in the major African malaria vector Anopheles funestus. Heredity (Edinb). 2020;124:621–32.
Black WC, Snell TK, Saavedra-Rodriguez K, Kading RC, Campbell CL. From Global to Local—New Insights into Features of Pyrethroid Detoxification in Vector Mosquitoes. Insects 2021, Vol 12, Page 276. 2021;12:276.
Faucon F, Gaude T, Dusfour I, Navratil V, Corbel V, Juntarajumnong W et al. In the hunt for genomic markers of metabolic resistance to pyrethroids in the mosquito Aedes aegypti: An integrated next-generation sequencing approach. PLoS Negl Trop Dis. 2017;11.
Neafsey DE, Taylor AR, MacInnis BL. Advances and opportunities in malaria population genomics. Nat Rev Genet 2021 228. 2021;22:502–17.
Mackenzie-Impoinvil L, Weedall GD, Lol JC, Pinto J, Vizcaino L, Dzuris N et al. Contrasting patterns of gene expression indicate differing pyrethroid resistance mechanisms across the range of the New World malaria vector Anopheles albimanus. PLoS ONE. 2019;14.
Messenger LA, Impoinvil LM, Derilus D, Yewhalaw D, Irish S, Lenhart A. A whole transcriptomic approach provides novel insights into the molecular basis of organophosphate and pyrethroid resistance in Anopheles arabiensis from Ethiopia. Insect Biochem Mol Biol. 2021;139:103655.
Hancock PA, Hendriks CJM, Tangena JA, Gibson H, Hemingway J, Coleman M, et al. Mapping trends in insecticide resistance phenotypes in African malaria vectors. PLOS Biol. 2020;18:e3000633.
Ochomo E, Bayoh NM, Kamau L, Atieli F, Vulule J, Ouma C, et al. Pyrethroid susceptibility of malaria vectors in four Districts of western Kenya. Parasites and Vectors. 2014;7:1–9.
Akogbéto MC, Dagnon F, Aïkpon R, Ossé R, Salako AS, Ahogni I, et al. Lessons learned, challenges and outlooks for decision-making after a decade of experience monitoring the impact of indoor residual spraying in Benin, West Africa. Malar J. 2020;19:1–16.
Assogba BS, Pasteur N, Makoundou P, Unal S, Baba-Moussa L, Labbé P, et al. Dynamic of resistance alleles of two major insecticide targets in Anopheles gambiae (s.l.) populations from Benin, West Africa. Parasites and Vectors. 2020;13:1–16.
Ye M, Nayak B, Xiong L, Xie C, Dong Y, You M et al. The Role of Insect Cytochrome P450s in Mediating Insecticide Resistance. Agric. 2022;12.
Berenbaum MR. Postgenomic chemical ecology: From genetic code to ecological interactions. J Chem Ecol. 2002;28:873–96.
Claudianos C, Ranson H, Johnson RM, Biswas S, Schuler MA, Berenbaum MR, et al. A deficit of detoxification enzymes: pesticide sensitivity and environmental response in the honeybee. Insect Mol Biol. 2006;15:615–36.
Zhang H, Zhao M, Liu Y, Zhou Z, Guo J. Identification of cytochrome P450 monooxygenase genes and their expression in response to high temperature in the alligatorweed flea beetle Agasicles hygrophila (Coleoptera: Chrysomelidae). Sci Reports 2018 81. 2018;8:1–13.
Sandeu MM, Mulamba C, Weedall GD, Wondji CS. A differential expression of pyrethroid resistance genes in the malaria vector Anopheles funestus across Uganda is associated with patterns of gene flow. PLoS ONE. 2020;15:e0240743.
Ibrahim SS, Muhammad A, Hearn J, Weedall GD, Nagi SC, Mukhtar MM, et al. Molecular drivers of insecticide resistance in the Sahelo-Sudanian populations of a major malaria vector Anopheles coluzzii. BMC Biol. 2023;21:1–23.
Ibrahim SS, Muhammad A, Hearn J, Weedall GD, Nagi SC, Mukhtar MM, et al. Molecular drivers of insecticide resistance in the Sahelo-Sudanian populations of a major malaria vector Anopheles coluzzii. BMC Biol. 2023;21:0–3.
Balabanidou V, Kampouraki A, Maclean M, Blomquist GJ, Tittiger C, Juárez MP, et al. Cytochrome P450 associated with insecticide resistance catalyzes cuticular hydrocarbon production in Anopheles gambiae. Proc Natl Acad Sci U S A. 2016;113:9268–73.
Vontas J, Grigoraki L, Morgan J, Tsakireli D, Fuseini G, Segura L, et al. Rapid selection of a pyrethroid metabolic enzyme CYP9K1 by operational malaria control activities. Proc Natl Acad Sci U S A. 2018;115:4619–24.
David JP, Ismail HM, Chandor-Proust A, Paine MJI. Role of cytochrome P450s in insecticide resistance: impact on the control of mosquito-borne diseases and use of insecticides on Earth. Philos Trans R Soc B Biol Sci. 2013;368.
Adebiyi MO, Olugbara OO. Homology Modeling of CYP6Z3 Protein of Anopheles Mosquito. Adv Sci Technol Eng Syst J. 2021;6:580–5.
Haruna AS, Vontas J, Eyo JE. First Incidence of CYP9K1, CYP6P4 and CYP6Z1 in Anopheles gambiae ss from Nigeria. 2020. https://doi.org/10.21203/rs.3.rs-19996/v1.
Irving H, Riveron JM, Ibrahim SS, Lobo NF, Wondji CS. Positional cloning of rp2 QTL associates the P450 genes CYP6Z1,CYP6Z3 and CYP6M7 with pyrethroid resistance in the malaria vectorAnopheles funestus. Heredity (Edinb). 2012;109:383.
Nardini L, Christian RN, Coetzer N, Koekemoer LL. DDT and pyrethroid resistance in Anopheles arabiensis from South Africa. Parasit Vectors. 2013;6:229.
Ingham VA, Wagstaff S, Ranson H. Transcriptomic meta-signatures identified in Anopheles gambiae populations reveal previously undetected insecticide resistance mechanisms. Nat Commun 2018 91. 2018;9:1–11.
Enayati AA, Ranson H, Hemingway J. Insect glutathione transferases and insecticide resistance. Insect Mol Biol. 2005;14:3–8.
Pavlidi N, Vontas J, Van Leeuwen T. The role of glutathione S-transferases (GSTs) in insecticide resistance in crop pests and disease vectors. Curr Opin Insect Sci. 2018;27:97–102.
Lucas ER, Nagi SC, Egyir-Yawson A, Essandoh J, Dadzie S, Chabi J, et al. Genome-wide association studies reveal novel loci associated with pyrethroid and organophosphate resistance in Anopheles gambiae and Anopheles coluzzii. Nat Commun. 2023;14:4946.
Xu Y, Xu J, Zhou Y, Li X, Meng Y, Ma L, et al. CPR63 promotes pyrethroid resistance by increasing cuticle thickness in Culex pipiens pallens. Parasites and Vectors. 2022;15:1–7.
Jacobs E, Chrissian C, Rankin-Turner S, Wear M, Camacho E, Broderick NA, et al. Cuticular profiling of insecticide resistant Aedes aegypti. Sci Rep. 2023;13:10154.
Chen EH, Hou QL, Dou W, Wei DD, Yue Y, Yang RL, et al. Genome-wide annotation of cuticular proteins in the oriental fruit fly (Bactrocera dorsalis), changes during pupariation and expression analysis of CPAP3 protein genes in response to environmental stresses. Insect Biochem Mol Biol. 2018;97:53–70.
Yahouédo GA, Chandre F, Rossignol M, Ginibre C, Balabanidou V, Mendez NGA, et al. Contributions of cuticle permeability and enzyme detoxification to pyrethroid resistance in the major malaria vector Anopheles gambiae. Sci Rep. 2017;2017 71:7:1–10.
Kefi M, Charamis J, Balabanidou V, Ioannidis P, Ranson H, Ingham VA et al. Transcriptomic analysis of resistance and short-term induction response to pyrethroids, in Anopheles coluzzii legs. BMC Genomics. 2021;22.
Matiya DJ, Philbert AB, Kidima W, Matowo JJ. Dynamics and monitoring of insecticide resistance in malaria vectors across mainland Tanzania from 1997 to 2017: A systematic review. Malar J. 2019;18:1–16.
Toé KH, N’Falé S, Dabiré RK, Ranson H, Jones CM. The recent escalation in strength of pyrethroid resistance in Anopheles coluzzi in West Africa is linked to increased expression of multiple gene families. BMC Genomics. 2015;16:1–11.
Ibrahim SS, Muhammad A, Hearn J, Weedall GD, Nagi SC, Mukhtar MM et al. Molecular drivers of insecticide resistance in the Sahelo-Sudanian populations of a major malaria vector. bioRxiv. 2022;:2022.03.21.485146.
Deng F, Wu S, Wu Y, Liu X, Wu P, Zhai Z. Identification of mucins and their expression in the vector mosquito Aedes albopictus. J Vector Ecol. 2020;45:297–305.
Korayem AM, Fabbri M, Takahashi K, Scherfer C, Lindgren M, Schmidt O, et al. A Drosophila salivary gland mucin is also expressed in immune tissues: Evidence for a function in coagulation and the entrapment of bacteria. Insect Biochem Mol Biol. 2004;34:1297–304.
Canesi L, Negri A, Barmo C, Banni M, Gallo G, Viarengo A et al. The Organophosphate Chlorpyrifos Interferes with the Responses to17β-Estradiol in the Digestive Gland of the Marine Mussel Mytilusgalloprovincialis. PLoS ONE. 2011;6.
Wipf NC, Duchemin W, Kouadio FPA, Fodjo BK, Sadia CG, Mouhamadou CS, et al. Multi-insecticide resistant malaria vectors in the field remain susceptible to malathion, despite the presence of Ace1 point mutations. PLOS Genet. 2022;18:e1009963.
Brogdon WG. (Centers for DC and P. Guideline for Evaluating Insecticide Resistance in Vectors Using the CDC Bottle Bioassay. 2010.
Vatandoost H, Abai MR, Akbari M, Raeisi A, Yousefi H, Sheikhi S, et al. Comparison of CDC Bottle Bioassay with WHO Standard Method for Assessment Susceptibility Level of Malaria Vector, Anopheles stephensi to Three Imagicides. J Arthropod Borne Dis. 2019;13:17.
LaMar D. FastQC. 2015.
Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34:i884–90.
Liao Y, GK S. W S. The Subread aligner: fast, accurate and scalable read mapping by seed-and-vote. Nucleic Acids Res. 2013;41.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25:2078–9.
Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30:923–30.
McCarthy DJ, Chen Y, Smyth GK. Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 2012;40:4288–97.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47.
Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, et al. Twelve years of SAMtools and BCFtools. Gigascience. 2021;10:1–4.
Klopfenstein DV, Zhang L, Pedersen BS, Ramírez F, Vesztrocy AW, Naldi A et al. GOATOOLS: A Python library for Gene Ontology analyses. Sci Rep. 2018;8.
Conesa A, Götz S, García-Gómez JM, Terol J, Talón M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinforma Appl NOTE. 2005;21:3674–6.
Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4:44–57.
Sherman BT, Hao M, Qiu J, Jiao X, Baseler MW, Lane HC, et al. DAVID: a web server for functional enrichment analysis and functional annotation of gene lists (2021 update). Nucleic Acids Res. 2022;50:W216–21.
Bonnot T, Gillard MB, Nagel DH. A Simple Protocol for Informative Visualization of Enriched Gene Ontology Terms. BIO-PROTOCOL. 2019. https://doi.org/10.21769/BioProtoc.3429.
Rao X, Huang X, Zhou Z, Lin X. An improvement of the 2ˆ(–delta delta CT) method for quantitative real-time polymerase chain reaction data analysis. Biostat Bioinforma Biomath. 2013;3:71.

No competing interests reported.

Additionalfile1.xlsx
Additional files File name: Additional file 1 File format: Additional_file1.xlsx Title of data: Insecticide intensity CDC bottle bioassay for An. gambiae collected from Bassila and Djougou
Additionalfile2.xlsx
File name: Additional file 2 File format: Additional_file2.xlsx Title of data: Descriptive statistics for the alphacypermethrin, deltamethrin and pyrimiphos-methyl read libraries and sequencing alignments.
Additionalfile3.pdf
File name: Additional file 3 File format: Additional_file3.pdf Title of data: Summary statistics of the read quantification.
Additionalfile4.xlsx
File name: Additional file 4 File format: Additional_file4.xlsx Title of data: Results of the differential gene expression analysis for all the comparisons.
Additionalfile5.pdf
File name: Additional file 5 File format: Additional_file5.pdf Title of data: Venn diagram showing differentially expressed genes among resistant, susceptible and unexposed mosquito populations. Description of data: Panel A represents gene differentially expressed in Bassila when mosquitoes are resistant to deltamethrin; in B when they are exposed to pirimiphos-methyl; in C when mosquitoes from Djougou are exposed to deltamethrin; in D when exposed to pirimiphos-methyl and in E when exposed to alphacypermethrin. Each Venn diagram section shows the number of differentially expressed genes meeting each set of conditions and the p values were adjusted for multiple testing based on FDR<0.01 and FC > 2.
Additionalfile6.xlsx
File name: Additional file 6 File format: Additional_file6.xlsx Title of data: Differentially expressed genes that overlap R-S and C-S comparisons in all experiments.
Additionalfile7.xlsx
File name: Additional file 7 File format: Additional_file7.xlsx Title of data: Differentially expressed genes that overlap R-S and R-C comparisons in all experiments.
Additionalfile8.xlsx
File name: Additional file 8 File format: Additional_file8.xlsx Title of data: Differentially expressed genes that overlap between R-S comparisons in pyrethroids and organophosphates
Additionalfile9.xlsx
File name: Additional file 9 File format: Additional_file9.xlsx Title of data: Differentially expressed genes involved in multiple-insecticide resistance
Additionalfile10.xlsx
File name: Additional file 10 File format: Additional_file10.xlsx Title of data: Gene ontology enrichment analysis
Additionalfile11.pdf
File name: Additional file 11 File format: Additional_file11.pdf Title of data: Gene ontology enrichments Description of data: Panel A shows enriched terms within biological process; B, the enriched terms in molecular functions and C, the cellular components.
Additionalfile12.xlsx
File name: Additional file 12 File format: Additional_file12.xlsx Title of data: All annotation results (VectorBase and Blast2GO)
Additionalfile13.xlsx
File name: Additional file 13 File format: Additional_file13.xlsx Title of data: List of genes to be validated using qPCR

Download PDF

Journal Publication

published 06 Apr, 2024

Read the published version in BMC Genomics →

Editorial decision: Revision requested
05 Nov, 2023
Reviews received at journal
23 Oct, 2023
Reviewers agreed at journal
06 Oct, 2023
Reviewers invited by journal
06 Oct, 2023
Editor assigned by journal
06 Oct, 2023
Editor invited by journal
06 Oct, 2023
Submission checks completed at journal
06 Oct, 2023
First submitted to journal
03 Oct, 2023

You are reading this latest preprint version

Transcriptomic analysis of Anopheles gambiae from Benin reveals overexpression of salivary and cuticular proteins associated with cross-resistance to pyrethroids and organophosphates

Status:

Journal Publication

Version 1

Abstract

Background

Results

Conclusion

Figures

Background

Results

RNA-Seq data quality control and mapping

Differential gene expression analysis

Differential gene expression associated with deltamethrin

Differential gene expression associated with pirimiphos-methyl

Genes associated with multiple insecticide resistance

Non-synonymous target site mutations

Gene ontology annotation and enrichment analysis

RNA-Seq data validation using quantitative PCR

Discussion

Conclusions

Methods

Study sites and Samples

Insecticide resistance intensity bioassays

Mosquito species identification

RNA extraction, RNA-Seq library preparation and sequencing

RNA-Seq data quality control and mapping

Detection of non-synonymous target site mutations

Gene ontology enrichment analysis

RNAseq data validation using quantitative PCR

Abbreviations

Declarations

Ethics approval and consent to participate

Consent for publication

Availability of data and materials

Competing interests

Funding

Authors' contributions

Acknowledgments

References

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1