全基因組重測序揭示青藏高原東緣及東南緣東方蜜蜂遺傳多樣性與適應(yīng)性進(jìn)化

唐相友, 宋華莉, 石 鵬, 張小燕, 唐紫寒,王文峰, 扎 羅, 陳新蘭, 周澤揚(yáng), 許金山,*

(1. 重慶師范大學(xué)生命科學(xué)學(xué)院, 重慶 401331; 2. 活性物質(zhì)生物技術(shù)教育部工程研究中心，重慶 401331;3. 重慶市媒介昆蟲重點實驗室，重慶 401331; 4. 西藏農(nóng)牧學(xué)院, 西藏林芝 860000)

青藏高原東緣及東南緣處于青藏高原邊緣地帶，因其地形復(fù)雜、生態(tài)境多樣化，為各類生物生態(tài)、繁衍發(fā)展提供了條件，即是全球生物多樣性熱點地區(qū)之一，也是研究物種形成和分化機(jī)制的理想?yún)^(qū)域(王鳳英等, 2005; 劉天猛, 2018; 張劍搏等, 2019)。青藏高原東緣及東南緣東方蜜蜂Apiscerana適應(yīng)了高原復(fù)雜多變的環(huán)境，其優(yōu)質(zhì)種質(zhì)資源對于蜜蜂遺傳改良具有重要意義。因此，開展高原東方蜜蜂遺傳多樣性及其環(huán)境適應(yīng)性的分子機(jī)制研究顯得尤為重要。此前，學(xué)者主要利用形態(tài)學(xué)、線粒體DNA和微衛(wèi)星分子標(biāo)記的方法探討了青藏高原東方蜜蜂的遺傳多樣性，表明在此區(qū)域存在阿壩蜜蜂Apisceranaabansis和藏南蜜蜂Apisceranaskorikoui兩個不同類型的高原東方蜜蜂遺傳群體(楊冠煌等, 1986; Zhuetal., 2017; Yuetal., 2019; 周丹銀等, 2019)。

隨著我國東方蜜蜂全基因組精細(xì)圖譜的完成與公開(Wangetal., 2020; Lanetal., 2021)，通過全基因組重測序鑒定群體單核苷酸多態(tài)性(SNP)變異位點，對特定區(qū)域的東方蜜蜂遺傳多樣性進(jìn)行更加全面的分析，彌補(bǔ)傳統(tǒng)形態(tài)學(xué)與單一分子標(biāo)記在東方蜜蜂遺傳分化研究中的不足，對挖掘和利用我國蜜蜂遺傳資源具有促進(jìn)作用(Chen NBetal., 2018; Jietal., 2020; Shietal., 2020)。當(dāng)前，尚缺乏對青藏高原東緣及東南緣東方蜜蜂的遺傳多樣性的深入研究，區(qū)域內(nèi)不同群體之間的近緣關(guān)系及擴(kuò)散來源也暫未厘清。因此，本研究采集了青藏高原東緣和東南緣及鄰近地區(qū)的東方蜜蜂樣本并進(jìn)行全基因組重測序，結(jié)合此前已公開的部分?jǐn)?shù)據(jù)，擬通過群體遺傳結(jié)構(gòu)及單倍型分析，研究青藏高原東緣和東南緣東方蜜蜂遺傳多樣性，并結(jié)合受選擇分析初步鑒定青藏高原東緣及東南緣東方蜜蜂高原適應(yīng)性相關(guān)基因，初步探討其潛在的分子適應(yīng)機(jī)制，以期為進(jìn)一步解析青藏高原東方蜜蜂的遺傳資源多樣性、種群擴(kuò)散規(guī)律以及適應(yīng)高原生境下的分子進(jìn)化機(jī)制提供參考。

1 材料與方法

1.1 樣本采集與收集

