Comparative proteomics reveals complex insecticides-associated resistance mechanism of Culex pipiens pallens Coquillett

Background: Mosquito control based on chemical insecticides is considered as an important element in the current global strategies for the control of mosquito-borne diseases. Unfortunately, the development of insecticide resistance of important vector mosquito species jeopardizes the effectiveness of insecticide-based mosquito control. As opposed to target site resistance, other mechanisms are far from being fully understood. Results: Susceptible strain of Cx. pipiens pallen showed elevated resistance levels to after 25 generations insecticide-selected, through bioinformatics analysis allowed detecting 2,502 proteins, of which 1513 were differentially expression in insecticide-selected strains as compared to the susceptible strain. Finally, midgut differential expression protein proles and 62 proteins were selected for verication of differential expression using parallel reaction monitoring strategy. Conclusions Signicant molecular resources were developed for Cx. pipiens pallen potential candidates involved in metabolic resistance as well as those participating in lower penetration or sequestration of insecticide. Global protein proles of change to three insecticide strains combined with midgut proles revealed multiple insecticide resistance mechanisms operate simultaneously in resistant insects of Cx. pipiens pallens. Future research that is targeted towards RNA interference on the identied metabolic targets such as cuticular, cytochrome P450s and glutathione S-transferase proteins could lay the foundation for a better understanding of the genetic basis of insecticide resistance in Cx. pipiens pallen.

Conclusions Signi cant molecular resources were developed for Cx. pipiens pallen potential candidates involved in metabolic resistance as well as those participating in lower penetration or sequestration of insecticide. Global protein pro les of change to three insecticide strains combined with midgut pro les revealed multiple insecticide resistance mechanisms operate simultaneously in resistant insects of Cx. pipiens pallens. Future research that is targeted towards RNA interference on the identi ed metabolic targets such as cuticular, cytochrome P450s and glutathione S-transferase proteins could lay the foundation for a better understanding of the genetic basis of insecticide resistance in Cx. pipiens pallen.

Background
Mosquitoes and mosquito-borne diseases continue to pose a huge public health burden worldwide, malaria, lymphatic lariasis, dengue, chikungunya, and west nile virus cause signi cant medical and economic impacts that disproportionately affect developing countries [1][2][3][4].
Since the 1950s, chemical insecticides have been massively used for controlling mosquito populations.
However, long-term intensive and widespread overuse or misuse use of insecticides induced intense selection pressure that has led to the development and subsequent intensi cation of various genetically modulated resistance mechanisms in mosquitoes [5][6][7], and has resulted in mosquito-borne diseases on the rise and outbreaks of mosquito-related diseases in recent years [8,9]. Today, with the widespread development of resistance in mosquitoes to the most commonly used insecticides, insecticide resistance was regarded as the most serious threat to the control of mosquitoes and mosquito-borne diseases [10].
Nowadays, vector control of mosquitoes remains crucial to reduce disease transmission, and has long been a critical part of the global strategy to manage mosquito-associated diseases, the application of chemical insecticides are the most important component in mosquitoes control [11], so, characterizing molecular mechanisms underlying resistance is a key step for improving resistance management strategies.
Proteins, as the primary functional molecules, the executor to perform the gene function in many different physiological processes, currently, large-scale expression pro ling and screening for resistance-related proteins in mosquitoes are limited, two-dimensional gel electrophoresis (2-DE) proteomics [32,33], and isobaric tags for relative and absolute quanti cation (iTRAQ) [34]. Today, with the advent of highly advanced proteomic platforms based on iTRAQ [35], such knowledge gap can be overcome by high throughput sequencing approaches, which can generates concomitantly protein expression, ultimately provides a direct insight into the activity of relevant proteins.
In the present study, iTRAQ labeling coupled with liquid chromatography/tandem mass spectrometric (LC-MS/MS) analysis was used to determine changes in protein amounts associated with adaptation to three insecticides (the pyrethroid cypermethrin, the propoxur, and dimethyl-dichloro-vinyl-phosphate) from distinct chemical families in the susceptible mosquito Cx. pipiens pallens. Results are discussed in regards of known and new putative adaptive mechanisms conferring insecticide resistance in mosquitoes, which improve our understanding of mechanisms developed by mosquitoes to resist insecticides.

