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 classification and functional enrichment of DEGs as follows: molecular function, 224 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.
Conclusions The results provide a reference for resistance mechanisms through the mosquito cuticle, furthermore, could provide a new perspective for disease vector control.
Figure 1
Figure 4
Figure 6
Figure 8
Figure 10
Overview of RNA-seq data
We selected 3 samples from different physiological stages of the sensitive and resistant strains, for a total of 6 samples, and each sample produced an average of 6.66 Gb of data. The sequenced clean reads were compared with the reference genome of Culex pipiens pallens, and the transcripts were reintegrated. A total of 13,517 new transcripts were detected, of which 8,653 were new alternative splicing isoforms of existing known protein-coding genes, 665 were transcripts of unknown protein-coding genes, and 4,199 were lncRNAs (see Tables 2, 3).
Table 2 Summary of differentially expressed genes
VS |
Upregulated |
Downregulated |
Cx_S_strain-VS-Cx_R_strain.DEseq2 |
167 |
145 |
Cx_S_strain-VS-Cx_R_strain.DEseq3 |
1035 |
2680 |
Cx_S_strain-VS-Cx_R_strain.DEseq4 |
944 |
1215 |
Cx_S_strain-VS-Cx_R_strain.DEseq5 |
657 |
975 |
Table 3 Summary of whole genome expression
Sample name |
Total gene number |
Number of known genes |
Number of novel genes |
Total transcript number |
Known transcript number |
Novel transcript number |
R_strain_1 |
14597 |
14036 |
561 |
19094 |
11413 |
7681 |
R_strain_2 |
14507 |
13921 |
586 |
19460 |
11415 |
8045 |
R_strain_3 |
14592 |
14016 |
576 |
19733 |
11533 |
8200 |
S_strain_1 |
14603 |
14040 |
563 |
19381 |
11440 |
7941 |
S_strain_2 |
14551 |
13980 |
571 |
19671 |
11562 |
8109 |
S_strain_3 |
14568 |
13999 |
569 |
19679 |
11501 |
8178 |
Prediction of new transcripts
After comparing the clean reads to the Culex pipiens pallens genome, we used StringTie [35] software to perform transcript reintegration for each sample and then used cuffmerge and cuffcompare software (both are packages in Cufflinks [36]) 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.
[Due to technical limitations, Table 4 could not be displayed here. Please see the supplementary files section to access the table.]
Detection of SNPs and INDELs
After comparing clean reads to the Culex pipiens pallens genome, we used Genome Analysis Toolkit (GATK) [37] software to call each chromosome, identify single nucleotide polymorphisms (SNPs) and insertion and deletion (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 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 included 3 major categories: molecular function, cellular composition and biological process. Subsequently, the 3 categories were separately classified. The results are provided in Figure 4.
The obtained DEGs were compared with the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, and the related metabolic pathways were identified. The results are provided in Figure 5.
Validation of the 3 target cuticle protein genes
To ensure 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 monomodal, and the electrophoresis results revealed specific target fragments. The average Ct value for the actin gene in the sensitive group was 25.905; the average Ct value of the actin gene in the resistant group was 26.227; the average Ct value of the XM_001863852 gene in the sensitive group was 26.813; the average Ct value of the XM_001863852 gene in the resistant group was 29.633; the average Ct value of the XM_001845883.1 gene in the sensitive group was 36.797; the average Ct value of the XM_001845883.1 gene in the resistant group was 35.93; the average Ct value of the XM_001845881 gene in the sensitive group was 32.647; and the average Ct value of the XM_001845881 gene in the resistant group was 33.837 (see Tables 6 and 7 for details).
Subsequently, we used the 2-△△Ct method to analyse the expression of target gene mRNA in the extracted RNA and normalized the result based on the housekeeping gene.
The specific calculation method is as follows:
△△Ct = (Ct target gene - Ct housekeeping gene) experimental group - (Ct target gene - Ct housekeeping gene) control group
The relative expression level of the target gene = 2^-△△Ct, which indicates the fold change in expression for the target gene in the experimental group relative to the control group.
The calculation results are provided in Tables 6 and 7.
Table 6 Quantitative PCR results
Target ID |
Target Name |
Sample ID |
Sample Name |
Ct Avg (SDM) |
Rel. Qty (SDM) |
T001 |
Actin |
S001 |
S (sensitive) 1 |
25.475 |
1.00E+00 |
T001 |
Actin |
S002 |
S 2 |
26.205 |
1.00E+00 |
T001 |
Actin |
S003 |
S 3 |
26.035 |
1.00E+00 |
T001 |
Actin |
S004 |
R (resistance) 1 |
25.37 |
1.00E+00 |
T001 |
Actin |
S005 |
R 2 |
26.19 |
1.00E+00 |
T001 |
Actin |
S006 |
R 3 |
27.12 |
1.00E+00 |
T002 |
XM_001863852 |
S001 |
S 1 |
31.79 |
1.00E+00 |
T002 |
XM_001863852 |
S002 |
S 2 |
24.00 |
3.67E+02 |
T002 |
XM_001863852 |
S003 |
S 3 |
24.65 |
2.08E+02 |
T002 |
XM_001863852 |
S004 |
R 1 |
30.235 |
2.73E+00 |
T002 |
XM_001863852 |
S005 |
R 2 |
30.61 |
3.72E+00 |
T002 |
XM_001863852 |
S006 |
R 3 |
28.055 |
4.16E+01 |
T003 |
XM_001845883.1 |
S001 |
S 1 |
36.19 |
1.00E+00 |
T003 |
XM_001845883.1 |
S002 |
S 2 |
37.07 |
9.01E-01 |
T003 |
XM_001845883.1 |
S003 |
S 3 |
37.13 |
7.68E-01 |
T003 |
XM_001845883.1 |
S004 |
R 1 |
35.43 |
1.58E+00 |
T003 |
XM_001845883.1 |
S005 |
R 2 |
35.94 |
1.95E+00 |
T003 |
XM_001845883.1 |
S006 |
R 3 |
36.42 |
2.67E+00 |
T004 |
XM_001845881 |
S001 |
S 1 |
36.48 |
1.00E+00 |
T004 |
XM_001845881 |
S002 |
S 2 |
29.33 |
2.36E+02 |
T004 |
XM_001845881 |
S003 |
S 3 |
32.13 |
3.01E+01 |
T004 |
XM_001845881 |
S004 |
R 1 |
33.04 |
1.01E+01 |
T004 |
XM_001845881 |
S005 |
R 2 |
35.15 |
4.13E+00 |
T004 |
XM_001845881 |
S006 |
R3 |
33.32 |
2.80E+01 |
Table 7 Calculation of the relative quantitative Ct values for target genes and the internal reference gene
|
β-actin |
XM_001863852 |
XM_001845883.1 |
XM_001845881 |
Average Ct value in the sensitive group |
25.905 |
26.813 |
36.797 |
32.647 |
Average Ct value in the resistant group |
26.227 |
29.633 |
35.93 |
33.837 |
△ Ct value in the sensitive group |
----- |
0.908 |
10.892 |
6.742 |
△ Ct value in the resistant group |
----- |
3.406 |
9.703 |
7.61 |
△△ Ct value |
----- |
2.498 |
-1.189 |
0.868 |
2^-△△ Ct value |
----- |
0.177 |
2.281 |
0.548 |
It can be seen from Table 7 that 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, and the fold change in expression was 2.281.
This is a list of supplementary files associated with this preprint. Click to download.
Posted 19 Dec, 2019
RNA-seq screening of the CPGs in Culex pipiens pallens among cypermethrin-resistant populations
Posted 19 Dec, 2019
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 classification and functional enrichment of DEGs as follows: molecular function, 224 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.
Conclusions The results provide a reference for resistance mechanisms through the mosquito cuticle, furthermore, could provide a new perspective for disease vector control.
Figure 1
Figure 4
Figure 6
Figure 8
Figure 10
Overview of RNA-seq data
We selected 3 samples from different physiological stages of the sensitive and resistant strains, for a total of 6 samples, and each sample produced an average of 6.66 Gb of data. The sequenced clean reads were compared with the reference genome of Culex pipiens pallens, and the transcripts were reintegrated. A total of 13,517 new transcripts were detected, of which 8,653 were new alternative splicing isoforms of existing known protein-coding genes, 665 were transcripts of unknown protein-coding genes, and 4,199 were lncRNAs (see Tables 2, 3).
Table 2 Summary of differentially expressed genes
VS |
Upregulated |
Downregulated |
Cx_S_strain-VS-Cx_R_strain.DEseq2 |
167 |
145 |
Cx_S_strain-VS-Cx_R_strain.DEseq3 |
1035 |
2680 |
Cx_S_strain-VS-Cx_R_strain.DEseq4 |
944 |
1215 |
Cx_S_strain-VS-Cx_R_strain.DEseq5 |
657 |
975 |
Table 3 Summary of whole genome expression
Sample name |
Total gene number |
Number of known genes |
Number of novel genes |
Total transcript number |
Known transcript number |
Novel transcript number |
R_strain_1 |
14597 |
14036 |
561 |
19094 |
11413 |
7681 |
R_strain_2 |
14507 |
13921 |
586 |
19460 |
11415 |
8045 |
R_strain_3 |
14592 |
14016 |
576 |
19733 |
11533 |
8200 |
S_strain_1 |
14603 |
14040 |
563 |
19381 |
11440 |
7941 |
S_strain_2 |
14551 |
13980 |
571 |
19671 |
11562 |
8109 |
S_strain_3 |
14568 |
13999 |
569 |
19679 |
11501 |
8178 |
Prediction of new transcripts
After comparing the clean reads to the Culex pipiens pallens genome, we used StringTie [35] software to perform transcript reintegration for each sample and then used cuffmerge and cuffcompare software (both are packages in Cufflinks [36]) 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.
[Due to technical limitations, Table 4 could not be displayed here. Please see the supplementary files section to access the table.]
Detection of SNPs and INDELs
After comparing clean reads to the Culex pipiens pallens genome, we used Genome Analysis Toolkit (GATK) [37] software to call each chromosome, identify single nucleotide polymorphisms (SNPs) and insertion and deletion (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 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 included 3 major categories: molecular function, cellular composition and biological process. Subsequently, the 3 categories were separately classified. The results are provided in Figure 4.
The obtained DEGs were compared with the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, and the related metabolic pathways were identified. The results are provided in Figure 5.
Validation of the 3 target cuticle protein genes
To ensure 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 monomodal, and the electrophoresis results revealed specific target fragments. The average Ct value for the actin gene in the sensitive group was 25.905; the average Ct value of the actin gene in the resistant group was 26.227; the average Ct value of the XM_001863852 gene in the sensitive group was 26.813; the average Ct value of the XM_001863852 gene in the resistant group was 29.633; the average Ct value of the XM_001845883.1 gene in the sensitive group was 36.797; the average Ct value of the XM_001845883.1 gene in the resistant group was 35.93; the average Ct value of the XM_001845881 gene in the sensitive group was 32.647; and the average Ct value of the XM_001845881 gene in the resistant group was 33.837 (see Tables 6 and 7 for details).
Subsequently, we used the 2-△△Ct method to analyse the expression of target gene mRNA in the extracted RNA and normalized the result based on the housekeeping gene.
The specific calculation method is as follows:
△△Ct = (Ct target gene - Ct housekeeping gene) experimental group - (Ct target gene - Ct housekeeping gene) control group
The relative expression level of the target gene = 2^-△△Ct, which indicates the fold change in expression for the target gene in the experimental group relative to the control group.
The calculation results are provided in Tables 6 and 7.
Table 6 Quantitative PCR results
Target ID |
Target Name |
Sample ID |
Sample Name |
Ct Avg (SDM) |
Rel. Qty (SDM) |
T001 |
Actin |
S001 |
S (sensitive) 1 |
25.475 |
1.00E+00 |
T001 |
Actin |
S002 |
S 2 |
26.205 |
1.00E+00 |
T001 |
Actin |
S003 |
S 3 |
26.035 |
1.00E+00 |
T001 |
Actin |
S004 |
R (resistance) 1 |
25.37 |
1.00E+00 |
T001 |
Actin |
S005 |
R 2 |
26.19 |
1.00E+00 |
T001 |
Actin |
S006 |
R 3 |
27.12 |
1.00E+00 |
T002 |
XM_001863852 |
S001 |
S 1 |
31.79 |
1.00E+00 |
T002 |
XM_001863852 |
S002 |
S 2 |
24.00 |
3.67E+02 |
T002 |
XM_001863852 |
S003 |
S 3 |
24.65 |
2.08E+02 |
T002 |
XM_001863852 |
S004 |
R 1 |
30.235 |
2.73E+00 |
T002 |
XM_001863852 |
S005 |
R 2 |
30.61 |
3.72E+00 |
T002 |
XM_001863852 |
S006 |
R 3 |
28.055 |
4.16E+01 |
T003 |
XM_001845883.1 |
S001 |
S 1 |
36.19 |
1.00E+00 |
T003 |
XM_001845883.1 |
S002 |
S 2 |
37.07 |
9.01E-01 |
T003 |
XM_001845883.1 |
S003 |
S 3 |
37.13 |
7.68E-01 |
T003 |
XM_001845883.1 |
S004 |
R 1 |
35.43 |
1.58E+00 |
T003 |
XM_001845883.1 |
S005 |
R 2 |
35.94 |
1.95E+00 |
T003 |
XM_001845883.1 |
S006 |
R 3 |
36.42 |
2.67E+00 |
T004 |
XM_001845881 |
S001 |
S 1 |
36.48 |
1.00E+00 |
T004 |
XM_001845881 |
S002 |
S 2 |
29.33 |
2.36E+02 |
T004 |
XM_001845881 |
S003 |
S 3 |
32.13 |
3.01E+01 |
T004 |
XM_001845881 |
S004 |
R 1 |
33.04 |
1.01E+01 |
T004 |
XM_001845881 |
S005 |
R 2 |
35.15 |
4.13E+00 |
T004 |
XM_001845881 |
S006 |
R3 |
33.32 |
2.80E+01 |
Table 7 Calculation of the relative quantitative Ct values for target genes and the internal reference gene
|
β-actin |
XM_001863852 |
XM_001845883.1 |
XM_001845881 |
Average Ct value in the sensitive group |
25.905 |
26.813 |
36.797 |
32.647 |
Average Ct value in the resistant group |
26.227 |
29.633 |
35.93 |
33.837 |
△ Ct value in the sensitive group |
----- |
0.908 |
10.892 |
6.742 |
△ Ct value in the resistant group |
----- |
3.406 |
9.703 |
7.61 |
△△ Ct value |
----- |
2.498 |
-1.189 |
0.868 |
2^-△△ Ct value |
----- |
0.177 |
2.281 |
0.548 |
It can be seen from Table 7 that 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, and the fold change in expression was 2.281.