采集青藏高原東緣和東南緣及鄰近地區(qū)的77群東方蜜蜂，保存至75%的酒精中。從每群中隨機(jī)抽取1頭工蜂，取頭、胸部組織提取總DNA，4℃保存待測序；從GenBank數(shù)據(jù)庫下載青藏高原東緣和東南緣及鄰近地區(qū)的90群東方蜜蜂全基因組重測序原始數(shù)據(jù)(GenBank登錄號: PRJNA418874)，上述167群分屬于11個地理種群，范圍覆蓋了青藏高原東緣(邛崍山脈東西、南北延伸以及四川九寨溝)及東南緣(西藏波密以及云南迪慶和西雙版納)區(qū)域(表1)，以西域黑蜂Apismelliferasinisxinyuan全基因組重測序數(shù)據(jù)(GenBank登錄號: PRJNA301648)作為外群。

1.2 全基因組重測序與SNP calling

委托公司使用Illumina HiSeq測序平臺進(jìn)行對1.1節(jié)采集的77群東方蜜蜂進(jìn)行全基因組重測序(GenBank登錄號: PRJNA827659)。使用FASTX-Toolkit軟件(http:∥hanno nlab.cshl.edu/fastx_toolk it/)去除低質(zhì)量(Q<20)序列和含有5%不可讀堿基(N堿基)以上的原始reads，獲得高質(zhì)量可用于后續(xù)遺傳分析的reads。利用比對軟件BWA-MEM，將樣本reads匹配到高質(zhì)量東方蜜蜂參考基因組上(Lanetal., 2021)，采用bamdst計算樣本的匹配率和測序深度。采用SAMtools(Li, 2011)軟件對比對BAM文件進(jìn)行排序處理，而后使用GATK4.0(do Valleetal., 2016)對比對BAM文件去除重復(fù)讀取和SNP calling，采用HaplotypeCaller和VariantFiltration程序包對SNPs進(jìn)行鑒定和過濾。采用VCFtools(Daneceketal., 2011)對SNP候選數(shù)據(jù)集進(jìn)行進(jìn)一步篩選，確定SNPs數(shù)據(jù)用于后續(xù)相關(guān)分析。

1.3 群體遺傳分析

基于高質(zhì)量SNPs數(shù)據(jù)，使用軟件ADMIXTURE(Alexanderetal., 2009)的基于塊松弛算法(block relaxation algorithm)分析群體遺傳結(jié)構(gòu)和個體祖先血緣混合情況。 使用EIGENSOFT(Priceetal., 2006)程序進(jìn)行主成分分析。采用PHYLIP(v3.697)鄰接法(neighbor-joining, NJ)基于SNPs構(gòu)建系統(tǒng)進(jìn)化樹，以西域黑蜂的SNPs作為外群，bootstrap值設(shè)置為1 000，使用在線工具iTOL(http:∥itol.embl.de/)進(jìn)行系統(tǒng)進(jìn)化樹可視化呈現(xiàn)。使用METAPOP2(López-Corteganoetal., 2019)計算種群間最小遺傳距離與種群內(nèi)成對遺傳分化指數(shù)(Fst)。將青藏高原東緣及東南緣高原東方蜜蜂分別與鄰近地區(qū)非高原東方蜜蜂組合，即：(1)藏南高原與滇北高原組：波密(BM)+迪慶(DQ)+西雙版納(XSBN)；(2)川西高原組：馬爾康(MEK)+丹巴(DB)+邛崍山(QLS)；(3)川北高原組：九寨溝(JZG)+神農(nóng)架(SNJ)。利用DnaSP6(Rozasetal., 2017)軟件進(jìn)行線粒體基因組單倍型分析，采用POPART(Leigh and Bryant, 2015)進(jìn)行圖形化展示。

1.4 選擇信號分析