Results
Insecticide resistance levels After 25 generations of larval selection with the insecticides cypermethrin, propoxur and DDV, bioassays revealed a constitutive increased resistance of each selected strain to its respective insecticide as compared to the parental susceptible strain (Table 1). Resistance levels of the Cx_pro and Cx_ddv strains were moderate but signi cant (11.34 fold and 8.17 fold respectively). Although signi cant, the resistance level of the Cx_cym strain to cypermethrin susceptible was considerably higher (243.00 fold).
Protein identi cation and differentially expressed proteins (DEP) screening A total of 2,502 protein were identi ed in two replicates of the LC-MS/MS experiments with a high con dence of peptide selection (FDR=0.01). Among them, 1513 protein obtained quantitative message, all the proteins were completed annotated with the accession number for each identi ed protein arranged by the Gene Ontology, KEGG pathway, Domain, Cluster of Orthologous Groups (COG) function catalogues and Subcellular Location analyses (Additional le 1: Table S1). Our analysis revealed signi cant changes in proteins in these Cx_cym, Cx_pro, Cx_ddv, and Cx_s strains, there were 164 proteins with signi cant differential expression (fold change > 1.2 and p value < 0.05) between the Cx_cym and Cx_s control groups, including 54 up-regulated and 110 down-regulated proteins. 156 proteins changed dynamically between the Cx_pro and Cx_s control groups, including 59 up-regulated and 97 down-regulated proteins, whereas 81 proteins changed between the Cx_ddv and Cx_s control groups, with 31 up and 50 downregulated proteins, respectively.

Bioinformatics analyses of the altered proteins
To better understand the characteristics of the altered proteins in response to insecticide selection, GO ontology was performed to analyze the upregulated and downregulated proteins among Cx_cym, Cx_pro, Cx_ddv, and Cx_s strains of Cx. pipiens pallens ( Figure 1). As shown, the most enriched GO terms regarding molecular function of the upregulated proteins were related to the structural constituent of cuticle, structural molecule activity, lipid transporter activity, chitin binding, transferase activity, transferring acyl groups and odorant binding, which had a higher fold enrichment (log2) and Fisher's exact test p-value (-log10). Whereas with the downregulated proteins, structural constituent of cuticle were the most enriched GO terms.
Hierarchical clustering was employed to characterize changes in the Cx_cym, Cx_pro, Cx_ddv, and Cx_s strains, all quanti ed proteins were classi ed into 7 clusters based on their expression level, revealed clearly distinguishable increasing or decreasing expression ( Figure 2). In KEGG enrichment analysis, drug metabolism -cytochrome P450 (map00982) and metabolism of xenobiotics by cytochrome P450 (map00980) pathways exist in response to insecticide selection conditions, implicating P450s are the important regulatory enzyme associated with metabolic resistance in Cx. pipiens pallens under insecticide selection. Furthermore, the differential changes in glutathione S-transferase (B0W6D0), microsomal glutathione s-transferase (B0X075) and glutathione-s-transferase theta (B0XGK3) proteins in both pathways suggest a critical regulatory role associated with P450, these ndings indicate that glutathione S-transferase may play an important role in the regulation of metabolic resistance pathway in Cx. pipiens pallens ( Figure 3A and Figure 3B).
Comparing the pattern of expression of different proteins among Cx_cym, Cx_pro, Cx_ddv, and Cx_s strains of Cx. pipiens pallens, all of these proteins are signi cantly differentially expressed in resistant strains when compared to their expression in the susceptible strain, almost all protein enrichment changes at the peak and the lowest valley were in Cx_cym, Cx_pro and Cx_ddv, indicating that most of these protein showed an increase or decrease in expression following the insecticides selection from susceptible to resistance strains in Cx. pipiens pallens, the insecticide selection alters the normal developmental pathway into an alternative one ( Figure 4).

Prediction of the PPI network under insecticide selection
Because the protein-protein interaction plays important roles in most cellular biological processes, the PPI network were predicted under insecticides selection conditions using the STRING (www.string-db.org) online PPI prediction software, to further analyze the functional correlation of differentially expressed proteins between each group ( Figure 6, Additional le 2: Figure S3-S4).

