全基因组重测序解析重庆武隆中华蜜蜂遗传变异及种群分化特征

任 勤, 姬聪慧, 龙小飞, 罗文华, 王瑞生, 陈 恒, 高丽娇, 曹 兰, 刘佳霖

(1.重庆市畜牧科学院,重庆 402460;2.重庆市武隆区畜牧技术推广站,重庆 408500)

中华蜜蜂(Apisceranacerana),简称“中蜂”,是我国最主要的本土蜜蜂,约有2 000多年的养殖历史,主要分布于南方山区和丘陵。中蜂不仅在维持野生植物多样性发展方面做出了巨大贡献,也为山区蜂农创造了较高的经济效益,助力推动乡村振兴和践行“绿水青山就是金山银山”的国家战略[1-3]。与西方蜜蜂(ApismelliferaLinnaeus)相比,中蜂具有独特的生物学特性,如较强的抗蜂螨能力、耐寒能力以及擅长采集零星蜜源[4-5]。在长期进化过程中,中蜂已适应我国山地和丘陵等环境,更适宜于为本土开花植物授粉,显著提高了我国农作物的产量和品质。然而,自1896年西方蜜蜂引入我国以来,由于资源的保护力度不足,中蜂的分布区域大幅缩减。2005年,中蜂种群数量已减少80%以上,导致我国野生植物和农作物得不到充分授粉,降低了中蜂的生态效益和经济效益[6]。2006年,中蜂归为畜禽遗传资源保护品种。2022年的研究显示,大部分地区中蜂的种群数量相对稳定,但西藏波密的种群数量仍持续下降[6]。因此,保护中蜂,特别是保护具有地方特色的优良中蜂种质资源,仍是我国蜂产业的重要任务之一。

武隆中蜂是以重庆市武隆区为核心的武陵山区自然形成的品种,且具有地方特色的优质中蜂种质资源。主要分布于重庆市东南部乌江下游,武陵山和大娄山峡谷地带。前期调研显示,武隆中蜂的良好抗逆性能和生产性能使其具有重要的研究价值和利用前景。然而,目前对于武隆中蜂的遗传变异和种群分化特征的研究较少,缺乏针对特色遗传变异位点研究导,限制了本地中蜂种质资源的保护和开发利用。

本研究将采集武隆中蜂进行全基因组重测序,并参考江西五府山中蜂基因组序列,检测单核苷酸多态性(single nucleotide polymorphism,SNP)和小片段的插入与缺失(small indel)位点,筛选出突变基因,以明确武隆中蜂的遗传变异特征,同时,结合已公布的我国9个地理型中蜂的全基因组重测序数据,基于SNP构建系统进化树,以期为绘制武隆中蜂的DNA指纹图谱提供基础,并为系统阐述我国中蜂种质资源的环境适应机制和遗传多样性提供依据。

1 材料与方法

1.1 样品采集与收集

2022年8—10月随机挑选武隆中蜂亲本蜂群和人工繁育的F1代蜂群各5群,每群采集100～200只工蜂保存于无水乙醇中。每群随机抽取3只工蜂,提取总DNA,并在4 ℃保存待测序。亲本蜂群为武隆本地野生中蜂种群,于2022年3—4月收集于重庆市武隆区鸭江乡,随后转移至武隆区凤山街道泡桐湾重庆市畜牧科学院中蜂良种选育场。2022年5月,通过人工育王的方式获得了F1代蜂王,并将其引入预先准备好的无王蜂群中,在交尾场内进行自然交尾。

1.2 全基因组重测序

委托北京百迈客生物科技有限公司利用Illumina/BGI测序平台对武隆中蜂样品进行全基因组重测序。样品包括亲本蜂群WL1(WL11、WL12、WL13、WL14和WL15)和人工繁育的F1代蜂群WL2(WL21、WL22、WL23、WL24和WL25)。测序数据已提交至GenBank,登录号为PRJNA973123。随后利用生物信息学分析平台BMKCloud进行后续的数据分析。

1.3 遗传变异分析

对原始数据进行处理(删除接头、去除ploy-N和低质量序列),获得了有效数据(clean data)。将clean data与参考基因组(江西五府山中蜂,GCA_002290385.1)进行比对。比对结果经Samtools软件[7]进行排序和去冗余处理,然后使用GATK软件[8]检测变异基因,检测结果过滤后获得高质量SNP和small indel数据,最终使用SNPsEff软件[9]进行SNP和small indel注释。基因功能通过序列比对,基于非冗余蛋白质序列数据库(non-redundant protein sequence database, Nr)、蛋白质直系同源簇(clusters of orthologous groups of proteins, COG)数据库、京都基因与基因组百科全书(Kyoto encyclopedia of genes and genomes, KEGG)数据库和基因本体(gene ontology, GO)数据库进行注释。