為挖掘青藏高原東方蜜蜂環(huán)境適應(yīng)性相關(guān)基因，將青藏高原東緣及東南緣高原東方蜜蜂分別與遺傳關(guān)系較近的非高原東方蜜蜂組合，即：(1)藏南高原組：波密(BM)+西雙版納(XSBN)；(2)滇北高原組：迪慶(DQ)+西雙版納(XSBN)；(3)川西高原組：丹巴(DB)+邛崍山(QLS)；(4)川北高原組：九寨溝(JZG)+神農(nóng)架(SNJ)。采用VCFTOOLS(Daneceketal., 2011)以滑動窗口大小50 kb，步長10 kb計算核苷酸多態(tài)性(θπ)和遺傳分化指數(shù)(Fst)，采用XPCLR(Chenetal., 2010)軟件計算XP-CLR值，滑動窗口大小設(shè)置為50 kb，步長10 kb，每代堿基突變率設(shè)置為5.3e-9。基于東方蜜蜂阿壩株的全基因組序列及注釋文件(Lanetal., 2021)提取上述組合中兩群體基因組中所有潛在受選擇的基因序列，通過Eggnog-mapper(Cantalapiedraetal., 2021)進(jìn)行GO注釋和KEGG通路分析，采用FDR進(jìn)行校正。

2 結(jié)果

2.1 SNP鑒定

與東方蜜蜂阿壩株全基因組精細(xì)圖譜序列(Lanetal., 2021)進(jìn)行比對發(fā)現(xiàn)，167群東方蜜蜂樣本全基因組重測序數(shù)據(jù)與參考基因組序列的平均匹配率為91.05%，平均測序深度為10.27×，數(shù)據(jù)質(zhì)量達(dá)成后續(xù)分析要求。通過全基因組SNP calling，共獲得137萬個高質(zhì)量的SNPs(表1)。

表1 基于全基因組重測序數(shù)據(jù)的青藏高原東緣、東南緣及鄰近地區(qū)東方蜜蜂地理種群的SNPs數(shù)量與核苷酸多態(tài)性(θπ)Table 1 Numbers of SNPs and nucleotide polymorphisms (θπ) of the geographic populations of Apis ceranaon the eastern and southeastern edges of the Qinghai-Tibet Plateau and adjacent areasbased on whole-genome resequencing data

2.2 遺傳分化

基于167群東方蜜蜂樣本基因組SNPs，當(dāng)假設(shè)樣本中存在3個群體時(K=3)，可區(qū)分出來自川西高原和藏南高原地理群體；當(dāng)K=7時，來自滇北高原和川北高原地理群體依次從混合地理群體中區(qū)分出來(圖1: A)。 這表明來自川西高原和川北高原的群體與藏南高原群體和滇北高原群體在遺傳結(jié)構(gòu)上具有明顯的群體分化。主成分結(jié)果(圖1: B)顯示，167群東方蜜蜂樣本被劃分為川西高原、川北高原、藏南高原、滇北高原和其他地理群體。以上結(jié)果表明青藏高原東緣及東南緣分布著高原東方蜜蜂4個遺傳分化群體，應(yīng)該與該區(qū)域廣闊且特殊的地理生境相關(guān)。

圖1 基于SNPs的青藏高原東緣、東南緣及鄰近地區(qū)東方蜜蜂遺傳結(jié)構(gòu)隨分組組數(shù)(K)的變化(A)和主成分分析(B)Fig. 1 Genetic structure change with grouping number (K) (A) and principal component analysis (B) of Apis ceranaon the eastern and southeastern edges of the Qinghai-Tibet Plateau and adjacent areas based on SNPs種群信息見表1。For population information, see Table 1. WSichPl: 川西高原Western Sichuan Plateau; STibetPl: 藏南高原Southern Tibet Plateau; NYnanPl: 滇北高原Northern Yunnan Plateau; NSichPl: 川北高原Northern Sichuan Plateau. 下同The same below.

基于東方蜜蜂成對種群的遺傳分化指數(shù)分析(表2)顯示，11個種群的成對Fst值在0.0199～0.1443之間，平均0.0758。來自川西高原、川北高原、藏南高原和滇北高原的4個高原地理群體之間的Fst值(Fst=0.0917～0.1443,平均0.1178)高于非高原群體之間的Fst值(Fst=0.0219～0.0603，平均0.0411)，這表明高原地理群體之間產(chǎn)生了較大的遺傳分化，這與遺傳結(jié)構(gòu)分析(圖1: A)和主成分分析(圖1: B)結(jié)果一致。