Validation of proteomics data by parallel reaction monitoring
To validate the accuracy and reproducibility of the proteomics analysis results, parallel reaction monitoring (PRM)-quanti able strategy analyses were performed, 65 proteins selected as those that may play important roles in pyrethroid resistance, the differentially expressed levels of these proteins were quanti ed in Cx_cym, Cx_pro, Cx_ddv and Cx_s strains of Cx. pipiens pallens were veri cation and validation, all proteins, peptides, and measured peptide peak areas are listed (Additional le 3: Table S2).

Differential protein expression in the midgut of Cx. pipiens pallen induced by Bti
A total of 564 protein were identi ed with 559 protein obtained quantitative message, there were 66 proteins with signi cant differential expression (fold change > 1.2 and p value < 0.05) between the Cx_bti and Cx_nbti control groups, including 33 up-regulated and 33 down-regulated proteins (Additional le 4).
Most of the proteins coding for cuticular, cytoskeleton related proteins, metabolic enzymes are enrichment differential changes in the digestive tract linings dissected from Cx_bti strain of Cx. pipiens pallens (Additional le 2: Figures S5-S8). Given our interest in the role of P450s in metabolic resistance, we further pro led the differential protein expression of Cytochrome P450 3A19 (B0X758) in the midgut of third instar larvae of Cx. pipiens pallen, P450 protein showed higher expression levels in Cx_bti strains.

