任 勤, 姬聰慧, 龍小飛, 羅文華, 王瑞生, 陳 恒, 高麗嬌, 曹 蘭, 劉佳霖

(1.重慶市畜牧科學院,重慶 402460;2.重慶市武隆區(qū)畜牧技術推廣站,重慶 408500)

中華蜜蜂(Apisceranacerana),簡稱“中蜂”,是我國最主要的本土蜜蜂,約有2 000多年的養(yǎng)殖歷史,主要分布于南方山區(qū)和丘陵。中蜂不僅在維持野生植物多樣性發(fā)展方面做出了巨大貢獻,也為山區(qū)蜂農(nóng)創(chuàng)造了較高的經(jīng)濟效益,助力推動鄉(xiāng)村振興和踐行“綠水青山就是金山銀山”的國家戰(zhàn)略[1-3]。與西方蜜蜂(ApismelliferaLinnaeus)相比,中蜂具有獨特的生物學特性,如較強的抗蜂螨能力、耐寒能力以及擅長采集零星蜜源[4-5]。在長期進化過程中,中蜂已適應我國山地和丘陵等環(huán)境,更適宜于為本土開花植物授粉,顯著提高了我國農(nóng)作物的產(chǎn)量和品質。然而,自1896年西方蜜蜂引入我國以來,由于資源的保護力度不足,中蜂的分布區(qū)域大幅縮減。2005年,中蜂種群數(shù)量已減少80%以上,導致我國野生植物和農(nóng)作物得不到充分授粉,降低了中蜂的生態(tài)效益和經(jīng)濟效益[6]。2006年,中蜂歸為畜禽遺傳資源保護品種。2022年的研究顯示,大部分地區(qū)中蜂的種群數(shù)量相對穩(wěn)定,但西藏波密的種群數(shù)量仍持續(xù)下降[6]。因此,保護中蜂,特別是保護具有地方特色的優(yōu)良中蜂種質資源,仍是我國蜂產(chǎn)業(yè)的重要任務之一。

武隆中蜂是以重慶市武隆區(qū)為核心的武陵山區(qū)自然形成的品種,且具有地方特色的優(yōu)質中蜂種質資源。主要分布于重慶市東南部烏江下游,武陵山和大婁山峽谷地帶。前期調研顯示,武隆中蜂的良好抗逆性能和生產(chǎn)性能使其具有重要的研究價值和利用前景。然而,目前對于武隆中蜂的遺傳變異和種群分化特征的研究較少,缺乏針對特色遺傳變異位點研究導,限制了本地中蜂種質資源的保護和開發(fā)利用。

本研究將采集武隆中蜂進行全基因組重測序,并參考江西五府山中蜂基因組序列,檢測單核苷酸多態(tài)性(single nucleotide polymorphism,SNP)和小片段的插入與缺失(small indel)位點,篩選出突變基因,以明確武隆中蜂的遺傳變異特征,同時,結合已公布的我國9個地理型中蜂的全基因組重測序數(shù)據(jù),基于SNP構建系統(tǒng)進化樹,以期為繪制武隆中蜂的DNA指紋圖譜提供基礎,并為系統(tǒng)闡述我國中蜂種質資源的環(huán)境適應機制和遺傳多樣性提供依據(jù)。

1 材料與方法

1.1 樣品采集與收集

2022年8—10月隨機挑選武隆中蜂親本蜂群和人工繁育的F1代蜂群各5群,每群采集100～200只工蜂保存于無水乙醇中。每群隨機抽取3只工蜂,提取總DNA,并在4 ℃保存待測序。親本蜂群為武隆本地野生中蜂種群,于2022年3—4月收集于重慶市武隆區(qū)鴨江鄉(xiāng),隨后轉移至武隆區(qū)鳳山街道泡桐灣重慶市畜牧科學院中蜂良種選育場。2022年5月,通過人工育王的方式獲得了F1代蜂王,并將其引入預先準備好的無王蜂群中,在交尾場內進行自然交尾。

1.2 全基因組重測序

委托北京百邁客生物科技有限公司利用Illumina/BGI測序平臺對武隆中蜂樣品進行全基因組重測序。樣品包括親本蜂群WL1(WL11、WL12、WL13、WL14和WL15)和人工繁育的F1代蜂群WL2(WL21、WL22、WL23、WL24和WL25)。測序數(shù)據(jù)已提交至GenBank,登錄號為PRJNA973123。隨后利用生物信息學分析平臺BMKCloud進行后續(xù)的數(shù)據(jù)分析。

1.3 遺傳變異分析