1.4 种群遗传分化分析

从GeneBank数据库中下载我国不同地区的57群中蜂的全基因组重测序原始数据。这57个群体分属于我国的12个地区和9个中蜂地理型。基于SNP数据并利用MEGA X软件[10]进行邻接法构建系统进化树,以西域黑蜂(Apismelliferasinisxinyuan)的全基因组重测序数据(GenBank登录号:PRJNA418015)作为外群。采用Admixture软件[11]进行群体结构分析,设定亚群数目(K值)为1～10个进行聚类,通过交叉验证确定最优的分群结果,并生成每个K值的堆叠图。此外,利用Eigensoft软件[12]进行主成分分析。最后,利用vcftools软件包[13]计算群体分化指数(F-statistics,Fst)。

表1 57群中华蜜蜂的全基因组重测序原始数据信息Table 1 Information of re-sequencing data for 57 A.cerana cerana colonies

2 结果与分析

2.1 数据质量分析

本试验共完成10个武隆中蜂样品的全基因组重测序(表2)。分析的clean data量达到33.53 Gb,Q30均超过94%。样品与参考基因组的平均比对率为91.02%,平均覆盖深度为12×,基因组覆盖度大于99.1%。这表明各试验组的数据完整性较好。

表2 样本重测序的数据质量分析Table 2 Quality statistics of samples for re-sequencing

2.2 SNP检测

与参考基因组(江西五府山中蜂,属华中中蜂)相比,各样品检测获得的SNP统计如表3所示,删除重复SNP,10个武隆中蜂样品共有3 886 321个高质量SNP,其中转换类型的SNP有1 302 218个,颠换类型的SNP有2 584 103个。

表3 样品与参考基因组比对检测得到的SNPTable 3 SNP of samples in comparison with reference genome

武隆中蜂样品中SNP的突变类型主要为C∶G>T∶A和T∶A>C∶G,占SNP总量的83.85%,其次是C∶G>G∶C和T∶A>A∶T,数量最少的突变类型为T∶A>G∶C和C∶G>A∶T。各样品的SNP主要集中在基因间区、内含子和基因上游区域,分别占SNP总量的70.83%、34.74%和11.62%。发生在编码区(coding sequence, CDS)的SNP占总量的7.40%。其中,79.01%为同义编码突变,20.58%为非同义编码突变。

2.3 Small indel检测

删除重复的small indel,本试验共获得1 392 028个small indel。其中,基因组插入类型占全部small indel的40.40%,基因组删除类型占59.57%;在CDS内的插入类型占small indel总量的0.24%,删除类型占0.31%(表4)。

表4 样品与参考基因组比对检测得到的small indelTable 4 Small indel of samples in comparison with reference genome

各样品的small indel在基因组上主要集中于基因间区、内含子和基因上游区域,分别占small indel总量的57.94%、27.08%和7.58%。发生在CDS上的small indel占总量的0.69%。其中,47.46%为移码突变,23.56%为密码子删除,14.43%为密码子改变+密码子删除。

2.4 突变基因分析

与参考基因组相比,样品中各种变异产生的突变基因统计信息如表5所示。删除重复基因后,武隆中蜂样品中共有9 579个突变基因。其中,发生非同义SNP突变的基因8 240个,CDS区发生small indel的基因有2 375个。

表5 各种变异产生的突变基因的分类统计Table 5 Classification and statistics of mutant genes produced by various variations

进一步筛选武隆中蜂样品中具有相同变异位点的突变基因,最终获得发生非同义SNP突变的基因589个,CDS区发生small indel的基因214个。根据GO分类和KEGG富集分析的结果显示,与参考基因组相比,武隆中蜂具有相同变异位点且发生非同义SNP突变的基因注释到生物过程、细胞组分和分子功能的分别有980、809和551个(图1A),这些基因主要富集在赖氨酸降解的KEGG通路(图1B)。而CDS区发生small indel的基因在生物过程、细胞组分和分子功能方面分别有299、311和176个(图2A),这些基因主要富集在轴突导向、细胞外基质受体互作和MAPK信号通路(图2B)。

图1 发生非同义SNP突变的基因的GO分类(A)和KEGG富集分析(B)Fig.1 GO classification and KEGG enrichment of genes with non-synonymous SNP

图2 CDS区发生small indel的基因的GO分类(A)和KEGG富集分析(B)Fig.2 GO classification and KEGG enrichment of genes with small indel in coding sequence

