RNA-seq screening of cuticle protein genes in Culex pipiens pallens among cypermethrin-resistant populations

Background A long-lasting overdependence on insecticides has led to the rapid spread of pyrethroid resistance in mosquito vectors, which poses a great risk to the general public. Although there are many studies on metabolic resistance and target resistance, few have investigated 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 . We 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 via 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 log 2 fold change > 3.0 and Padj < 0.05, and applied real-time fluorescence quantitative polymerase chain reaction (PCR) to verify the DEGs.

respectively; the expression of the XM_001845883.1 in the resistant strain was higher than that in the susceptible strain, with a 2.281-fold change in expression.

Conclusions
Our results provide a reference for resistance mechanisms via the mosquito cuticle as well as a new perspective for disease vector control.

Background
Culex pipiens pallens is the most common mosquito in northern urban areas and townships in China.
In addition to stings and bites, it is also the main vector of several arboviruses, such as West Nile virus (WNV), St. Louis encephalitis (SLE), Sindbis virus (SINV), Rift Valley fever (RVFV), and Japanese encephalitis (JEV), and is the main vector for lymphatic filariasis [1,2]. Historically, chemical control is the most commonly used measure to control vector mosquitoes [3,4]. Insecticides, particularly pyrethroids, due to their low mammalian toxicity, high insecticidal activity, fast action, and ease of decomposition in the environment, remain a mainstay for the control of mosquito vectors [5,6].
However, with the long-term use of pyrethroids, the resistance of mosquitoes to pyrethroids is increasing [7]. Many studies have found that pyrethroid resistance in Culex pipiens larvae is a global problem, with resistance ratios of up to 7000, 710, 370, and 18 for permethrin, deltamethrin, cypermethrin, and λ-cyhalothrin, respectively [8]. Culex pipiens pallens/Cx. quinquefasciatus in southern China show different levels of resistance to pyrethroid insecticides. In Hainan and other provinces, the resistance level has reached several-thousand-fold [9]. The JPal-per strain of Cx.
quinquefasciatus showed a marked resistance to permethrin (2500-fold compared to that of an insecticide susceptible strain) [ 10] and to other pyrethroids, such as phenothrin (2460-fold) and etofenprox (4160-fold) [11] during the larval stage. A recent study found that a cypermethrinresistant strain of Culex pipiens pallens, Coq, showed a 283.06-and 80.68-fold resistance to cypermethrin and permethrin compared to that of susceptible strains [12]. Shi et al. found that the resistance of Culex pipiens pallens to insecticide increased from generation to generation with consistent exposure to insecticides. They selected a mild selection strain and administered deltamethrin at a constant concentration of 0.05 ppm for 24 generations, and found that the level of resistance grew exponentially, with an increase in the resistance ratio of over 8-fold [13].
Mosquito resistance mechanisms include metabolic resistance, target resistance, cuticle resistance, and behavioural resistance [14]. Metabolic resistance refers to the degradation, isolation, or transportation/excretion of insecticides from cells prior to binding to 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. P450monooxygenases, glutathione S-transferases, and carboxy/cholinesterases are the main enzymes involved in this process [15][16][17][18]. Target-site resistance, or mutations in target binding sites for insecticides, is caused by a modification of the chemical insecticide site of action, which reduces or prevents insecticide binding at that site. Mutations in the voltage sensitive sodium channel (Vssc) gene are the most common causes of target-site resistance [19]. Behavioural resistance results from selection pressure of mosquitoes under long-term exposure to pesticides, such that mosquitoes show a series of behavioural changes to avoid pesticides. For example, the long-term application of indoor residual spraying (IRS) and insecticide-treated nets (ITNs) was found to cause mosquitoes to change from endophagy to exophagy and from endophily to exophily, resulting in the peak bloodsucking time throughout the day to change from late at night to dusk [20].
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 the chemical composition of the cuticle [21].
However, this mechanism remains poorly understood, and its importance in the Aedes species is yet to be confirmed [16,22,23]. A previous study found that this mechanism may play a major role in the development of resistance, where it normally occurs simultaneously with other mechanism(s) [24], causing resistance to single or multiple insecticides [25]. It has been reviewed elsewhere that cuticle thickening is associated with metabolic detoxification, whereby a thicker cuticle leads to a gradual insecticide absorption rate that will increase the effectiveness of metabolic detoxification in Anopheles funestus [26]. Moreover, it is crucial to note that insects with cuticular resistance will display a resistance level of not more than 3-fold in comparison to that of susceptible insects, but the co-occurrence of other resistance mechanisms will lead to a marked surge in the insecticide resistance level [27]. This is demonstrated by Anopheles gambiae in Benin [28], in which the overexpression of cuticular genes and P450 genes gave rise to a relatively high resistance level.
In order to screen the cuticle protein genes (CPGs) responsible for cypermethrin resistance in Culex pipiens pallens, we sequenced the transcriptome of cypermethrin-resistant strains and cypermethrinsusceptible strains. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed, as well as the functional analysis and real-time polymerase chain reaction (PCR) validation of three important differentially expressed genes (DEGs).
Our aim was to provide a reference for the cuticle resistance mechanism of Culex pipiens pallens and provide a novel perspective on mosquito control and management.

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 immediately stored in a -80 °C freezer, and the RNA was extracted and sent to the BGI group (Shenzhen, China) for transcriptome sequencing. The resistant strain was screened from sensitive strains in our laboratory in accordance with the larvae dipping method recommended by the World Health Organization (WHO).