對原始數(shù)據(jù)進行處理(刪除接頭、去除ploy-N和低質量序列),獲得了有效數(shù)據(jù)(clean data)。將clean data與參考基因組(江西五府山中蜂,GCA_002290385.1)進行比對。比對結果經(jīng)Samtools軟件[7]進行排序和去冗余處理,然后使用GATK軟件[8]檢測變異基因,檢測結果過濾后獲得高質量SNP和small indel數(shù)據(jù),最終使用SNPsEff軟件[9]進行SNP和small indel注釋?；蚬δ芡ㄟ^序列比對,基于非冗余蛋白質序列數(shù)據(jù)庫(non-redundant protein sequence database, Nr)、蛋白質直系同源簇(clusters of orthologous groups of proteins, COG)數(shù)據(jù)庫、京都基因與基因組百科全書(Kyoto encyclopedia of genes and genomes, KEGG)數(shù)據(jù)庫和基因本體(gene ontology, GO)數(shù)據(jù)庫進行注釋。

1.4 種群遺傳分化分析

從GeneBank數(shù)據(jù)庫中下載我國不同地區(qū)的57群中蜂的全基因組重測序原始數(shù)據(jù)。這57個群體分屬于我國的12個地區(qū)和9個中蜂地理型?；赟NP數(shù)據(jù)并利用MEGA X軟件[10]進行鄰接法構建系統(tǒng)進化樹,以西域黑蜂(Apismelliferasinisxinyuan)的全基因組重測序數(shù)據(jù)(GenBank登錄號:PRJNA418015)作為外群。采用Admixture軟件[11]進行群體結構分析,設定亞群數(shù)目(K值)為1～10個進行聚類,通過交叉驗證確定最優(yōu)的分群結果,并生成每個K值的堆疊圖。此外,利用Eigensoft軟件[12]進行主成分分析。最后,利用vcftools軟件包[13]計算群體分化指數(shù)(F-statistics,Fst)。

表1 57群中華蜜蜂的全基因組重測序原始數(shù)據(jù)信息Table 1 Information of re-sequencing data for 57 A.cerana cerana colonies

2 結果與分析

2.1 數(shù)據(jù)質量分析

本試驗共完成10個武隆中蜂樣品的全基因組重測序(表2)。分析的clean data量達到33.53 Gb,Q30均超過94%。樣品與參考基因組的平均比對率為91.02%,平均覆蓋深度為12×,基因組覆蓋度大于99.1%。這表明各試驗組的數(shù)據(jù)完整性較好。

表2 樣本重測序的數(shù)據(jù)質量分析Table 2 Quality statistics of samples for re-sequencing

2.2 SNP檢測