在已有数据中筛选出具有相同SNP和small indel的突变基因16个(表6)。其中包括2个嗅觉受体基因(gene-APICC_02321和gene-APICC_02315)、1个细胞色素P450基因(gene-APICC_00983)以及1个内质网氨肽酶基因(gene-APICC_07171)。

2.5 种群遗传分化

基于武隆和我国12个地区中蜂种群的SNP构建的系统进化树显示,长白山中蜂、阿坝中蜂、海南中蜂和华南中蜂(福建漳州)最早从进化树中分离出来,随后分为两个大的分支。其中,一个分支由武隆中蜂、北方中蜂和华中中蜂组成,在这个分支中,武隆中蜂与其他种群(北方中蜂和华中中蜂)形成姐妹支;另一支由云贵高原中蜂、华南中蜂(广西河池)、西藏中蜂和滇南中蜂组成(图3)。

图3 基于67群中华蜜蜂的SNP数据构建系统进化树(1 000次重复)Fig.3 Phylogenetic tree of 67 A.cerana cerana colonies by neighbor-joining method based on SNP(1 000 replicates)

遗传结构分析(图4A)显示,假设存在5～10个样品(K=5～10)时,可清楚区分西藏中蜂、长白山中蜂、阿坝中蜂和海南中蜂。在K=10时,武隆中蜂也无法被完全区分出来。主成分分析显示,来自我国不同地区的67群中蜂被划分为阿坝中蜂、海南中蜂、长白山中蜂、西藏中蜂和其他(包括武隆中蜂、华南中蜂、北方中蜂、华中中蜂、云贵高原中蜂和滇南中蜂)。前两个主成分共解释了整体变异量的10.05%(图4B)。剔除阿坝中蜂、海南中蜂、长白山中蜂和西藏中蜂,主成分分析显示其他中蜂种群可划分为武隆中蜂和其他(包括华南中蜂、北方中蜂、华中中蜂、云贵高原中蜂和滇南中蜂)(图4C),PC1和PC2分别解释了整体变异量的3.75%和2.57%。

A中不同颜色代表样品在划分为K个亚群时的遗传结构不同。

Fst分析结果显示(图5),全体两两比较中,Fst最大值(0.424 8)出现于长白山中蜂和海南中蜂群体之间,最小值(0.024 1)出现于武隆中蜂亲本蜂群WL1与人工繁育的F1代蜂群WL2群体之间。武隆中蜂WL1和WL2与长白山中蜂、海南中蜂、阿坝中蜂和西藏中蜂群体间的Fst最大,分别为0.231 4、0.191 6、0.158 2、0.138 6和0.231 8、0.191 7、0.156 6、0.141 5,与云贵高原中蜂和华中中蜂群体间的Fst最小,分别为0.029 0、0.028 0和0.028 7、0.029 9。其次,武隆中蜂与北方中蜂(陕西安康)、华南中蜂(广西河池)维持较小的Fst。

图5 武隆中华蜜蜂与其他12个地区中华蜜蜂的Fst热图Fig.5 Heatmap of Fst of Wulong A.cerana cerana with 12 other geographic populations

3 讨论

近年来,全基因组重测序技术已广泛应用于蜜蜂种质资源的遗传多样性研究,在挖掘特色分子标记和重要功能基因、解析地方种群适应性进化的分子机制等方面发挥着重要作用[14-16]。本研究采集武隆中蜂进行全基因组重测序,分析量clean data达33.53 Gb,以同属华中中蜂的江西五府山中蜂基因组为参考,筛选出了3 886 321个高质量的SNP和1 392 028个small indel。其中,武隆中蜂具有相同变异位点且发生非同义SNP突变的基因有589个,CDS区发生small indel的基因214个。另外,还发现16个具有相同SNP和small indel的变异基因,表明这些突变基因可能在武隆中蜂适应本地自然环境的过程中发挥着重要作用。

武隆中蜂具有相同变异位点且发生非同义SNP突变的基因主要富集在赖氨酸降解的KEGG通路,CDS区发生small indel的基因富集在轴突引导、细胞外基质受体互作和MAPK信号通路。赖氨酸在蜜蜂生长发育、脂肪代谢和免疫力调控等方面发挥着重要作用[17-18]。在机体内,赖氨酸的降解主要通过酵母氨酸途径和哌啶酸途径转换为乙酰辅酶A,进而在TCA循环中参与营养物质代谢,为机体供能并增强免疫力[19-20]。轴突导向是神经元发育中的基本过程,在神经网络的形成中发挥着重要作用[21-22],细胞外基质受体互作与昆虫脑部的发育相关[23],轴突导向和细胞外基质受体互作通路的变化可能涉及神经系统的进化,进而调节昆虫的运动、飞行和趋避等行为。MAPK信号通路是最重要的信号转导系统之一,介导细胞增殖和机体生长发育等多种生理过程[24]。近期研究表明,蜂王长寿可能与MAPK信号通路增强细胞活性相关[25],此外,该通路还参与蜜蜂对病原微生物和低温胁迫的响应[26-27],说明武隆中蜂优良的抗逆性能可能与MAPK信号通路相关基因的突变有关。