Discussion
Cx. pipiens pallenst is the primary vector of the Japanese encephalitis virus and lariasis in China, that has developed resistance to several insecticides used for mosquitoes control. To date, the lack of genomics data available for this species has hampered characterisation of the molecular mechanisms underlying resistance. A profound understanding of the molecular mechanisms regulating insecticides resistance and development may aid in the control of this mosquitoes by facilitating the development of more sustainable and environmentally-friendly approaches. In this study, the proteomics of three distinct chemical families insecticides selected strains (the Cx_cym, Cx_pro, Cx_ddv groups) and susceptible strain of mosquitoe Cx. pipiens pallens were compared.
The insect cuticle, also known as the exoskeleton, is the outermost part of the insect body, it serves a variety of functions such as sensory perception of the environment, means of locomotion [36], maintain the physical structure of the organism, the foregut, hindgut, tracheal system and apodemes [37], and protection from desiccation, moreover, the cuticle also serves as the rst and major barrier against external adverse compounds penetration [38]. So, the cuticle is a major route of insecticide penetration in insects, insect cuticles are composed of cuticular protein, chitin and lipids, in fact, the properties of insect cuticle including permeability to pyrethroids are in uenced not only by the chitin sclerotization, construction, and hydration, but also by the regular combinations of different cuticular proteins and their arrangements [38,39], recent protein-level and gene-level analyses have demonstrated a surprising diversity of cuticular proteins within and amongst species [40].
In previous studies, some cuticle genes have been identi ed overexpressed in pyrethroid resistant strains of several mosquito species in An. stephensi [41], An. funestus [14] and Cx. pipiens pallens [34] and bed bug, Cimex lectularius [39,42,43]. The abundance of the protein in limb cuticle correlates nicely with the > 2-fold increased abundance of their transcript in pyrethroid resistance An. gambiae [44]. The uptake of organophosphate can be reduced by a thickening or a change in the chemical composition of the cuticle in organophosphate-resistant strains of Culex quinquefasciatus [45] and Culex tarsalis [46], pyrethroid resistant An. funestus and An. gambiae do indeed have thicker cuticle than sensitive forms [37,47], in the German cockroach, Blattella germanica decreased penetration of labeled insecticides in (1R)-transpermethrin-selected strain of the house y and fenvalerate resistance [48,49], low level of cross resistance to many insecticides observed in the An. stephensi strain might be partially explained by reduced penetration [50], the epicuticle hydrocarbon lipids were enriched in pyrethroid resistance An. gambiae [51] and An. arabiensis [52] female mosquitoes.
Our experiments indicate that GO terms representing proteins of structural constituent of cuticle, structural molecule activity, chitin binding proteins were affected by insecticide selection in Cx. pipiens pallens, as a relatively large number of cuticle-associated proteins with a diversity of functions showed highly differential abundance in cypermethrin resistant strain, propoxur resistant strain, DDV resistant strain and susceptible strain (Fig. 1). And in our study, chitin-binding type R&R consensus, insect cuticle protein, actin, conserved site, actin/actin-like conserved site, and actin family were enriched domains for downregulated proteins in Cx_cym vs. Cx_s. In Cx_pro vs. Cx_s, chitin binding domain was upregulated proteins, whereas insect cuticle protein, chitin-binding type R&R consensus, actin, conserved site, actin/actin-like conserved site, actin family were downregulated proteins, respectively. In Cx_ddv vs. Cx_s, the enriched upregulated proteins was ferredoxin-like fold domain, chitin-binding type R&R consensus, whereas downregulated proteins were insect cuticle protein, chitin-binding type R&R consensus, adult cuticle protein 1 (Fig. 5).
Change of cuticular protein determines the cuticle composition and performance is the basis for explaining the performance of cuticle-based insecticide resistance. Therefore, other alterations of cuticle structure and composition such as cuticular remodelling may be also involved in cuticular resistance. Taken together, our observations con rm the notion that cuticle proteins show altered expression and cellular cytoskeletal rearrangement, suggesting forti cation of structural components in response to insecticide selection and leading to increased insecticide resistance in Cx. pipiens pallens. Although how cuticle proteins are involved in the process of cuticle alterations for slowing down the penetration of insecticides is puzzling, this result provides useful information for us to conduct further studies on the cuticular protein and cuticle.
Metabolic resistance occurs via subtle alterations in either the protein levels or the activity of detoxi cation enzyme that resist insecticides [53,54], such as P450s, glutathione S-transferases (GSTs), P450s possess catalytic roles in the metabolism of many essential endogenous molecules, and allows insects to metabolize insecticides at a higher rate [31,55], up-or down-regulation of P450 genes may be responsible for the detoxi cation of insecticides and the homeostatic response of mosquitoes to changes in cellular environment [56]. Cx. pipiens pallens, like most insect species, metabolises xenobiotics such as insecticides using a suite of detoxi cation enzymes, in this study functional annotation of cytochrome P450 role was identi ed by proteomic studies directly, several detoxi cation proteins cytochrome P450 (B0XDA9), cytochrome P450 3A19(B0 × 758), cytochrome P450 4g15 (B0WQV0, B0XHB7), cytochrome P450 9b1 (B0W6Y4), cytochrome P450 monooxygenase CYP9M10 (E5E7H3) and cytochrome P450 12b1 (B0WTS3) have been identi ed that are differentially expressed in insecticide resistant populations in Cx. pipiens pallens. In this study, P450 proteins whose expression levels changed signi cantly belonged primarily to the different families. In addition, data from many other studies are highlighting a key role for P450 in insecticide detoxi cation in Anopheles funestus [14], Culex quinquefasciatus [57,58]. Elevated levels of P450 activity are frequently observed in pyrethroid-resistant malaria vectors in Africa [59][60][61][62][63][64][65][66][67], in line with other insect pests such as the white y [68] and cockroach [69], over-expression of cytochrome P450s contribute to resistance to neonicotinoid insecticides in the hemipterans Bemisia tabaci and Myzus persicae [70,71]. Furthermore, it has been suggested that these elevated P450s also confer crossresistance to carbamates [60].
Beyond insecticide metabolism and resistance, P450 enzymes playing essential roles in hormone metabolism or cuticular hydrocarbon synthesis [72]. Recently, the P450 enzyme responsible for the oxidative decarbonylation of long-chain fatty aldehydes [73] was shown to be a member of the CYP4G subfamily [74], which are restricted to insects, with clearly identi able sequence peculiarities [74], interesting, in this present study, cytochrome P450 4g15 (B0WQV0, B0XHB7) right belong to CYP4G family.
The tissue speci c expression of a protein normally is related to its function in that tissue [75], the functional signi cance of the digestive tract linings in regards to Cx. pipiens pallens is especially important, all these tissues contain epidermis and cuticle, as the digestive tract lining constitutes the major part of their body, when Cx. pipiens pallens larvae come in contact with or consume insecticides, they may develop resistance by modi cation in the insect cuticle or digestive tract linings that prevent or reduce the rate of penetration. In this study, cuticular, cytoskeleton related proteins are enrichment differential expressions in the digestive tract linings compared from Cx_bti strain of Cx. pipiens pallens (Additional le 3, Additional le 4: Figures S5-S8), higher levels of proteins when compared to the levels were detected for cytochrome P450 3A19 (B0 × 758) of the midgut in Cx_bti vs. Cx_nbti, suggesting a potential role in metabolic resistance of P450. Once insecticide enters the Cx. pipiens pallens larvae, enhanced metabolic detoxi cation could decrease the concentration of insecticides before they reach the target site, the midgut-associated Cx. pipiens pallens P450s consistent with their role in xenobiotic metabolism.

