Non-Coding RNAs Potentially Involved in Pyrethroid Resistance of Anopheles funestus Population in Western Kenya

Backgrounds The resurgence of Anopheles funestus, a dominant vector of human malaria in western Kenya was partly attributed to insecticide resistance. However, evidence on the molecular basis of pyrethroid resistance in western Kenya is limited. Noncoding RNAs (ncRNAs) form a vast class of RNAs that do not code for proteins and are ubiquitous in the insect genome. Here, we demonstrated that multiple ncRNAs could play a potential role in An. funestusresistance to pyrethroid in western Kenya. Materials and Methods Anopheles funestus mosquitoes were sampled by aspiration methods in Bungoma, Teso, Siaya, Port Victoria and Kombewa in western Kenya. The F1 progenies were exposed to deltamethrin (0.05%), permethrin (0.75%), DDT (4%) and pirimiphos-methyl (0.25%) following WHO test guidelines. A synergist assay using piperonyl butoxide (PBO) (4%) was conducted to determine cytochrome P450s’ role in pyrethroid resistance. RNA-seq was conducted on a combined pool of specimens that were resistant and unexposed, and the results were compared with those of the FANG susceptible strain. This approach aimed to uncover the molecular mechanisms underlying pyrethroid resistance. Results Pyrethroid resistance was observed in all the sites with an average mortality rate of 57.6%. Port Victoria had the highest level of resistance to permethrin (MR=53%) and deltamethrin (MR=11%) pyrethroids. Teso had the lowest level of resistance to permethrin (MR=70%) and deltamethrin (MR=87%). Resistance to DDT was observed only in Kombewa (MR=89%) and Port Victoria (MR=85%). A full susceptibility to P-methyl (0.25%) was observed in all the sites. PBO synergist assay revealed high susceptibility (>98%) to the pyrethroids in all the sites except for Port Victoria (MR=96%, n=100). Whole transcriptomic analysis showed that most of the gene families associated with pyrethroid resistance comprised non-coding RNAs (67%), followed by imipenemase (10%), cytochrome P450s (6%), cuticular proteins (5%), olfactory proteins (4%), glutathione S-transferases (3%), UDP-glycosyltransferases (2%), ATP-binding cassettes (2%) and carboxylesterases(1%). Conclusions This study unveils the molecular basis of insecticide resistance in An. funestus in western Kenya, highlighting for the first time the potential role of non-coding RNAs in pyrethroid resistance. Targeting non-coding RNAs for intervention development could help in insecticide resistance management.