RNA Sequencing
Three RNA-seq libraries per population (3 biological replicates) were prepared using TRIzol reagent (Thermo Fisher Scientific, Inc., Waltham, MA, USA). 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. After the second strand of cDNA was synthesized, the ends of the double-stranded cDNA (dscDNA) were blunt-ended using the EcoRI restriction sequence. After terminal phosphorylation and XhoI digestion, the dscDNA was recovered using a recovery kit. An Agilent 2100 Bioanalyzer and ABI StepOnePlus Real-Time PCR System were then used for quality tests. The Illumina HiSeq platform was used for RNA-seq after the quality of the dscDNA was confirmed.
Sequenced reads were assigned to each sample (unplexing), and adaptors were removed. The read quality was assessed for each sample using FastQC. Reads were then filtered based on their length, pairing, and quality using Trimmomatic [29] with the following parameters: Leading, 25; Trailing, 25; Minlen, 60; Slidingwindow, 4-25. Only paired reads were kept. The reads were then mapped to the Culex pipiens quinquefasciatus genome using Tophat2 [30] with the following parameters: do not report discordant pair alignments; final read mismatches = 3; intron length = 45-300000; use coverage search. Only read pairs mapping at a unique location (mapQ > 50) were retained. The quantification of the transcription levels was performed using the Cuffdiff2 module of Cufflinks The DESeq2 and PossionDis algorithms were used to perform DEG detection. DESeq2 is a differential analysis software based on the principle of negative binomial distribution, and the analysis was conducted according to the method reported by Michael et al. [31]. The PossionDis difference analysis algorithm is based on the Poisson distribution model, and the analysis was conducted according to the method described in Audic et al. [32].

Screening and verification of DEGs between sensitive and resistant strains
Based on the results of the 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 (-log 2 Ratio ≥ 1) and FDR ≤ 0.001 for subsequent real-time PCR verification.
Referring to the Culex pipiens pallens genome data in the gene library, real-time quantitative PCR primers were designed using Primer Premier 5 software and the nucleotide sequences of the selected mRNAs. β-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. The RNA was extracted using TRIzol reagent, and the cDNA was synthesized using a RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, Inc.) using oligo(dT)18. Quantitative reverse transcription PCR (qRT-PCR) was performed using a CFX96 Touch (Bio-Rad Laboratories, Inc., Hercules, CA, USA).
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 the resistant and susceptible mosquitoes was calculated using the ΔΔCt method [33] with actin (ADIR001186-RA) as a control. The 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 used to perform comparisons between groups. p < 0.05 was used as a basis for determining statistical significance.  Tables 2, 3).  [35]) to compare the reintegrated transcripts with the annotation information for the Culex pipiens genome. We selected transcripts with a class code type of u, i, o, and j as candidates for novel transcripts.
A total of 13,517 new transcripts were detected. Detailed statistical information is provided in Table 4.

Detection of single nucleotide polymorphisms (SNPs) and insertions and deletions (INDELs)
After comparing the clean reads to the Culex pipiens pallens genome, we used Genome Analysis Toolkit (GATK) [36] software to call each chromosome, identify SNPs and INDEL sites for each sample, and store the final results in variant call format (VCF). The SNP statistical information for all samples is provided in Table 5. We then analysed the site information for each SNP and INDEL, as shown in Fig. 1 and Fig. 2. Table 5 Single nucleotide polymorphisms (SNPs) Transition: Substitution of a purine with a purine or substitution of a pyrimidine with a pyrimidine; Transversion: Substitution of a purine with a pyrimidine.