Conclusions
In summary, these results signi cantly increase the molecular resources available for the study of mosquitoes and we have identi ed several candidate proteins potentially involved in different phases of insecticide metabolism as well as those participating in other modes of insecticide resistance in Cx. pipiens pallens. Signi cant differentially expressed proteins across biologically variant samples were revealed by iTRAQ proteomics approach. Speci cally, the high occurrence of differential expression cuticular, cytoskeleton related proteins, and the expression patterns of P450s in midgut(cuticular) tissue support previous studies suggesting represent unique sites for penetration resistance (cuticle) as well as metabolic resistance (P450s) in Cx. pipiens pallens.
A common phenomenon of insecticide resistance is that multiple mechanisms operate simultaneously in resistant mosquotes of Cx. pipiens pallens. Typically, a combination of diverse mechanisms provides signi cantly higher levels of resistance than one individual mechanism. In addition it has provided further support for the role of P450 genes in insecticides resistant while suggesting that other gene families such as cuticular genes could also be involved. Future functional studies on cuticular and P450s proteins could lay the foundation for identifying hot-spots for insecticide resistance in Cx. pipiens pallens, which could provide the basis for developing effective management strategies.

Mosquito selection with insecticide and larvicidal bioassay
The laboratory Tangkou strain of Cx. pipiens pallens, originating from Tangkou Village (Jining prefecture, Shandong Province), has been colonized in standard insectary conditions without exposure to any insecticides since 1960, fully susceptible to insecticides, all the susceptible mosquito populations were reared at a constant room temperature of approximately 28 ℃ and 75% relative humidity with a photophase of 14 h and a scotophase of 10 h, adult mosquitoes were provided with 10% sugar solution, was used as a parental strain to repeatedly selected for 25 generations to produce three independent resistant strains at the larval stage with the pyrethroid insecticide (cypermethrin), the carbamate insecticide (propoxur) and the organophosphates insecticide (DDV) at the median lethal concentration (LC50) [34], a detailed selection procedure was described previously [76,77]. Technical-grade formulations of cypermethrin, propoxur and DDV were employed for the larval bioassay study, different concentrations of these insecticides were prepared using denatured alcohol as solvent (98 ml of absolute alcohol + 2 ml of methyl ethyl ketone) as described by Shetty [78].The cypermethrin, propoxur and DDV selected strains were de ned as cypermethrin-resistant (Cx_cym), propoxur-resistant (Cx_pro) and DDV-resistant (Cx_ddv) strains.
For validation of differently protein expression in the midgut of Cx. pipiens pallen, larvae were reared in round plastic tubs (diameter 0.6 m) lled with water with sh food twice daily, experimental larvae were randomly collected from several tubs to compensate for size differences and feeding history which are known to be in uenced by larval density, all laboratory experiments were carried out with third instar larvae of laboratory-reared, 50 third instar larvae were exposed for 5 h to 0.021 mg/l of Bacillus thuringiensis var. israelensis (Cx_bti, Shandong Lukang Sheryl Pharmaceutical Co., Ltd, 3000 ITU/mg), the dosage following the standard testing procedures for microbial tests [79], we chose a 5 h treatment time to allow larvae to ingest the Bacillus thuringiensis var. israelensis without severely affecting their behavior and without causing evident damage to the intestinal tissue at the microscopic histological level. In the control group (Cx_nbti), larvae were exposed to water under the same conditions. Larvae exposed to the insecticide were washed in water and midgut pools were dissected and transferred to 1.5-mL microcentrifuge tubes containing 1mm phenylmethylsulfonyl uoride (PMSF) and stored at −80 ℃ before use.

