在重庆、云南、安徽3个省(市)的稻田中采集4龄按蚊幼虫^[27]，并在当地养蚊室内饲养至成虫。室内饲养条件为27 ℃、相对湿度80%。随后，提取蚊虫足部进行物种的分子鉴定^[28]。雌性成蚊羽化72 h后，采用世界卫生组织(WHO)推荐的药膜生物测定法进行抗性检测，将雌性成蚊暴露在含有0.05%溴氰菊酯药膜接触桶中60 min，随后将其转移到恢复桶中，若蚊虫在60 min内死亡或者机械刺激无反应，则判定为拟除虫菊酯敏感个体；24 h后统计蚊虫的死亡率，存活个体判定为拟除虫菊酯抗性个体。敏感品系成蚊的相对敏感基线LC₅₀(半致死浓度)为0.006 7 mg/L，云南抗性品系(YN-FR)、重庆抗性品系(CQ-FR)、安徽抗性品系(AH-FR)的抗性倍数分别为500倍、800倍、800倍，均达到国家标准的高抗水平^[29]。最终，从重庆、云南、安徽3个省(市)分别收集6只拟除虫菊酯抗性成蚊(Field Pyrethroid-Resistant Strain，FR)和6只拟除虫菊酯敏感成蚊(Field Pyrethroid-Susceptible Strain，FS)，共36个样本用于后续分析。

将36个中华按蚊野外样本送往华大公司进行测序，并采用Illumina HiSeq2000平台测序获得重测序数据，中华按蚊的参考基因组(GCA-000441895.2)由重庆师范大学昆虫与分子生物学研究所提供^[30-31]。

使用fastp软件^[32]对原始测序数据(Raw reads)进行质量控制与过滤，获得高质量测序读段(Clean reads)。然后，使用BWA软件^[33]将36个中华按蚊样本的高质量测序读段与参考基因组进行比对。通过SAMtools(处理高通量测序比对数据的软件集)^[34]将比对生成的SAM(存储测序序列与参考基因组比对结果的文本格式)文件转换为排序后的BAM(SAM的二进制压缩格式)文件，并去除PCR(聚合酶链式反应)重复序列，最终获得用于后续分析的比对结果文件。

使用GATK^[35]对36个中华按蚊样本的比对结果进行SNPs和InDels变异检测。为确保数据质量，对原始变异位点进行过滤，并使用ANNOVAR软件^[36]进行注释分析。

前期基于中华按蚊基因组注释蛋白库，使用BLASTP和Hmmsearch两种方法，获得中华按蚊P450基因^[17]。再使用VCFtools软件^[37]对GATK过滤后的变异结果进一步筛选，剔除次要等位基因频率(maf)＜0.05的位点、去除基因组型缺失率超过20%的位点，同时仅保留二等位SNPs变异位点。基于中华按蚊基因组注释文件，使用Python脚本提取P450基因家族成员的注释信息。最后，根据P450基因在基因组中的起始和终止位置，使用VCFtools(参数：-bed)提取All.filtered.vcf文件中的变异位点^[38]，作为P450基因家族成员的SNPs数据，用于后续群体结构分析。

使用PHYLIP软件中的邻近法(Neighbor-Joining Methods)构建系统发育树(F84模型，100次Bootstrap重复)，并使用Ggtree进行可视化展示。

基于36个中华按蚊样本的遗传变异数据，使用GCTA软件^[39]进行主成分分析，并根据主成分聚类结果将36个中华按蚊样本划分为不同的遗传亚群。

使用Admixture软件^[40]进行群体遗传结构分析。通过最大似然估计法确定最优群体数目K(预设范围2-N)，并使用贝叶斯算法计算每个K值下的群体分群。

36个样本分为6个组，即安徽抗性组AH-FR、安徽敏感组AH-FS、重庆抗性组CQ-FR、重庆敏感组CQ-FS、云南抗性组YN-FR、云南敏感组YN-FS。每组6个样本，且相同地区的抗性组和敏感组比较，如安徽抗性组AH-FR与安徽敏感组AH-FS比较，另外两个地区也是如此。将组内基因型一致，但组间基因型存在差异的位点定为差异位点。筛选差异位点的标准为：在安徽抗性组(AH-FR)中，6个样本中有80%的样本(即5个样本)具有一致的基因型(记为该组的优势基因型A)；同时在安徽敏感组(AH-FS)中，6个样本中也要有80%的样本(即5个样本)具有一致的基因型(记为该组的优势基因型B)。仅当A与B不同时，该位点才被定义为差异位点。然后，对这些位点进行注释，确定其所在基因位置及突变类型。