Sample A-G C-T Transition A-C A-T C-G G-T Transversion Total
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 Fig. 3.

Validation of the 3 target cuticle protein genes
To ensure the amplification of the target genes and the housekeeping gene, we performed primer verification.
The results showed that the amplification curve for the primers was good and that the melting curve was  Tables 6 and 7 for details).  The calculation results are presented in Tables 6 and 7.
As shown in Table 7, the expression levels of the target genes XM_001863852 and XM_001845881 were similar between the sensitive and resistant strains of Culex pipiens pallens; the fold changes in expression were 0.177 and 0.548, respectively. The expression level of the target gene XM_001845883.1 in the resistant strain was higher than that in the sensitive strain, with a fold change in expression of 2.281.

Discussion
In 1963, a survey of houseflies (Fannia canicularis) identified the cause of insecticide resistance: the penetration of chemical pesticides in resistant lines was slower than in sensitive lines, suggesting that a slower penetration could be the cause for dichlorodiphenyltrichloroethane (DDT) and pyrethroid resistance [37]. The cuticular protein (CP) family was first discovered in 2007 by the tandem mass spectrometry analysis of epidermal exfoliation from Anopheles gambiae [38]. The majority of the gene family members have the prefix CPLC (cuticular protein of low complexity) and often play an important role in protein-protein interaction networks [39,40]. Studies have shown that pyrethroid-resistant female Anopheles sinensis have thicker cuticles than pyrethroid-sensitive female Anopheles sinensis, and that female mosquitoes also have thicker cuticles their male counterparts [26]. In addition, Lily et al. confirmed that cuticle thickening is present in the pyrethroid-resistant strains of Cimex lectularius [41]. Cuticle analysis by electron microscopy and the characterization of lipid extracts showed that resistant mosquitoes had a thicker outer skin layer and a higher hydrocarbon content (approximately 29%) [42].
These findings suggest 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 resistant strain was lower than that in the sensitive strain, with fold changes in expression of 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, with a 2.281-fold change in expression. GO function analysis indicated that all 3 were genes responsible for cuticle structural components, and a non-redundant (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 series: XM_001845880, XM_001845881, and XM_001845883.1. However, when designing the primers, the ones designed for these 3 were not ideal, with varying degrees of non-specific amplification. We compared the designed primer fragments to the BLAST database and found that the primers were highly consistent with an unknown conserved hypothetical protein. We wondered whether the 3 different genes in the XM_001845 series were different splicing isoforms of the same gene, and if so, whether such frequent splicing could promote resistance in mosquitoes.
After considering the experimental cost, experimental significance, and feasibility of the experiment, we selected 2 of the genes to conduct further validation, and the experimental results for XM_001845881 and XM_001845883.1 were opposite to that of each other. Therefore, we questioned whether the 2 genes had antagonistic effects in the genetic pathways leading to the formation of cuticle resistance in mosquitoes; for example, XM_001845883.1 is responsible for promoting the formation of cuticle resistance, while XM_001845881 is responsible for regulating the expression of other upstream or downstream cuticle protein genes, preventing the overexpression of other related cuticle proteins.
Since the 3 identified target genes in our study are novel, the original identification process was relatively complicated and innovative; however, the specific regulatory networks and the in vivo function of the respective genes will need to be further investigated.

Conclusions
These data provide transcriptomic information related to the resistance of Culex pipiens pallens and preliminarily verify the relationship between the identified cuticle protein genes and mosquito drug resistance, partially explaining the specific mechanism of mosquito cuticle resistance, providing a scientific basis for the study of new target insecticides and a novel perspective for the mosquitoes control and management.

Availability of data and materials
The datasets generated and analysed during the current study are available in the NCBI Sequence Read Archive (SRA) repository, the accession number is: PRJNA601003.

Authors' contributions
QQS, PC, CXZ, IJL, XS, XXG, HFW, YW, and HML conducted the sample collection and wrote the manuscript. HWW and MQG reviewed and edited the manuscript. All authors read and approved the final manuscript.

Ethics approval and consent to participate
Not applicable.

Consent for publication
Not applicable.      Pathway classification of DEGs. The X axis represents the number of DEGs, and the Y axis represents the KEGG pathway.