发达的嗅觉系统是昆虫适应环境的重要体现,它介导蜜蜂的信息交流、定位和采集等行为,是中蜂擅长采集零星蜜源的基础[28-29]。嗅觉受体蛋白参与昆虫的嗅觉识别过程,它们与外界气味分子和气味结合蛋白的复合体结合,将化学信号转换为电信号并传递到大脑中枢神经系统[28,30]。细胞色素P450是蜜蜂解毒代谢外源物质的关键酶,能在温和条件下催化多种氧化反应,改变外源物质的结构,降低其毒性,并使其从亲脂性转变为亲水性[31-32]。研究表明,P450参与蜜蜂对植物次生物质和多种农药的解毒代谢[33-34]。此外,细胞色素P450基因(CYP6a17)还参与调控昆虫对温度的偏好,这种偏好是果蝇(Drosophilamelanogaster)等变温动物优化体温的重要行为策略[35]。本试验筛选出了16个在武隆中蜂中具有相同SNP和small indel的突变基因,其中包括2个嗅觉受体基因和1个细胞色素P450基因(CYP6a17),这暗示嗅觉进化、解毒代谢和对温度偏好性的转变可能是武隆中蜂适应重庆本地自然环境的重要机制。

目前,我国中蜂分为9个地理型,包括北方中蜂、华南中蜂、华中中蜂、云贵高原中蜂、长白山中蜂、海南中蜂、滇南中蜂、阿坝中蜂和西藏中蜂[36]。由于长期的地理隔离,海南中蜂、阿坝中蜂、西藏中蜂和长白山中蜂与其他地理型中蜂的亲缘关系较远[15,37-38]。本研究通过对系统进化树、遗传结构和遗传分化的分析,成功区分了这4种中蜂与其他样本,说明长期的地理隔离是推动中蜂种群分化的关键,能有效防止外来遗传资源的干扰和污染。此外,阿坝中蜂、西藏中蜂和长白山中蜂长期生活在高原地区,这些地区的生境、地形复杂,植物多样,为中蜂的繁衍和发展提供了重要条件,也促进了中蜂逐渐适应高原生境。唐相友等[14]的研究表明,在我国的一些高原地区,中蜂种群存在较高的遗传分化水平,这暗示在青藏高原等地区可能存在更多具有独特地理遗传特性的中蜂群体。

重庆是我国中蜂的主要养殖区。本研究结果显示,在系统进化树和遗传结构的主成分分析中,武隆中蜂与我国其他地理型中蜂有显著的遗传分化,说明武隆中蜂是重庆地区独特的中蜂种质遗传资源。根据第3次全国畜禽遗传资源普查,重庆中蜂属于华中中蜂。本研究中,武隆中蜂与湖北神农架中蜂群体的遗传分化水平最低,且在进化树中归为同一支,说明两个同属华中中蜂的种群亲缘关系最近。其次,武隆中蜂与陕西安康的北方中蜂在进化树和Fst分析中也显示出较近的亲缘关系。尽管武隆中蜂与云贵高原中蜂和广西河池的华南中蜂的遗传分化水平最低,但在进化树中属于不同的大支,因此,推测武隆中蜂与云贵高原中蜂和华南中蜂的亲缘关系远于华中中蜂和北方中蜂。武隆中蜂亲本蜂群与人工繁育的F1代蜂群群体之间的Fst最小,且两个群体在进化树中归为同一分支,表明初步的人工繁育并未对武隆中蜂的遗传多样性产生影响。在群体遗传结构分析中,在K=10时,仅能将西藏中蜂、长白山中蜂、海南中蜂和阿坝中蜂区分出来,暗示其他地区中蜂可能存在自然或人为干预的遗传资源交流。多年来,重庆本地蜂农为提高养殖效益,随意从其他地区引进中蜂种群,污染了本地遗传资源,导致品种杂化和生产性能快速下降。

综上所述,本研究揭示了重庆武隆中蜂的遗传变异和种群分化特征,并推测了与适应性进化相关的功能通路和基因,为进一步挖掘、保护和开发利用重庆本地中蜂种质资源提供了重要的理论依据。同时,在引入外来中蜂种群时要注意避免对本地遗传资源的干扰和污染,以保护和提升中蜂的生产性能和适应性。