Introduction
Vector control particularly the use of bed nets treated with pyrethroids has had an impact on entomological parameters, such as reducing infection rate in the vector population, vector abundance, and parity rate [1,2] leading to a decline in malaria morbidity and mortality in sub-Saharan Africa as a result of a decline in vectorial capacity [3][4][5].Despite these successes, malaria resurgence and outbreaks have been reported in various transmission settings in sub-Saharan Africa where ITNs were deployed [6][7][8][9].Hence, the effectiveness of the primary vector control methods with regard to insecticide resistance needs continuous monitoring and probing of the resistance mechanisms.
In contrast to other major vectors, Anopheles funestus sensu stricto (hereafter An. funestus) has received very scant attention owing to the di culty in colonizing this species under laboratory conditions.An. funestus is distributed throughout Africa similar to the distributed union of An. gambiae.After developing resistance and exhibiting behavioural adaptability, An. funestus has a higher ability to colonize a niche [10,11].It is one of the most ubiquitous and e cient malaria vectors in the world; highly susceptible to the P. falciparum parasite, highly anthropophagic and endophilic [12][13][14].
The signi cance of studying this mosquito is highlighted by its versatility in ecological adaptation and the emergence of resistance to recommended public health insecticides for vector control [10,15].
Increased resistance to pyrethroids used for bed net impregnation has led to low e cacy of conventional LLINs against An.funestus [16].Resistance monitoring focuses on transmission foci, hotspots of localized outbreaks, or after spikes in disease cases in pre-elimination and elimination settings [17].For effective insecticide resistance management, it is essential to genetically characterize insecticide resistance pro les and mechanisms in the vector populations.Metabolic resistance poses the biggest threat to the control of malaria vectors [18].Cytochrome P450s, Glutathione S-transferases (GSTs) and carboxylesterases (COEs) are well-established enzyme families in malaria vectors known to confer resistance to pyrethroids [19,20].These detoxi cation genes are pivotal in the molecular mechanism of insecticide resistance.
Non-coding RNAs (ncRNAs) form a vast class of RNAs that do not code for protein.Examples of ncRNAs include transfer RNA (tRNA), ribosomal RNA (rRNA), small nuclear RNA (snRNA), small nucleolar RNA (snoRNA), microRNA (miRNA), PIWIinteracting RNA (piRNA), endogenous small interfering RNA (siRNA), circular RNA (circRNA), long non-coding RNA (lncRNA), protein functional effector small ncRNA (pfeRNA), and other ncRNAs whose functions remain unknown [21,22].They can control the expression of genes at the chromosomal, transcriptional, post-transcriptional, and translational levels and play a role in the entire developmental process.ncRNAs have been demonstrated in studies on arthropods to be essential for several physiological and developmental processes, including molting, reproduction, immunity, wing development, and insecticide resistance [23].ncRNAs can modify signalling pathways involved in these biological processes by targeting both DNA and RNA substrates.Sequences of regulatory ncRNAs can also help establish epigenetic alterations such as histone acetylation/deacetylation, DNA/histone methylation, etc. within the nucleus by bringing in chromatin remodelling agents that are known to change transcriptional activity [24,25].Based on their length, ncRNAs are arbitrarily divided into two groups: small ncRNAs (scnRNAs, < 200 nts) and long ncRNAs (lncRNAs, > 200 nts) [26].Depending on where they are in relation to genes that code for proteins, lncRNAs can also be categorized as sense, antisense, intronic, or intergenic [27].With regards to insecticide resistance in insects, lncRNAs that were found to be differentially expressed during the larval stage development of resistant Plutella xylostella genotypes [28] and uniquely differentially expressed during the egg to adult moth stages in Bt-toxin resistant strains of the same insect [29].Similarly, the expression of the lncRNAs in P. xylostella was linked to the expression of the cytochrome P450, the ATP-binding cassette (ABC) transporter and the esterase genes involved in resistance to chlorantraniliprole insecticide [27].Moreover, some long intergenic non-coding RNAs were overexpressed in deltamethrin-resistant larvae of Plutella xylostella exposed to deltamethrin [28].ncRNAs are intriguing candidates to study when organisms are exposed to insecticides and other toxicants since they are involved in pathways linked to responses to cellular stress [28,30].The genes for ribosomal proteins, such as L39 [31], S4 [32], L22 [33], and S29 [34], have been found to be associated with the resistance mechanism of Culex mosquitoes.
In the malaria-endemic region of western Kenya, there has been a resurgence of endophilic An. funestus and increased 20fold over a decade ago [9,35].The resurgence of this vector was partly attributed to resistance to pyrethroids used in ITN impregnation [36].As the country is aiming to achieve the malaria elimination goal by 2030, it is very crucial to have a comprehensive understanding of the resistant pro le of this important, re-emerged vector to inform stakeholders of the right choice of control strategy to adopt.To date, there have been few investigations on An. funestus susceptibility to insecticides in Kenya.The initial study on An. funestus susceptibility to insecticides from two study areas in western Kenya was reported in 2007 [37] and even though the species were not identi ed using molecular techniques, previously identi ed Anopheles species from the same areas revealed that only An. funestus was present [38].Later in western Kenya, seven adults An. funestus were sampled and their F1 progenies' susceptibility to insecticides revealed that they were susceptible to DDT but resistant to permethrin [39].Further study in Kisumu in the lowland area of western Kenya has shown that An. funestus is resistant to pyrethroids (deltamethrin and permethrin) with overexpression CYP6P9a and CYP6P9b responsible for pyrethroid resistance [15].A recent study in the same Kisumu, using microarray for transcriptome analysis has revealed that overexpression of cytochrome P450s notably, CYP4H18, CYP6M7, CYP9K1, CYP4C36 and CYP4H17 in pyrethroidresistant An. funestus population [40].The use of microarrays can only be used with the gene families that have been identi ed on the array, and they only give information on relative expression levels.The RNA-seq technology offers single nucleotide level resolution, absolute rather than relative gene expression pro le, and a comprehensive view of the transcriptome in a speci c state [41].
In this study, we examined the insecticide resistance pro le of An. funestus across ve sites in four counties in western Kenya and elucidated the molecular mechanisms of resistance using RNA-seq.Our results provide new novel insights into insecticide resistance at the molecular level in this important malaria vector, which has received limited attention, and could help in designing effective control strategies.

