RNA-seq screening of the CPGs in Culex pipiens pallens among cypermethrin-resistant populations

Background Long-lasting overdependence on insecticides has led to the rapid spread of pyrethroid resistance in mosquito vectors, which is of great concern to the general public. There are many studies on metabolic resistance and target resistance, but fewer studies have been conducted on cuticle resistance and behaviour resistance. The cuticle of mosquitoes has been hypothesized to play a role in insecticide resistance by reducing penetration or sequestering insecticides. Methods We used RNA sequencing (RNA-seq) to analyse the transcriptome of cypermethrin-resistant and cypermethrin-susceptible strains of Culex pipiens pallens . Sequenced 6 samples using an Illumina HiSeq platform, and generated approximately 6.66 Gb bases from each sample on average. Mapping the sequenced reads to a reference genome and reconstructing the transcripts, through gene expression analysis, we detected differentially expressed genes (DEGs) among the samples. Followed Gene Ontology (GO) classification and functional enrichment. Finally, we screened the genes of cuticle proteins associated with drug resistance throughout the genome, selected the significant DEGs with a log2 fold change>3.0 and Padj<0.05, and applied real-time fluorescence quantitative PCR to verify the DEGs. Results We obtained 13,517 novel transcripts, of which 8,653 were previously unknown splicing events for known genes, 665 were novel coding transcripts without any known features, and 4,199 were long noncoding RNA. A total of 1035, 944, and 657 genes were upregulated in comparisons between samples, and 2680, 1215, and 975 genes were downregulated in comparisons between samples. Finally, among all samples, 167 genes upregulated, and 145 genes downregulated. The GO genes; cellular component, 149 genes; and biological process, 272 genes. The expression of XM_001863852 and XM_001845881 in resistant strains of Culex pipiens pallens was lower than that in the laboratory sensitive strain, with fold changes in expression of 0.177 and 0.548, respectively; the expression of the XM_001845883.1 in the resistant strain was higher than that in the susceptible strain, and a 2.281-fold change in expression. The results provide a for the furthermore, could provide new for vector

. Because most mosquito-borne diseases (except JEV) currently do not have effective vaccines or therapeutic drugs, effective control of vector mosquitoes is the main measure to prevent mosquito-borne diseases. Among the measures, chemical control is the primary method of vector management [3,4].Pyrethroid insecticides are synthetic analogs of naturally occurring pyrethrins from Chrysanthemum spp. [5]. Due to their low mammalian toxicity, high insecticidal activity, fast action and ease of decomposition in the environment, pyrethroids are currently a dominant class of insecticides used globally against mosquitoes and other human disease vectors [6].However, with the long-term use of insecticides, mosquito resistance is becoming increasingly more serious, which not only greatly reduces the control effect of existing insecticides but also makes the development of new drugs more difficult. The emergence of insecticide resistance has become a major obstacle for controlling mosquito-borne diseases [7].
Insecticides, pyrethroids in particular, remain a mainstay for the control of these important vectors. In this paper, we review what is known about the levels, mechanisms and fitness costs of pyrethroid resistance in Cx. pipiens. Pyrethroid resistance in Cx. pipiens is a global problem, and resistance ratios of up to 7000fold have been found in larvae of field collected mosquitoes [8].
High levels of resistance to pyrethroids in Culex mosquitoes have been widely reported [9,10]. Previous surveys have shown that Culex pipiens pallens/Cx. quinquefasciatus in southern China have different levels of resistance to pyrethroid insecticides. In Hainan and other provinces, the resistance level has reached several thousand fold [11]. The development of mosquito resistance has had an important impact on the control of Culex pipiens pallens.
Mosquito resistance mechanisms include metabolic resistance, target resistance, cuticle resistance, and behavioural resistance [12]. Metabolic resistance refers to the degradation, isolation, or transportation/excretion of insecticides from cells prior to binding the target. Metabolic resistance results from increased detoxification caused by the overexpression of or conformational changes in the enzymes involved in chemical insecticide metabolism, sequestration, and excretion.
P450-monooxygenases, glutathione S-transferases, and carboxy/cholinesterases are the main enzymes involved in this process [13][14][15][16].Target-site resistance, or mutations in target binding sites for insecticides, iscaused by a modification of the chemical insecticide site of action, reducing or preventing insecticide binding at that site. Mutations in the voltage sensitive sodium channel (Vssc) gene are the most common causes of target-site resistance [17]. Behavioural resistance results from selection pressure of mosquitoes under long-term exposure to pesticides, and mosquitoes thus show a series of behavioural changes to avoid pesticides. For example, long-term application of indoor residual spraying (IRS) and insecticidetreated nets (ITNS) caused mosquitoes to change from endophagy to exophagy and from endophily to exophily and caused peak bloodsucking to change from late night to dusk [18].
Cuticle thickening is implicated in insecticide resistance by reducing the uptake of the insecticide that reaches the target site in response to the modification of chemical composition of the cuticle [19]. However, the mechanism remains poorly understood, and its importance in Aedes species is yet to be confirmed [14,20,21].A study revealed that this mechanism may play a major role in the development of resistance where it normally happens simultaneously with other mechanism(s) [22], causing resistance to single or multiple insecticides [23]. It has been reviewed elsewhere that cuticle thickening is associated with metabolic detoxification whereby thicker cuticle causes gradual insecticide absorption rate that will increase the effectiveness of metabolic detoxification in Anopheles funestus [24]. Moreover, it is crucial to take note that insects with cuticular resistance will display resistance level of not more than 3-fold in comparison to susceptible insects, but the cooccurrence of other resistance mechanism will lead to a surge in insecticide resistance level markedly [25]. This is demonstrated by Anopheles gambiae in Benin [26] through which overexpression of cuticular genes and P450 genes gave rise to a relatively high resistance level.
High-throughput sequencing, i.e., RNA-seq, has become the main choice for measuring expression levels [27]. RNA-seq can be performed without prior knowledge of the reference or sequence of interest and allows a wide variety of applications such as the 'de novo' reconstruction of a transcriptome (without a reference genome), the evaluation of nucleotide variations, and the evaluation of methylation patterns [28]. RNA-seq technology has some advantages over the cDNA microarrays, such as a high level of data reproducibility, which reduces the number of technical replicates for experiments. In addition, RNA-seq allows identification and quantification of the expression of isoforms and unknown transcripts [29]. We used RNA-seq to analyse the transcriptome of cypermethrin-resistant and cypermethrin-susceptible strains of Culex pipiens pallens and obtained several cuticle protein genes that were differentially expressed in laboratory sensitive and resistant lines. However, these genes were unverified in the genome. We determined the sensitivity and resistance of mosquitoes to cypermethrin by realtime PCR and Centers for Disease Control and Prevention (CDC)Bottle Bioassay. In addition to the identification of the cypermethrin cuticle resistance genes, the relationship between the identified cuticle protein genes and mosquito resistance was verified to establish the specific mechanism of cypermethrin resistance in Culex pipiens pallens and provide new ideas for mosquito control and treatment.