序列比对结果如表 1所示。由表 1可知，36个样本的比对率为78.26%~98.64%，平均测序深度为8.769×~27.817×，说明样本建库良好；参考基因组的覆盖度为89.73%~97.31%，说明比对情况良好，可用于下一步群体遗传结构分析。

P450基因的变异数量统计结果如表 2所示。36个中华按蚊样本的P450基因家族共检测到59 070个变异位点，其中SNPs数量为46 952个，平均每kb区间内鉴定出0.68个SNPs位点；InDels数量为12 118个，平均每kb区间内鉴定出0.19个InDels位点。

P450基因SNPs和InDels的注释统计结果如表 3所示。SNPs被注释到基因间区(intergenic regions)的数量为62 011个，内含子区域数量为45 750个，下游区域SNPs数量为11 416个，上游区域SNPs数量为11 795个，外显子区域SNPs数量为13 971个，5'端非翻译区(5'UTR)SNPs数量为1 085个，3'端非翻译区(3'UTR)SNPs数量为1 478个，剪切位点区域(splicing)SNPs数量为25个。InDels被注释到基因间区(intergenic regions)的数量为9 238个，内含子区域InDels数量为7 487个，上游区域InDels数量为1 976个，3'端非翻译区(3'UTR)InDels数量为248个，外显子区域InDels数量为224个，下游区域InDels数量为1 831个，5'端非翻译区(5'UTR)InDels数量为152个，剪切位点区域(splicing)InDels数量为13个。由此可见，中华按蚊P450基因的SNPs和InDels主要分布在基因间区(intergenic regions)，其次是内含子区。

基于P450-SNP的中华按蚊系统进化树分析结果如图 1所示。36个中华按蚊样本按照地理种群聚类为3组，种群间的遗传差异和地理区域相对应。安徽12个样本中AH-FS-13、AH-FR-2自展值较低(＜70%)；而其他样本以较高的自展值聚类(＞80%)，说明这些样本可能来自共同祖先。云南12个样本中YN-FR-5、YN-FR-4、YN-FS-4、YN-FR-3自展值较低(＜70%)；而其他样本都以较高的自展值聚类(＞80%)，说明这些样本可能来自共同祖先。重庆12个样本中CQ-FR-3、CQ-FS-2、CQ-FR-6、CQ-FS-1、CQ-FR-5、CQ-FS-4自展值较低(＜70%)；而其他样本都以较高的自展值聚类(＞80%)，说明这些样本可能来自共同祖先。CQ-FS-5与CQ-FS-3以自展值100%聚类在一起，说明它们可能来自共同祖先。CQ-FR-2和CQ-FR-3与重庆种群其他样本显著分离，表明重庆种群的群内遗传差异较大，含有不同的祖先来源。

主成分分析结果如图 2所示。36个中华按蚊样本聚类与地理分布基本一致，其中云南(YN)和安徽(AH)种群各自形成独立聚类，而重庆(CQ)种群内部个体分散，表明其遗传多样性较高。此外，图 2中PC1、PC2和PC3分别代表第一、第二和第三主成分，其后方括号内的百分比(16.422%、5.991%、4.661%)为各主成分的方差贡献率，由eigenval.xls特征值文件计算得出。由此可知，主成分分析结果与系统进化树的分析结果一致。

群体遗传结构分析结果如图 3所示。随着设定的祖代群体数K从2逐步提高至4，中华按蚊36个样本的群体结构呈现相应的动态变化。当K=2时，36个中华按蚊可划分为两个群体；当K=3时，重庆种群个体从云南和安徽种群中分化出来，36个中华按蚊样本可划分为3个群体；当K=4时，重庆种群的群内差异较大，产生了一个亚群。综上，群体遗传结构分析结果与系统进化树分析和主成分分析结果一致。

中华按蚊拟除虫菊酯抗性相关的P450基因在重庆、云南、安徽抗性种群中的差异SNPs分别如表 4、表 5、表 6所示。结合实验室前期研究发现的3个(CYP9K1、CYP6P3v1、CYP9J10) P450基因被RNAi干扰后，降低了中华按蚊对拟除虫菊酯类杀虫剂的抗性。基于此，本研究重点分析这3个基因的遗传变异。CYP9K1基因在重庆、云南、安徽抗性种群中检测到21个差异SNPs位于基因间区，1个同义替换(L493L)位于外显子区域。CYP6P3v1基因在云南抗性种群中检测到1个同义替换(L118L)和1个差异SNPs，分别位于外显子区域和内含子区域。CYP9J10基因在重庆抗性种群中检测到1个非同义替换(S162N)，位于外显子区域。综上，CYP6P3v1存在外显子区域同义替换(L118L)和内含子区域SNPs；CYP9J10存在外显子区域非同义替换(S162N)；而CYP9K1主要存在基因间区SNPs。此外，这些基因均属于CYP6和CYP9家族成员。