不同種群間的最小遺傳距離分析結(jié)果(表2)顯示，相較于其他種群，DB和MEK與QLS之間的最小遺傳距離值最小；BM和DQ與XSBN之間最小遺傳距離最??；JZG與SNJ之間最小遺傳距離最小。結(jié)合這些種群所屬的地理位置，可以表明川西高原群體與川西山地群體具有更近的遺傳關(guān)系；藏南高原、滇北高原群體與滇南群體具有更近的遺傳關(guān)系；川北高原群體與秦巴群體具有更近的遺傳關(guān)系。

表2 基于SNPs的青藏高原東緣、東南緣及鄰近地區(qū)東方蜜蜂地理種群遺傳分化指數(shù)(Fst)(上三角)和最小遺傳距離(Nm)(下三角)Table 2 Pairwise genetic differentiation index (Fst)(above diagonal) and the minimum genetic distance (Nm) (below diagonal) of the geographic populations of Apis cerana on the eastern and southeastern edgesof the Qinghai-Tibet Plateau and adjacent areas based on SNPs

2.3 系統(tǒng)進(jìn)化

基于SNPs構(gòu)建的系統(tǒng)進(jìn)化樹(圖2)顯示，東方蜜蜂的青藏高原東緣及東南緣的高原群分為兩個姐妹支，即東緣支和東南緣支，其中東南緣支包含藏南高原與滇北高原地理群體，東緣支則包含川西高原與川北高原地理群體。由此表明，對于東方蜜蜂高原群體而言，在地理位置上更接近的群體在親緣關(guān)系上也更加接近。

圖2 鄰接法構(gòu)建的基于SNPs的青藏高原東緣、東南緣的高原東方蜜蜂群的系統(tǒng)進(jìn)化樹(1 000次重復(fù))Fig. 2 Phylogenetic tree of the plateau colonies of Apis cerana on the eastern and southeastern edges of the Qinghai-TibetPlateau constructed by neighbor-joining method based on SNPs (1 000 replicates)以西域黑蜂的SNPs為外群。The SNPs of Apis mellifera sinisxinyuan were used as the outgroup.

2.4 線粒體基因組單倍型

基于東方蜜蜂線粒體基因組單倍型分析，分別在川西高原組(MEK+DB+QLS)(圖3: A)、川北高原組(JZG+SNJ) (圖3: B)和藏南高原與滇北高原組(BM+DQ+XSBN) (圖3: C)線粒體基因組中發(fā)現(xiàn)24, 16和21個單倍型。其中，川西高原組發(fā)現(xiàn)祖先單倍型Hap_15分布于QLS中，MEK與DB的單倍型為衍生單倍型；川北高原組中祖先單倍型Hap_6分布于SNJ中，JZG的單倍型為衍生單倍型；而在藏南高原與滇北高原組中XSBN的單倍型為主要單倍型，BM和DQ的單倍型為衍生單倍型。這為初步推斷青藏高原東緣及東南緣高原東方蜜蜂來源提供了有效證據(jù)。

圖3 青藏高原東緣、東南緣及鄰近地區(qū)東方蜜蜂地理種群線粒體基因組單倍型網(wǎng)絡(luò)圖Fig. 3 Haplotype network of the mitochondrial genome of the geographic populations of Apis ceranaon the eastern and southeastern edges of the Qinghai-Tibet Plateau and adjacent areasA: 川西高原Western Sichuan Plateau; B: 川北高原Northern Sichuan Plateau; C: 藏南高原與滇北高原Southern Tibet Plateau and Northern Yunnan Plateau. Hap: 單倍型Haplotype.

2.5 環(huán)境適應(yīng)性相關(guān)基因挖掘