與參考基因組(江西五府山中蜂,屬華中中蜂)相比,各樣品檢測獲得的SNP統(tǒng)計如表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,數(shù)量最少的突變類型為T∶A>G∶C和C∶G>A∶T。各樣品的SNP主要集中在基因間區(qū)、內含子和基因上游區(qū)域,分別占SNP總量的70.83%、34.74%和11.62%。發(fā)生在編碼區(qū)(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在基因組上主要集中于基因間區(qū)、內含子和基因上游區(qū)域,分別占small indel總量的57.94%、27.08%和7.58%。發(fā)生在CDS上的small indel占總量的0.69%。其中,47.46%為移碼突變,23.56%為密碼子刪除,14.43%為密碼子改變+密碼子刪除。

2.4 突變基因分析

與參考基因組相比,樣品中各種變異產(chǎn)生的突變基因統(tǒng)計信息如表5所示。刪除重復基因后,武隆中蜂樣品中共有9 579個突變基因。其中,發(fā)生非同義SNP突變的基因8 240個,CDS區(qū)發(fā)生small indel的基因有2 375個。

表5 各種變異產(chǎn)生的突變基因的分類統(tǒng)計Table 5 Classification and statistics of mutant genes produced by various variations

進一步篩選武隆中蜂樣品中具有相同變異位點的突變基因,最終獲得發(fā)生非同義SNP突變的基因589個,CDS區(qū)發(fā)生small indel的基因214個。根據(jù)GO分類和KEGG富集分析的結果顯示,與參考基因組相比,武隆中蜂具有相同變異位點且發(fā)生非同義SNP突變的基因注釋到生物過程、細胞組分和分子功能的分別有980、809和551個(圖1A),這些基因主要富集在賴氨酸降解的KEGG通路(圖1B)。而CDS區(qū)發(fā)生small indel的基因在生物過程、細胞組分和分子功能方面分別有299、311和176個(圖2A),這些基因主要富集在軸突導向、細胞外基質受體互作和MAPK信號通路(圖2B)。

圖1 發(fā)生非同義SNP突變的基因的GO分類(A)和KEGG富集分析(B)Fig.1 GO classification and KEGG enrichment of genes with non-synonymous SNP

圖2 CDS區(qū)發(fā)生small indel的基因的GO分類(A)和KEGG富集分析(B)Fig.2 GO classification and KEGG enrichment of genes with small indel in coding sequence

在已有數(shù)據(jù)中篩選出具有相同SNP和small indel的突變基因16個(表6)。其中包括2個嗅覺受體基因(gene-APICC_02321和gene-APICC_02315)、1個細胞色素P450基因(gene-APICC_00983)以及1個內質網(wǎng)氨肽酶基因(gene-APICC_07171)。

2.5 種群遺傳分化

基于武隆和我國12個地區(qū)中蜂種群的SNP構建的系統(tǒng)進化樹顯示,長白山中蜂、阿壩中蜂、海南中蜂和華南中蜂(福建漳州)最早從進化樹中分離出來,隨后分為兩個大的分支。其中,一個分支由武隆中蜂、北方中蜂和華中中蜂組成,在這個分支中,武隆中蜂與其他種群(北方中蜂和華中中蜂)形成姐妹支;另一支由云貴高原中蜂、華南中蜂(廣西河池)、西藏中蜂和滇南中蜂組成(圖3)。

圖3 基于67群中華蜜蜂的SNP數(shù)據(jù)構建系統(tǒng)進化樹(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)時,可清楚區(qū)分西藏中蜂、長白山中蜂、阿壩中蜂和海南中蜂。在K=10時,武隆中蜂也無法被完全區(qū)分出來。主成分分析顯示,來自我國不同地區(qū)的67群中蜂被劃分為阿壩中蜂、海南中蜂、長白山中蜂、西藏中蜂和其他(包括武隆中蜂、華南中蜂、北方中蜂、華中中蜂、云貴高原中蜂和滇南中蜂)。前兩個主成分共解釋了整體變異量的10.05%(圖4B)。剔除阿壩中蜂、海南中蜂、長白山中蜂和西藏中蜂,主成分分析顯示其他中蜂種群可劃分為武隆中蜂和其他(包括華南中蜂、北方中蜂、華中中蜂、云貴高原中蜂和滇南中蜂)(圖4C),PC1和PC2分別解釋了整體變異量的3.75%和2.57%。

A中不同顏色代表樣品在劃分為K個亞群時的遺傳結構不同。

Fst分析結果顯示(圖5),全體兩兩比較中,Fst最大值(0.424 8)出現(xiàn)于長白山中蜂和海南中蜂群體之間,最小值(0.024 1)出現(xiàn)于武隆中蜂親本蜂群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個地區(qū)中華蜜蜂的Fst熱圖Fig.5 Heatmap of Fst of Wulong A.cerana cerana with 12 other geographic populations

3 討論

近年來,全基因組重測序技術已廣泛應用于蜜蜂種質資源的遺傳多樣性研究,在挖掘特色分子標記和重要功能基因、解析地方種群適應性進化的分子機制等方面發(fā)揮著重要作用[14-16]。本研究采集武隆中蜂進行全基因組重測序,分析量clean data達33.53 Gb,以同屬華中中蜂的江西五府山中蜂基因組為參考,篩選出了3 886 321個高質量的SNP和1 392 028個small indel。其中,武隆中蜂具有相同變異位點且發(fā)生非同義SNP突變的基因有589個,CDS區(qū)發(fā)生small indel的基因214個。另外,還發(fā)現(xiàn)16個具有相同SNP和small indel的變異基因,表明這些突變基因可能在武隆中蜂適應本地自然環(huán)境的過程中發(fā)揮著重要作用。

武隆中蜂具有相同變異位點且發(fā)生非同義SNP突變的基因主要富集在賴氨酸降解的KEGG通路,CDS區(qū)發(fā)生small indel的基因富集在軸突引導、細胞外基質受體互作和MAPK信號通路。賴氨酸在蜜蜂生長發(fā)育、脂肪代謝和免疫力調控等方面發(fā)揮著重要作用[17-18]。在機體內,賴氨酸的降解主要通過酵母氨酸途徑和哌啶酸途徑轉換為乙酰輔酶A,進而在TCA循環(huán)中參與營養(yǎng)物質代謝,為機體供能并增強免疫力[19-20]。軸突導向是神經(jīng)元發(fā)育中的基本過程,在神經(jīng)網(wǎng)絡的形成中發(fā)揮著重要作用[21-22],細胞外基質受體互作與昆蟲腦部的發(fā)育相關[23],軸突導向和細胞外基質受體互作通路的變化可能涉及神經(jīng)系統(tǒng)的進化,進而調節(jié)昆蟲的運動、飛行和趨避等行為。MAPK信號通路是最重要的信號轉導系統(tǒng)之一,介導細胞增殖和機體生長發(fā)育等多種生理過程[24]。近期研究表明,蜂王長壽可能與MAPK信號通路增強細胞活性相關[25],此外,該通路還參與蜜蜂對病原微生物和低溫脅迫的響應[26-27],說明武隆中蜂優(yōu)良的抗逆性能可能與MAPK信號通路相關基因的突變有關。

發(fā)達的嗅覺系統(tǒng)是昆蟲適應環(huán)境的重要體現(xiàn),它介導蜜蜂的信息交流、定位和采集等行為,是中蜂擅長采集零星蜜源的基礎[28-29]。嗅覺受體蛋白參與昆蟲的嗅覺識別過程,它們與外界氣味分子和氣味結合蛋白的復合體結合,將化學信號轉換為電信號并傳遞到大腦中樞神經(jīng)系統(tǒng)[28,30]。細胞色素P450是蜜蜂解毒代謝外源物質的關鍵酶,能在溫和條件下催化多種氧化反應,改變外源物質的結構,降低其毒性,并使其從親脂性轉變?yōu)橛H水性[31-32]。研究表明,P450參與蜜蜂對植物次生物質和多種農(nóng)藥的解毒代謝[33-34]。此外,細胞色素P450基因(CYP6a17)還參與調控昆蟲對溫度的偏好,這種偏好是果蠅(Drosophilamelanogaster)等變溫動物優(yōu)化體溫的重要行為策略[35]。本試驗篩選出了16個在武隆中蜂中具有相同SNP和small indel的突變基因,其中包括2個嗅覺受體基因和1個細胞色素P450基因(CYP6a17),這暗示嗅覺進化、解毒代謝和對溫度偏好性的轉變可能是武隆中蜂適應重慶本地自然環(huán)境的重要機制。