Mosquito sorting and identi cation
Live mosquitoes were sorted by separating male mosquitoes from the females and Culex spp from Anopheles.Later Anopheles mosquitoes were morphologically identi ed as An.funestus s.l and An.gambiae s.l following morphological and taxonomic keys [42,43].

Raising of F1 progenies.
An. funestus blood-fed, gravid and half gravid were put into cages to lay eggs on wet lter papers.The F 0 females in these physiological states were fed on 10% sugar solutions soaked in cotton wool and laying pads/Petri dishes.After laying the eggs, the eggs were allowed to hatch into larvae.The larvae were put in a pan containing spring water and were fed with the larvae feed, tetramin until they matured to a pupal stage where they were transferred into cages to emerge into adults.

Insecticide susceptibility tests
The F1 adult female mosquitoes aged between 3-5 days old were used for the bioassay.The insecticide susceptibility test was carried out following the standard insecticide tube test method developed by the WHO [17].The mortality was scored 24 hours post-exposure after maintaining under standard laboratory conditions at a temperature of 27 ˚C ± 2 ˚C and a relative humidity of 75% ± 10%.

PBO synergist bioassays
After establishing pyrethroid resistance in the An.funestus population, a synergist bioassay was conducted with PBOimpregnated papers to determine the role of P450 monooxygenases in pyrethroid resistance.The PBO inhibits these enzymes' activity in insects including mosquitoes.The female mosquito samples (F1 progenies) were pre-exposed to 4% PBO for 1 hour before they were immediately exposed to the pyrethroids (0.05% deltamethrin and 0.75% permethrin).

Preparation of samples for molecular and transcriptome analysis
Surviving resistant An. funestus samples after exposure to the insecticides (permethrin and deltamethrin) and unexposed (An.funestus F1 progenies samples that were not exposed to any insecticides) were killed immediately by keeping them in a deep freezer for about 10 minutes until they were completely knockdown.Samples were immediately stored in 0.5 ml Eppendorf tubes with RNALater and were immediately frozen at -80˚C for subsequent molecular and whole transcriptome analysis.

DNA extraction and molecular identi cation of species
DNA was extracted from the legs of each stored mosquito specimen using the Chelex®-100 method [44] and was transferred into pre-labelled 1.5ml storage vials and stored at -20˚C for molecular analysis.An. funestus-speci c PCR was conducted to con rm species using the species-speci c primers (ITS2A/FUN) in the internal transcribed spacer region (ITS2) on the ribosomal DNA [45,46].Species-speci c primers for An.funestus ( 5) and universal primer ( 5) were used.A nal volume of 12.5 µl of PCR mixture containing 1µl of genomic DNA, 6.5 µl DreamTaq Green PCR Master Mix (2x), 0.5 µl of each of the primers and 4.0 µl of PCR water.Genomic DNA ampli cation was performed using a T100 thermal cycler (Biorad).The PCR conditions include initial denaturation at 95 ˚C for 3 seconds, denaturation of 94 ˚C for 30 seconds, annealing at 55 ˚C for 30 seconds for 34 cycles, extension at 72 ˚C for 45 seconds and nal extension at 72 ˚C for 6 seconds.The DNA bands were visualized using the agarose gel electrophoresis.