本研究通过全基因组重测序分析发现，中华按蚊P450基因的SNPs和InDels主要分布于基因间区，这可能与基因间区在基因组中所占比例较高有关。通过进化树分析、主成分分析、群体遗传结构分析，发现36个中华按蚊的聚类关系与其地理分布基本一致。种群间的遗传差异与地理区域相对应，这与以往的研究结果一致，即中华按蚊的遗传距离和地理距离呈正相关，可能存在地理隔离现象^[41]。

本研究首次发现，CYP6P3v1存在外显子区域同义替换(L118L)和内含子区域SNPs，CYP9J10存在外显子区域非同义替换(S162N)，而CYP9K1主要存在基因间区SNPs。此外，这些基因均属于CYP6和CYP9家族成员，该结果与埃及伊蚊差异SNPs检测结果一致^[26]。

Dusfour等^[42]报道埃及伊蚊(Ae.aegypti) CYP9J10、CYP9J9、CYP6BB2、CYP6M11和CYP6N12基因的表达与拟除虫菊酯类杀虫剂抗性相关。Yan等^[17]通过RNA-seq转录组测序和RT-qPCR验证发现，中华按蚊CYP9K1基因在重庆和安徽2个地理种群的拟除虫菊酯抗性品系中上调表达；中华按蚊CYP6P3v1基因在重庆、云南、安徽3个地理种群的拟除虫菊酯抗性品系中均呈现上调表达。已有报道催命按蚊CYP9K1与拟除虫菊酯类杀虫剂抗性相关^[43]。中华按蚊的CYP6P3v1已被证实与拟除虫菊酯类杀虫剂抗性有关^[18]。

以往的研究报道，P450基因的代谢抗性主要是由基因上调表达和解毒酶活性增加所致^[44-45]。蚊虫对拟除虫菊酯的抗性是由解毒酶基因编码区变异产生的^[46]。昆虫P450基因编码区的遗传变异可导致杀虫剂抗性，例如埃及伊蚊(Ae.aegypti)抗性品系中的CYP6、CYP9和CYP12基因家族均存在与抗性相关的突变位点^[47]。非翻译区的突变能调控基因表达，进而导致杀虫剂抗性^[48]。研究表明，同义替换SNP可通过影响mRNA稳定性、翻译效率或剪接过程参与基因表达调控^[49]。内含子区域的突变可能改变内含子的剪切形式，产生不同的剪切体，由于这些剪切体稳定性不同，从而影响基因表达水平^[23]。Amichot等^[50]在黑腹果蝇(Drosophila melanogaster)的双对氯苯基三氯乙烷(DDT)中发现3个非同义替换(R335S、L336V、V476L)位于CYP6A2活性位点附近，这些突变位点与DDT的代谢抗性有关。本研究发现重庆和安徽抗性种群的CYP9K1仅存在一种变异类型(即基因间区SNPs)；云南抗性种群的CYP6P3v1存在两种变异类型(即外显子区域同义替换、内含子区域SNPs)；重庆抗性种群的CYP9J10存在一种变异类型(即外显子区域非同义替换)。

本研究外显子区域同义替换(L118L)、内含子区域及基因间区SNPs可能正向调控CYP6P3v1和CYP9K1的表达，进而介导拟除虫菊酯抗性；而外显子区域非同义替换(S162N)可能通过修饰CYP9J10基因编码的蛋白结构，增强其解毒酶活性，进而介导拟除虫菊酯抗性。以上这些P450基因的遗传变异与拟除虫菊酯杀虫剂抗性的关联性，还需进一步通过荧光素酶基因报告实验和解毒酶活性代谢实验来验证SNPs位点是否参与基因表达调控或者介导解毒酶活性。此外，本研究可能遗漏低频变异，后续需进一步扩大野外样本采集数量，以验证遗传变异的可靠性和发生频率。

本研究综合进化树分析、主成分分析、群体遗传结构分析结果，发现36个中华按蚊样本的聚类结果与地理分布一致，不同种群间可能存在地理隔离现象。CYP9K1、CYP6P3v1、CYP9J10的遗传变异可能介导拟除虫菊酯抗性。上述P450基因遗传变异与拟除虫菊酯抗性的关联，仍需通过荧光素酶报告基因实验和解毒酶活性代谢实验进一步验证，以明确相关SNPs能否调控基因表达或影响解毒酶活性。本研究结果为深入解析中华按蚊对拟除虫菊酯的抗性机制，特别是P450基因的表达调控与解毒酶活性提供了关键的数据基础。