在藏南高原組(BM+XSBN)基因組中發(fā)現(xiàn)71個潛在受選擇區(qū)域和169個潛在受選擇基因，主要富集在細(xì)胞脂質(zhì)代謝過程、硫化合物生物合成等GO功能條目(圖4)，參與Wnt, FoxO, Notch和Hippo等KEGG通路(圖5)；在滇北高原組(DQ+XSBN)基因組中發(fā)現(xiàn)22個潛在受選擇區(qū)域和53個潛在受選擇基因，主要富集在防御反應(yīng)調(diào)節(jié)、免疫反應(yīng)調(diào)節(jié)等GO功能條目，參與FoxO等KEGG通路；在川西高原組(DB+QLS)基因組中發(fā)現(xiàn)40個潛在受選擇區(qū)域和40個潛在受選擇基因，主要富集在防御反應(yīng)調(diào)節(jié)、DNA修復(fù)等GO功能條目，參與Toll, Imd, FoxO, MAPK和Hedgehog等KEGG通路；在川北高原組(JZG+SNJ)基因組中發(fā)現(xiàn)4個潛在受選擇區(qū)域和10個潛在受選擇基因，主要富集在神經(jīng)肌肉的過程、大腦皮層發(fā)育等GO功能條目，參與線粒體自噬等KEGG通路。總體來說，這些受選擇基因參與脂肪酸代謝、光轉(zhuǎn)導(dǎo)、溫度適應(yīng)、卵巢發(fā)育等信號通路，但在功能富集區(qū)域又明顯不同，這顯示出不同高原東方蜜蜂在適應(yīng)不同生境具有較大的基因選擇差異。

圖4 藏南高原組(BM+XSBN)東方蜜蜂基因組中受選擇基因的GO功能富集Fig. 4 GO functional enrichment of selected genes in the genomes of Apis ceranain the plateau group (BM+XSBN) of southern Tibet

圖5 藏南高原組(BM+XSBN)東方蜜蜂基因組中受選擇基因KEGG通路富集Fig. 5 KEGG pathway enrichment of selected genes in the genomes of Apis ceranain the plateau group (BM+XSBN) of southern Tibet

在滇南的藏南高原與滇北高原東方蜜蜂基因組中發(fā)現(xiàn)了5個共同的潛在受選擇基因，分別為突觸囊泡樣糖蛋白2B基因(GenBank登錄號: 107993575), UNC93樣蛋白MFSD11基因(GenBank登錄號: 107994554), 叉頭盒蛋白O(FOXO)基因(GenBank登錄號: 107993192), 胰島素基因增強(qiáng)蛋白-1(ISL-1)基因(GenBank登錄號: 107995439)以及一個未注釋基因LOC107993577 (GenBank登錄號: 107993577)，其中ISL-1和FOXO，主要參與胰島素分泌以應(yīng)對細(xì)胞壓力。

3 討論

在本研究對青藏高原東緣及東南緣及鄰近地區(qū)的11個地理種群的167群東方蜜蜂進(jìn)行全基因組重測序數(shù)據(jù)分析，基于群體遺傳結(jié)構(gòu)和主成分分析發(fā)現(xiàn)，高原東方蜜蜂群分屬川西高原、藏南高原、滇北高原、川北高原(圖1: A, B)，分布在青藏高原東緣及東南緣帶上。遺傳分化指數(shù)表明(表2)，高原東方蜜蜂的遺傳分化指數(shù)明顯高于非高原地區(qū)的，且部分地區(qū)遺傳分化指數(shù)也高于西方蜜蜂Apismellifera亞種間分化水平(平均Fst值約為0.1)(Wallbergetal., 2014)，表明部分高原東方蜜蜂存在較大的遺傳分化水平，達(dá)到了亞種分化水平。我們的研究結(jié)果也暗示著在青藏高原很可能存在更多的獨特的地理遺傳地理群體。