RNA extraction
Total RNA was extracted from a pool of ten mosquitoes for each group (pyrethroids-resistant group and unexposed to pyrethroids).Total RNA was isolated and puri ed from the whole mosquito using the ZYMO Quick-RNA miniprep kit [47].
The details of the pooled samples is shown in Table 1.

Bioassay data analysis
The mortality rate of the sample tested was expressed as the total number of dead An. funestus mosquitoes of all the replicates exposed to a particular insecticide and expressed this as the percentage of all the population exposed to that insecticide.Abbott's formula was used to correct mortality if the mortality at 24 hours in the control tube was between 5% and 20%.Following the WHO criteria [17] for determining insecticide resistance in the malaria vector population, a population is classi ed as susceptible when the mortality is between 98-100%, resistant when mortality is less than 90% and suspected resistance when mortality is between 90-97%.

Quality control and differential expression analysis
Upon obtaining paired-end sequence reads from the sequencing centre (range 46,981,760-96640964 total reads), they were checked for quality using FASTQC (v0.11.5) [48] and cleaned to remove adapters (Table 1).Trimmomatic module (v.0.39) [49] was used to remove the Illumina adapters (TruSeq3-PE-2) that were used to construct the library, which resulted in the selection of reads that were more than 50 bp and a Phred-Quality-Score greater than 20 for downstream analysis.The resulting reads were con rmed to be of acceptable quality by running FASTQC (v0.11.3) before they were aligned to the reference genome.The Anopheles funestus FUMOZ genome from VectorBase (AfunF3.53) was used as the reference genome and this was aligned using the HISAT2 v2.2.1, which involved building the HISAT2 index le for the genome [50].To reduce the size of the output SAM tools from the alignment output, they were converted to BAM les, sorted, and indexed using SAMtools v1. 10 (Li et al., 2009).The sorted and indexed les were used as input for the Htseqcount reads, which were created using the module htseq-count (v.0.6.1) as described [52].A reference gene transfer format was used to count the number of alignment mapping to each gene based on union and intersection-strict [52].Gene expression values were normalized using Relative Log Expression (RLE) from DESeq2.Expression abundance between different treatments (pyrethroids-resistant group vs pyrethroid susceptible FANG colony of An. funestus mosquitoes [53] and unexposed vs the susceptible FANG colony) for the study sites (Teso, Port Victoria, Siaya and Kombewa) was determined using DESeq2 (v.1.18.0) [54].The susceptible FANG colony raw sequence data was retrieved from the GenBank

(accession: ERR981209-ERR981211). A correlation of gene expression between biological replicates was calculated by
Pearson's correlation as suggested before [55], while the Benjamini-Hochberg method was applied in calibrating the p-value to decrease chances of false positives [56].Therefore, differential expression between treatment groups was considered signi cant if the p-value < 0.05 and the fold change (FC) > 1.5 [57].Hierarchical clustering analysis was applied to cluster genes exhibiting similar expression patterns/levels, while the Gene Ontology (GO) from the GO database (http://geneontology.org/) was utilized for functional analysis of the differentially expressed genes to establish their biological pro les ( 52).