目前,我國中蜂分為9個地理型,包括北方中蜂、華南中蜂、華中中蜂、云貴高原中蜂、長白山中蜂、海南中蜂、滇南中蜂、阿壩中蜂和西藏中蜂[36]。由于長期的地理隔離,海南中蜂、阿壩中蜂、西藏中蜂和長白山中蜂與其他地理型中蜂的親緣關系較遠[15,37-38]。本研究通過對系統(tǒng)進化樹、遺傳結構和遺傳分化的分析,成功區(qū)分了這4種中蜂與其他樣本,說明長期的地理隔離是推動中蜂種群分化的關鍵,能有效防止外來遺傳資源的干擾和污染。此外,阿壩中蜂、西藏中蜂和長白山中蜂長期生活在高原地區(qū),這些地區(qū)的生境、地形復雜,植物多樣,為中蜂的繁衍和發(fā)展提供了重要條件,也促進了中蜂逐漸適應高原生境。唐相友等[14]的研究表明,在我國的一些高原地區(qū),中蜂種群存在較高的遺傳分化水平,這暗示在青藏高原等地區(qū)可能存在更多具有獨特地理遺傳特性的中蜂群體。

重慶是我國中蜂的主要養(yǎng)殖區(qū)。本研究結果顯示,在系統(tǒng)進化樹和遺傳結構的主成分分析中,武隆中蜂與我國其他地理型中蜂有顯著的遺傳分化,說明武隆中蜂是重慶地區(qū)獨特的中蜂種質遺傳資源。根據(jù)第3次全國畜禽遺傳資源普查,重慶中蜂屬于華中中蜂。本研究中,武隆中蜂與湖北神農(nóng)架中蜂群體的遺傳分化水平最低,且在進化樹中歸為同一支,說明兩個同屬華中中蜂的種群親緣關系最近。其次,武隆中蜂與陜西安康的北方中蜂在進化樹和Fst分析中也顯示出較近的親緣關系。盡管武隆中蜂與云貴高原中蜂和廣西河池的華南中蜂的遺傳分化水平最低,但在進化樹中屬于不同的大支,因此,推測武隆中蜂與云貴高原中蜂和華南中蜂的親緣關系遠于華中中蜂和北方中蜂。武隆中蜂親本蜂群與人工繁育的F1代蜂群群體之間的Fst最小,且兩個群體在進化樹中歸為同一分支,表明初步的人工繁育并未對武隆中蜂的遺傳多樣性產(chǎn)生影響。在群體遺傳結構分析中,在K=10時,僅能將西藏中蜂、長白山中蜂、海南中蜂和阿壩中蜂區(qū)分出來,暗示其他地區(qū)中蜂可能存在自然或人為干預的遺傳資源交流。多年來,重慶本地蜂農(nóng)為提高養(yǎng)殖效益,隨意從其他地區(qū)引進中蜂種群,污染了本地遺傳資源,導致品種雜化和生產(chǎn)性能快速下降。

綜上所述,本研究揭示了重慶武隆中蜂的遺傳變異和種群分化特征,并推測了與適應性進化相關的功能通路和基因,為進一步挖掘、保護和開發(fā)利用重慶本地中蜂種質資源提供了重要的理論依據(jù)。同時,在引入外來中蜂種群時要注意避免對本地遺傳資源的干擾和污染,以保護和提升中蜂的生產(chǎn)性能和適應性。