Mosquito sample collection
We collected laboratory sensitive and resistant strains of Culex pipiens pallens in the following developmental stages: I, II, III, and IV instar larvae, pupa and female Culex pipiens pallens 3 days after hatching that had not fed on blood. A total of 200 mg of Culex pipiens pallens was collected at each developmental stage and placed into a 1.5 ml Eppendorf (EP) tube, to which 150 µl of TRIzol lysis buffer was added to soak the mosquitoes. The samples were quickly stored in a -80 °C freezer, and RNA was extracted and sent to the BGI group for transcriptome sequencing. The resistant strain was screened from sensitive strains in our laboratory in accordance with the larvae dipping method recommended by the WHO.

RNA sequencing
Three RNA-seq libraries per population (3 biological replicates) were prepared using TRIzol reagent (Thermo Fisher Scientific, Inc., Waltham, MA, U.S.A.). The extracted RNA was treated with RNase-free DNase (Qiagen GmbH, Hilden, Germany) and purified using an RNeasy MinElute Cleanup Kit (Qiagen GmbH) to remove DNA.
The amount of total RNA was measured in a NanoDrop 2000 (Thermo Fisher Scientific, Inc.). The quality of the extracted RNA was verified by agarose gel electrophoresis. Subsequently, the reaction systems to synthesize the first and second strands of cDNA were constructed, and after the second strand of cDNA was synthesized, the ends of the double-stranded cDNA were blunt ended with the EcoRI restriction sequence. After terminal phosphorylation and XhoI digestion, dscDNA was recovered using a recovery kit. An Agilent 2100 Bioanalyzer and ABI StepOnePlus Real-Time PCR System were then used for quality tests, and the Illumina HiSeq platform was used for RNA-seq after dscDNA quality was confirmed.
Sequenced reads were assigned to each sample (unplexing), and adaptors were removed. Read quality was assessed for each sample using FastQC. Reads were then filtered based on their length, pairing and quality using Trimmomatic [30] with the following parameters: Leading, 25; Trailing, 25; Minlen, 60; and Slidingwindow, 4-25. Only paired reads were kept. Reads were then mapped to the Culex pipiens quinquefasciatus genome using Tophat2 [31] with the following parameters: don't report discordant pair alignments; final read mismatches = 3; The DEseq2 and PossionDis algorithms were used to perform DEG detection.
DEseq2 is differential analysis software based on the principle of negative binomial distribution, and the analysis was conducted according to the method used in Michael I et al. [32]. The PossionDis difference analysis algorithm is based on the Poisson distribution model, and the analysis was conducted according to the method described in Audic S et al. [33].