Phenotypic resistance pro le in western Kenya
Pyrethroid resistance was observed in all the sites with an average mortality rate (MR) of 57.6%.Port Victoria had the highest level of resistance to permethrin (MR = 53%) and deltamethrin (MR = 11%) pyrethroids.Teso had the lowest level of resistance to permethrin (MR = 70%) and deltamethrin (MR = 87%).Resistance to DDT was observed only in Kombewa (MR = 89%) and Port Victoria (MR = 85%).However, after samples were pre-exposed to the synergist, PBO, high susceptibility (> 98%) to the pyrethroids (deltamethrin and permethrin) was observed in all the sites except for Port Victoria where suspected resistance (96%) was observed for PBO + deltamethrin.In addition to the pyrethroid resistance, resistance to DDT was observed in Kombewa (89%) and Port Victoria (85%) (Table 2).Notwithstanding, a suspected resistance to DDT was observed in Siaya (93%) and Teso (92%).An. funestus was, however, fully susceptible to pirimiphos-methyl (0.25%) in all the sites.N: number of mosquitoes exposed to the insecticide 3.2 Differentially expressed genes between groups.
After quality control and the elimination of genes with low read counts, differential expression analysis was carried out on the transcripts.Three pairwise comparisons were performed: resistant versus susceptible, resistant versus unexposed (control) and unexposed versus susceptible.The resistant versus unexposed comparison helps to account for the induction of transcription during the pyrethroid exposure; genes were ltered by analyzing their expression pro les in the susceptible An. funestus population, under the assumption that constitutive resistance genes will be signi cantly differentially expressed between both survivors of the bioassay and the unexposed eld F1 progenies when compared to the susceptible FANG colony.The volcano plots (Fig. 2) showed the downregulation and upregulation of the genes between resistant vs susceptible (Fig. 2A), resistant vs unexposed (Fig. 2B) and unexposed vs susceptible (Fig. 2C).There was a clear distinction between overexpression of genes in resistant vs susceptible (Fig. 2A) and unexposed vs susceptible (Fig. 2C) group comparisons.However, there was no difference in the gene overexpression between the resistant and unexposed group comparisons.To summarize the expression pattern of genes in the sample groups, principal component analysis (PCA) and heatmap were used.The PCA plots indicate a representation of differences in the sample groups (resistant, unexposed and susceptible).The samples in the resistant and unexposed groups clustered together to the left-hand side away from the susceptible counterparts indicating the similarity between them (Fig. 3).The susceptible FANG group clusters towards the left side away from the resistant and the unexposed groups (Fig. 3).The heatmap revealed that there was an obvious grouping of the samples into resistant, unexposed and susceptible.
Dissimilarity in the gene expression levels was noticed between the groups.The overall gene expression pro le indicates a higher level of expression in the resistant and unexposed sample groups compared to the susceptible (Fig. 4).Moreover, most of the genes were highly expressed in the Kombewa-resistant (Kr01, Kr02 and Kr03) samples.This was followed by the resistant samples from Port Victoria, Siaya and Teso.However, low levels of gene expression were observed in the susceptible samples (ERR981209, ERR981210 and ERR981211).Comparison using the Venn diagrams, 33 genes (n = 14176) were differentially expressed in all the comparisons [resistant vs susceptible (R-S), resistant vs unexposed (R-C) and vs unexposed and susceptible (C-S)] (Fig. 5A).However, 953, 35 and 455 common genes were differentially expressed in only R-S, R-C and C-S respectively.More genes (1597) were signi cantly differentially expressed between R-S and C-S comparisons compared to the other comparisons.This was followed by 87 differentially expressed genes observed between C-S and R-C comparisons and 43 differentially expressed genes between R-S and R-C comparisons (Fig. 5A).Most of the downregulated genes were found in the R-S and this was followed by C-S (Fig. 5B).Nine Hundred and forty-nine (949) genes were downregulated between R-S and C-S comparisons (Fig. 5B).
Figure 5: Venn diagram comparing upregulated and downregulated genes between-group comparisons.A indicates upregulated genes between the groups and B indicates downregulated genes between groups.R-S: eld-resistant population that survived pyrethroid exposure vs susceptible colony, R-C: eld-resistant population that survived pyrethroid exposure vs unexposed (control) eld population and C-S: unexposed (control) eld population vs susceptible colony.