Extraction of protein and proteolysis
Each sample (the Cx_cym, Cx_pro, Cx_ddv and Cx_s strains) were extracted separately from 5 egg rafts, 5 fourth-instar larvae, 5 pupae, 2-3 days old 5 adult males, 5 adult females (without blood feeding) and midgut pools. Brie y, each samples were completely homogenized in protein extraction buffer (8 M urea, 2 mM EDTA, 10 mM dithiothreitol (DTT) and 1% Amresco Protease Inhibitor Cocktail), then samples were lysed at room temperature for 3 min of sonication, and centrifuged for 10 min at 13000 rpm, the supernatants was subjected to precipitated with chilled acetone for 3 h at -20 ℃. The resulting protein concentration of each sample was measured using Bradford Protein Assay Kit. Then, a total of 100 μg protein/sample was combined from egg, larvae, pupae, adult male and female stage (each for 20 μg), reduced with 10 mM DTT for 60 min at 37 ℃, followed by adding 25 mM iodoacetamide (IAA) to each protein samples, then alkylation by 45 min in darkness at room temperature, subsequently, 100 mM tetraethylammonium bromide (TEAB) was added to reduce the concentration of urea to < 2 M. Finally, samples were digested with trypsin (enzyme-to-substrate ration of 1:50) overnight at 37 ℃, and an additional second digestion was done for 4 h to ensure complete cleavage.
After trypsin digestion, each peptides samples desalting was performed by using Strata X SPE column, dried, and re-suspended in 25 μL 500 mM TEAB and labeled with 8-plex iTRAQ kit. Each dried and labeled peptide sample was reconstituted by using HPLC solution A (2% ACN, pH 10) and fractionated by high-pH reverse-phase HPLC on a Waters Bridge Peptide BEH C18 (130 Å, 3.5 μm, 4.6 × 250 mm). Loaded peptides were eluted with 2% to 98% acetonitrile gradient buffer solution at pH 10 in 60 fractions at a speed of 0.5 ml/min over 88 min. A total of 20 fractions were combined and each fraction was desalted by using Ziptip C18 (Merck Millipore, Ziptip Pipette Tips 10μL). Samples fractions were dried on vacuum concentrator and stored at -20 ℃ pending MS analyses.
High-resolution LC-MS/MS analysis LC-MS/MS experiments were performed on a NanoLC 1000 LC-MS/MS liquid chromatography system (Thermo) equipped with a Proxeon EASY-nLC 1000 coupled to LTQ-Orbitrap Elite. Brie y, peptides were acidi ed in 0.1% formic acid, were trapped in a 3 μm 100Å Acclaim PepMap®100C18 (Thermo) 75 μm×2 cm in 100% solvent A (0.1 M acetic acid in water), at a ow rate of 5 μL/min was used. Then, the analytical separation was performed on 2 μm 100 Å Acclaim PepMap® RSLC C1850 μm × 15 cm using a linear gradient of solvent B (98% ACN with 0.1% formic acid). The gradient ran at a ow rate of 250 nL /min as follows: 0-60 min from 10% to 35% solvent B, 35% to 50% in 10min and ramping to 100% solvent B. Mass spectrometry was operated in a data-dependent mode, automatically switching between MS and MS/MS (at mass range of m/z 350 to 1800 and resolution of 60,000) using repetitively full MS spectra scans followed by collision induces dissociation (HCD, at 38 normalized collision energy). The 20 most intense precursors were selected for subsequent decision tree-based ion trap HCD fragmentation in the MS survey scan and the following settings were used: collision energy: 38% above threshold, ion count: 300 with 30.0s dynamic exclusion. Full width at half maximum (FHMW) at 400 m/z using an AGC setting of 1e6 ions, xed rst mass was set as 100 m/z.

Mass spectrometry data analysis
The mass spectrometry MS/MS raw data were searched against all Culex quinquefasciatus protein in the Culex Taxonomy database, downloaded from Uniprot Database using Sequest software integration in Proteome Discoverer (version 1.3, Thermo Scienti c) with the following parameters: mass tolerance set to 20 ppm, product ion tolerance set to 0.02 Da, only trypsin sequences allowed. Oxidation, and protein Nterminus acetylation were accepted as variable modi cations and carbamidomethylation was accepted as static modi cation, a maximum of two miscleavages was permitted, peptide-and protein-level false discovery rates (FDRs) were ltered to 1%.