鑒于青藏高原地區(qū)多高山與峽谷，地處于這些地方的東方蜜蜂與鄰近地區(qū)的東方蜜蜂的聯(lián)系尚且未知。我們的研究基于各地理種群間遺傳距離暗示了青藏高原東緣及東南緣東方蜜蜂與鄰近地區(qū)東方蜜蜂的關(guān)系，研究表明藏南高原、滇北高原和滇南的東方蜜蜂具有更近的遺傳關(guān)系，川西高原與川西山地的東方蜜蜂具有更近的遺傳關(guān)系，而川北高原與秦巴神農(nóng)架林區(qū)的東方蜜蜂具有更近的遺傳關(guān)系(表2)。因此，結(jié)合本研究結(jié)果與此前相關(guān)研究(周丹銀等, 2019; Jietal., 2020; Shietal., 2020)，我們初步推測了青藏高原東緣及東南緣不同東方蜜蜂地理群體來源，即處于青藏高原東南緣的藏南高原、滇北高原東方蜜蜂來源于滇南，東緣的川西高原、川北高原東方蜜蜂分別來源于川西山地、秦巴，暗示了青藏高原東緣及東南緣東方蜜蜂是由鄰近非高原東方蜜蜂的種群擴(kuò)散后的地理隔離形成而產(chǎn)生的獨立分化起源。

川西高原平均海拔約為3 740 m，主要屬于溫帶-寒溫帶季風(fēng)氣候；滇北高原平均海拔約為3 380 m，氣候類型與川西高原相似；藏南高原平均海拔超過4 500 m，主要屬于亞熱帶山地濕潤季風(fēng)氣候。氣候多樣性與地貌差異形成了青藏高原東緣及東南緣不同地理型的東方蜜蜂分化，而這種分化又與適應(yīng)性進(jìn)化具有一定的關(guān)聯(lián)。我們基于對青藏高原東緣及東南緣高原東方蜜蜂與鄰近非高原東方蜜蜂的選擇性分析，發(fā)現(xiàn)潛在受選擇基因參與多個信號通路(圖5)，其中，Wnt, MAPK, Notch以及Hippo信號通路在蜜蜂等昆蟲的群體分化、胚胎形成、形態(tài)發(fā)生、成形盤發(fā)育和器官大小調(diào)節(jié)等發(fā)育過程中發(fā)揮重要作用(Deardenetal., 2006; Bolósetal., 2007; Komiya and Habas, 2008; Halder and Johnson, 2011; Wilsonetal., 2011; Yinetal., 2018)。Wnt, Notch和Hedgehog信號通路被認(rèn)為與卵巢活動有關(guān)(Duncanetal., 2016; Chenetal., 2017; Chen Cetal., 2018)，這可能顯示出高原型較非高原型東方蜜蜂的卵巢活動更加強(qiáng)烈。Toll和Imd信號通路主要涉及包括蜜蜂在內(nèi)的昆蟲的免疫能力(Aneteetal., 2021)。Hippo信號通路在西方蜜蜂中被認(rèn)為是適應(yīng)寒冷氣候的重要通路(Chen Cetal., 2016, 2018; Shietal., 2020)，這顯示出這些基因可能與高原東方蜜蜂耐寒性相關(guān)。除此之外，F(xiàn)oxO信號通路也被認(rèn)為參與溫度適應(yīng)(Chen Cetal., 2016, 2018)。同時，我們還發(fā)現(xiàn)了潛在受選擇基因參與代謝通路、光轉(zhuǎn)導(dǎo)等相關(guān)通路，這些通路可能為高原東方蜜蜂適應(yīng)不同生境提供了基礎(chǔ)?；趤碓从诘崮系牟啬细咴c滇北高原東方蜜蜂的潛在受選擇基因，我們發(fā)現(xiàn)了5個共同受選擇基因。有研究表明ISL-1在調(diào)節(jié)細(xì)胞生命和胚胎發(fā)育中發(fā)揮重要作用(Pfaffetal., 1996)，Evm.model.16.590 (ISL-1)和FOXO均參與調(diào)節(jié)胰島素分泌以應(yīng)對細(xì)胞壓力(Jüngeretal., 2003; Guoetal., 2021)，我們推測其為青藏高原東南緣的蜜蜂適應(yīng)當(dāng)?shù)丨h(huán)境所需要。這對我們初步了解東方蜜蜂適應(yīng)不同生境下的適應(yīng)性基因情況提供了基礎(chǔ)，為后續(xù)進(jìn)行蜜蜂種質(zhì)資源挖掘與遺傳改良提供了參考。