Differentially expressed ncRNAs linked to pyrethroid Resistance.
Differentially expressed ncRNAs between resistant vs susceptible (R-S) and unexposed vs susceptible (C-S) were determined by a fold change FC > 1.5 and FDR < 0.05.The whole transcriptome analysis shows that ncRNAs constituted 67%, the highest proportion of the gene families involved in pyrethroid resistance (Fig. 6).This was followed by IMPs (10%), CYPs (6%), CPs (5%), OPs (4%), GSTs (3%), UGTs (2%), ABCs (2%) and COEs (1%) (Fig. 6).IMPs: Imipenemase, CYPs: Cytochrome P450s, CPs: cuticular proteins, OPs: olfactory proteins, GSTs: Glutathione Stransferases, UGTs: UDP-glycosyltransferases, ABCs: ATP-binding cassettes, COEs: carboxylesterases and ncRNA: noncoding RNA The main ncRNAs that were overexpressed are in the resistant vs susceptible (R-S) and unexposed vs susceptible (C-S) are Metazoa_SR, RNaseP_nu, U3_1, Arthropod_7S, LSU_rRNA_eukarya_, SSU_rRNA_eukarya_2, LSU_rRNA_eukarya_13, SSU_rRNA_eukarya_46, LSU_rRNA_eukarya_2, LSU_rRNA_eukarya_3, SSU_rRNA_eukarya_15, LSU_rRNA_eukarya_5, LSU_rRNA_eukarya_6, SSU_rRNA_eukarya_164, SSU_rRNA_eukarya_19, SSU_rRNA_eukarya_200, LSU_rRNA_eukarya_155, LSU_rRNA_eukarya_17, LSU_rRNA_eukarya_17, LSU_rRNA_eukarya_214 and RNase_MRP (Table 3).Similarly, to identify the main genes in the enzyme families responsible for high pyrethroid metabolic resistance, FDR < 0.05 and a fold change FC > 1.5 were used.The main enzyme families identi ed are the cytochrome P450s, GSTs, salivary gland proteins, Peptidase S1 domain-containing proteins, UGTs and sulfotransferases (Table 4).However, most of these genes were moderately differentially expressed.The ndings indicate that in western Kenya, different genes within these enzyme families were responsible for resistance (Table 4).The top cytochrome P450 enzymes are moderately overexpressed in the An.funestus in western Kenya were CYP6P9, CYP6P9, CYP6N, CYP6N, CYP9J, CYP49A, CYP6P, AFUN02089, AFUN01936, CYP9K, CYP304B.These genes were overexpressed in resistant vs susceptible and unexposed/control vs susceptible group comparisons.However, CYP304C1 and CYP315A1 were overexpressed only in the resistant vs susceptible comparison (Table 4).Among the GSTs, the overexpressed genes in resistant vs susceptible and unexposed/control vs susceptible groups are GSTD, GSTT, GSTE, GSTD, and GSTD3 were overexpressed only in the resistant vs susceptible comparison (Table 4).AFUN02142, AFUN021428 and AFUN019106 were the only cuticular proteins that were overexpressed in the resistant vs susceptible comparison.The differential expression analyses revealed that some of these UGTs were overexpressed in the resistant vs susceptible and unexposed/control vs susceptible groups comparisons (Table 4).These include UGT302A, UGT310B, UGT308D, UGT306A3 and AFUN003620.AFUN016205 and AFUN016207 were the sulfotransferases that were overexpressed in the An.funestus population from western Kenya (Table 4).The summary of the RNA-seq data set for the FC and P-values of each gene is presented in Additional le 1. R: resistant eld mosquito population that survived the pyrethroid exposure, S: susceptible FANG colony, C: unexposed/control eld mosquito population, FC: fold change, NS: not signi cant 3.5 Gene ontology analysis of the differentially expressed genes GO term annotation pathways analysis was employed to elucidate the biological functions and signalling pathways that may be regulated by the differentially expressed genes in An. funestus.Our ndings revealed that these genes were engaged in a wide variety of biological functions and signalling pathways.Detailed GO enrichment for the differentially expressed genes in ontologies of cellular components, biological processes, and molecular function is represented in Fig. 7.
By GO annotation, the differentially expressed genes in pyrethroid-resistant An. funestus were enriched mostly in cellular macromolecule metabolic processes, cytoplasm, cellular protein metabolic processes and gene expression (Fig. 7).
Figure 7: Gene Ontology (GO) enrichment analysis of the differentially expressed genes.
The x-axis indicates the gene count/number of genes while the y-axis indicates the enriched terms.The colour is used to distinguish at different levels.