Functional classi cation of proteins
Functional annotation and classi cation of all identi ed proteins were determined by using the Blast2GO program against the non-redundant NCBI protein database. The gene function databases complied for the Gene Ontology (GO), the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, and Cluster of Orthologous Groups (COG) of proteins were used to infer putative functions of identi ed proteins. Detailed information can be found on http://www.geneontology.org. Signi cant GO terms (p < 0.05) were mapped with the REViGO online tool (http://revigo.irb.hr), which removes redundant GO terms and visualizes the semantic similarity of remaining terms [80]. KEGG Pathway is part of the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, which is a reference database for pathway mapping, for KEGG Pathway analyses, the main biochemical metabolism and signal transduction pathways of identi ed proteins were extracted using the Search pathway tool in the Kegg Mapper platform (http://www.genome.jp/kegg/mapper.html) [81][82][83]. COG is the database for orthologous protein classi cation. Identi ed proteins were compared with the COG protein database to predict the function of proteins and conduct statistical analyses.
Functional enrichment, protein expression pattern and function clustering Enrichment statistics of GO, KEGG pathway and Domain was done using Fisher's exact test, correction for multiple hypothesis testing was carried out under standard false discovery rate control methods and the pathways with a corrected p value < 0.05 were considered as the most signi cant pathways.
Expression-based and functional enrichment-based clustering for different protein groups was used to explore potential relationships between different protein groups at special protein function. Firstly, we collated all the protein groups obtained after functional enrichment analysis along with their P values.
Secondly, we sorted those categories enriched in at least one of the protein groups with a P value <0.05. This ltered P value matrix was transformed by the function x = −log10 (P value). Thirdly, z-transformed applies on x values for each functional category and z scores were clustered by one-way hierarchical clustering (Euclidean distance, average linkage clustering). Finally, cluster membership was visualized by a heat map using the "pheatmap" function from the R-package. Prediction of proteine-protein interaction (PPI) was done using the STRING database (version 10.0) and visualized using Cytoscape software (version 3.4.0) [84,85].
Parallel reaction monitoring analysis A label-free targeted parallel reaction monitoring (PRM) method[86] was adopted to verify the reliability of our label-based proteomics. A total of 62 differentially expressed proteins, including one internal standard housekeeping protein, were chosen. For the Cx_cym, Cx_pro, Cx_ddv and Cx_s strains, equal amounts of protein samples were analyzed for semi-quantitative measurements, and each strain underwent three replications. Peak areas were extracted from PRM mass spectrum data using Skyline software[87].

Declarations
Ethics approval and consent to participate

Consent for publication
All authors have given consent for publication.

Availability of data and materials
All data generated or analyzed during this study were included in this published article and its additional les.

Competing Interests
Author Tao Li was employed by the company Nanning MHelix ProTech Co., Ltd., the remaining authors declare that the research was conducted in the absence of any commercial or nancial relationships that could be construed as a potential con ict of interest.

Description Of Additional Files
Additional le1 Table S1: Global protein pro les among insecticide selected of Culex pipiens pallens using iTRAQ.
Additional le 2 Figure S1-S8 Figure S1: Annotated KEGG map for oxidative phosphorylation, map00190 The blue border belongs to the background protein, and the black border indicates the protein was not identi ed in this experiment. The red/green colour in the gure marks the differentially expressed proteins detected in this experiment, in which red represents upregulated proteins and green represents downregulated proteins. Half red and half green indicates the protein is both upregulated and downregulated (the same meaning as in this manuscript).   Clustering of protein expression patterns and functional enrichment. Differential protein relative expression matrix heatmap of the three insecticide resistance groups; the color in the gure indicates the relative expression level of the protein in the sample. The cluster of more abundantly expressed proteins is highlighted in red, and the green color represents a lower expression level, color bars represent speci c Page 24/27 expression abundance. The total number of observed spectra assigned to each protein in each insecticide selection was used as the basis for clustering. background proteins, and the white color indicates proteins not identi ed in this experiment. The red/green colors in the gure indicate to the differentially expressed proteins detected in this study, with red representing upregulated proteins and green downregulated proteins. Half red and half green indicates both upregulated and downregulated proteins for that gene product.

Figure 4
Comparison of protein expression pattern clusters between Cx_cym, Cx_pro, Cx_ddv and Cx_s in Cx.
pipiens pallens For the expression trend line graph of each sub-cluster, the x-axis is the comparison sample group, and the y-axis is the relative expression level of the protein in the group of samples. Each line in the gure represents a protein, and the different color representations show the relationship between relative expression and the mean value. Each graph shows one type of expression pattern, a trend that re ects changes in the expression of this group of proteins.3