Screening and verification of DEGs between sensitive and resistant strains
Based on the results of the abovementioned transcriptome sequencing analysis, we analysed the genes related to cuticle proteins among the susceptible and resistant strain DEGs and identified 3 mRNAs with a fold change greater than 2 (-log2Ratio ≥ 1) and FDR ≤ 0.001 for subsequent real-time PCR verification.
Referring to the Culex pipiens pallens genome data in the gene library, realtime quantitative PCR primers were designed using Primer Premier 5 software and the nucleotide sequences of the selected mRNAs, and β-actin was used as the quantitative internal mRNA reference. The primers were synthesized by Shenzhen BGI. The base sequences of the specific primers are provided in Table 1. RNA was extracted using TRIzol reagent, and cDNA was synthesized. cDNA was synthesized using a RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, Inc.) using oligo(dT)18. qRT-PCR was performed using a CFX96 Touch (Bio-Rad Laboratories, Inc., Hercules, CA, U.S.A.). Twenty-five nanograms of cDNA and 500 nM of each forward and reverse primer were used in each reaction. The relative expression of each gene in resistant and susceptible mosquitoes was calculated by the ΔΔCt method [34] using actin (ADIR001186-RA) as a control. Real-time PCR data were analysed using reliability simulation tool (REST) software and hypothesis testing; other data were expressed as the mean ± standard deviation (XJ ± S) and analysed with STATA7.0 software. Student's T test was performed for comparisons between groups, and p<0.05 was used as a basis for determining statistical significance.  Tables 2, 3).  statistical information for all samples is provided in Table 5. We then analysed the site information for each SNP and INDEL, as shown in Figure 1 and Figure 2.
[See supplementary files for Table 5.] Numbers, functional categorization and pathway analysis of DEGs DEGs were obtained by comparing the expression levels of differential genes between the sample groups. The results are shown in Figure 3.
DEGs were compared with the GO library for functional classification. GO The calculation results are provided in Tables 6 and 7.  It can be seen from Table 7

Discussion
It has been reported that at least 504 insect species worldwide exhibit chemical insecticide resistance [38], that more than 109 (subspecies) vector mosquitoes are resistant to one or more insecticides [39], and that the number of species and the Based on the high chemical resistance of Culex pipiens pallens in Shandong Province [40], this study used high-throughput transcriptome sequencing technology to qualitatively and quantitatively study cypermethrin-resistant strains of Culex pipiens pallens and identified the cuticle protein genes responsible for the drugresistant phenotype. Real-time quantitative PCR were used to verify the identified cuticle protein genes and to confirm whether there is a high correlation between cuticle protein genes and the drug resistance phenotype.
Mosquito resistance generally includes 4 types: metabolic resistance, target resistance, cuticle resistance, and behavioural resistance. In 1963, a survey of houseflies (Fannia canicularis) identified the cause of insecticide resistance. It was found that the penetration of chemical pesticides in resistant lines was slower than in sensitive lines, suggesting that slower penetration could be the cause for dichlorodiphenyltrichloroethane (DDT) and pyrethroid resistance [41]. Studies have shown that pyrethroid-resistant female Anopheles sinensis have thicker cuticles than do sensitive female Anopheles sinensis and that female mosquitoes also have thicker cuticles than do male mosquitoes [42]. In addition, Lily et al. [43] confirmed that cuticle thickening is present in pyrethroid-resistant strains of Cimex lectularius.
Cuticle analysis by electron microscopy and characterization of lipid extracts showed that resistant mosquitoes had a thicker outer skin layer and a higher hydrocarbon content (approximately 29%) [44]. Multiple studies have confirmed that the cuticle of resistant lines is much thicker than that of sensitive lines [43,44].
The cuticular protein (CP) family was first discovered in 2007 by tandem mass spectrometry analysis of epidermal exfoliation from Anopheles gambiae [ 45]. Most gene family members have the prefix CPLC (cuticular protein of low complexity), and such members often play an important role in protein-protein interaction networks [46,47].
The above evidence suggests that insect cuticle proteins play an indispensable role in mosquito resistance. The expression of the target genes XM_001863852 and XM_001845881 in the Culex pipiens pallens laboratory resistant strain was lower than that in the laboratory sensitive strain, and the fold changes in expression were 0.177 and 0.548, respectively, while the expression of XM_001845883.1 in the resistant strain was higher than that in the sensitive strain, and with a 2.281-fold change in expression. GO function analysis indicated that all 3 were genes responsible for cuticle structural components, and an NR database comparison also showed that these genes code cuticle proteins in Culex pipiens. In fact, when we screened the genes, we initially identified 3 different genes in the XM_001845    Pathway classification of DEGs The X axis represents the number of DEGs, and the Y axis rep Figure 5 Pathway classification of DEGs The X axis represents the number of DEGs, and the Y axis rep

Supplementary Files
This is a list of supplementary files associated with the primary manuscript. Click to download.