Discussion
The successful implementation and development of insecticide resistance management measures depends on elucidating the mechanisms underlying resistance in malaria vectors.In this study, we have characterized the phenotypic resistance pro le of An. funestus and the molecular basis of pyrethroid resistance in western Kenya.This is one of the most comprehensive studies on the An.funestus susceptibility status to pyrethroids and DDT in western Kenya.
Our study revealed a high level of pyrethroid resistance across western Kenya although resistance levels vary from site to site.In addition, resistance to DDT has been detected in Kombewa and Port Victoria.This con rmed a previous study in East Africa including western Kenya which reported widespread pyrethroid resistance in the An.funestus population [15].The rise of multiple resistance of An. funestus was also con rmed in a previous study in western Kenya [15], Benin, west Africa [59] and Malawi, southern Africa [60].An. funestus was, however, fully susceptible to pirimiphos methyl, the organophosphate in all the study sites.This is congruent with a previous study in Tanzania where a full susceptibility of this vector to pirimiphos methyl in Tanzania was reported [12].This is an indication that this insecticide can still be maintained for IRS programs in western Kenya.The preexposure of samples to the synergist, PBO has shown that An. funestus was fully susceptible to the pyrethroids in all the study sites except Port Victoria where 96% mortality was observed for the PBO + deltamethrin.This implies that the metabolic resistance mechanism (cytochrome P450 monoxygenases) was fully involved in insecticide resistance in the An.funestus in these sites but partially involved in Port Victoria [17].Other mechanism(s) might be contributing to pyrethroid resistance in Port Victoria leading to that site having the highest level of resistance compared to the other sites.
An. funestus has no kdr markers for resistance [61] hence metabolic resistance mechanism through overexpression of detoxi cation genes plays a crucial role in insecticide resistance [62,63].In this study, we have identi ed the top twenty ncRNAs that were differentially expressed in resistant and unexposed eld populations of An. funestus from western Kenya.Although their mechanisms of pyrethroid resistance in An. funestus is unknown, they could be playing a role in regulating the expression of pyrethroid-resistant metabolic genes in the An.funestus resistant populations.Our ndings add up to a body of evidence which hypothesised that ncRNAs play roles in insecticide resistance development [28].In general, the biological roles of ncRNAs in detoxi cation and insecticide resistance pathways are poorly understood.
However, few studies have reported that some ncRNAs (notably microRNAs) interfered with the expression of insecticidedetoxifying enzymes.For instance, MiR-2b-3p has been proposed to potentially suppress the cytochrome P450 9f2 .Non-coding rRNAs constitute over 80% of the total cellular RNA in mosquitoes [70].In this study, the majority of rRNA was removed using the RiboZero Plus kit.Typically, the literature indicates that rRNA comprises anywhere from 1 to 20% of the nal rRNA-depleted sequencing libraries [71].Our results demonstrate that our RNA-seq method can effectively detect and quantify both coding and non-coding RNA.
Two different ncRNA-based insect management approaches have been proposed following these ndings: (a) using biodegradable ncRNA-insecticide solutions to control insects, [72,73] and (b) using metabolic engineering techniques to nd and take advantage of target species' ncRNA-associated signalling pathways [74].However, the molecular mechanisms behind the functioning of ncRNAs in detoxi cation and insecticide resistance signalling pathways are still not clear, despite the mounting body of evidence suggesting these molecules are signi cant regulators of insect development [28,69,75].This is because research in this area is still in its early stages.The straightforward CRISPR-Cas9 genome editing technique has the potential to generate novel understandings of the roles of regulatory ncRNA sequences as well as ncRNA-based techniques targeted at managing insects including disease vectors.Moreover, using inhibitors to target speci c ncRNAs might interfere with the expression levels and reduce or reverse insecticide resistance.Thus, ncRNAs could be potential targets for vector control in the future.Recently, Oberemok et al. [76] proposed an innovative strategy to tackle insecticide resistance and create safer compounds.Their method employs synthetic DNA oligomers to disrupt gene expression by targeting ribosomal RNA (rRNA) rather than messenger RNA (mRNA).Because rRNA makes up 80% of cellular RNA and is more plentiful and stable than mRNA, it represents a promising target for DNA antisense oligonucleotide (ASO) interventions.This strategy seeks to offer a more e cient and enduring approach to combating insecticide resistance.
An. funestus is a notorious vector of human malaria in Africa and has contributed to over 90% of all malaria transmission in some parts of eastern and southern African regions [12,13].It is noteworthy that the outcome of our study represents an advancement in the molecular basis of insecticide resistance in An. funestus population.Our comprehension of the functional importance of ncRNAs in insecticide resistance pathways could enable the creation of ncRNA-based vector control techniques to control An. funestus.

Conclusions
An. funestus population is highly resistant to pyrethroids in western Kenya with Port Victoria recording the highest levels of resistance to the type I and type II pyrethroids.However, preexposure to PBO synergists recorded high susceptibility to the pyrethroids except in Port Victoria.We have shown for the rst time that insecticide resistance in An. funestus is linked to the expression of ncRNAs hence a better understanding of these molecular events could help to develop resistance management strategies for future malaria control.

Figure 1 :
Figure 1: Map of study sites where mosquitoes were sampled in western Kenya.The software ArcGIS Pro 2.6 was used to create the map.Map sources: USGS, ESRI, and CGIAR (www.esri.com)

Figure 2 .
Figure 2. Volcano plot indicating upregulation and downregulation for resistant vs susceptible (A), resistant vs unexposed (B) and unexposed (control) vs susceptible (C).The X-axis indicates the log2 fold-change-positive and negative values are up and down-regulated respectively relative to the susceptible group in A and C. The Y-axis indicates -log 10 of the adjusted P-value (FDR) (-log 10 FDR values > 200 for A, > 9 for B and > 80 for C).In each volcano plot, genes that are overexpressed in the population are > 0 on the x-axis.P-values of < 0.05 are indicated by the horizontal line, while 2-fold expression differences are indicated by vertical dotted lines.

Figure 3 :
Figure 3: A principal component analysis showing the gene expression pattern of the sample groups relative to the susceptible group.

Figure 4 .
Figure 4. Heatmap indicating the expression of genes in the sample groups relative to the susceptible group.

Figure 6 :
Figure 6: Pie chart showing the proportion of gene family involving pyrethroid resistance.

Figures
Figures

Figure 3 A
Figure 3

Table 1
RNA-seq read ltering and mapping statistics.
The quality of each total RNA sample was assessed using High Sensitivity RNA Tapestation (Agilent Technologies Inc., California, USA) and quanti ed with Qubit 3.0 RNA HS assay (ThermoFisher, Massachusetts, USA).Ribosomal RNA was depleted with Ribo-Zero Plus rRNA Removal Kit (Illumina Inc., California, USA).Samples were then heated, fragmented and randomly primed according to the manufacturer's recommendation.The rst strand was synthesized with the Protoscript II Reverse Transcriptase with a longer extension period, approximately 30 minutes at 42 C.All remaining steps for library construction were performed according to the NEBNext® Ultra™ II Directional RNA Library Prep Kit for Illumina® (New England BioLabs Inc., Massachusetts, USA).Final libraries quantity was assessed by Qubit 3.0 (ThermoFisher, Massachusetts, USA) and quality was assessed by TapeStation D1000 ScreenTape (Agilent Technologies Inc., California, USA).The nal library size was about 350bp with an insert size of about 200bp.Illumina® 8-nt dual indices were used.

Table 2
Mortality rate of An. funestus exposed to different insecticides and synergists (24-hr post-exposure).

Table 3 :
List of the top non-coding RNA (ncRNA) involving pyrethroid resistance 3.4 Differentially expressed metabolic genes associated with pyrethroid resistance.

Table 4
List of the top genes of immunity, metabolic, cuticle and olfactory R: resistant eld mosquito population that survived the pyrethroid exposure, S: susceptible FANG colony, C: unexposed/control eld mosquito population, FC: fold change, NS: not signi cant This study was approved by Maseno University's Ethics Review Committee (MUERC/00778/19).Verbal consent was sought from owners of households before mosquitoes were collected inside the living rooms.Consent for publication: Not applicable.Tables Table 3 is available in the Supplementary Files section.