基于 28S rDNA D2區(qū)序列的陜南凹緣菱紋葉蟬分子系統(tǒng)學(xué)與遺傳多樣性研究

楊金宏，孔衛(wèi)青

(安康學(xué)院，陜西省蠶桑重點(diǎn)實(shí)驗(yàn)室, 陜西安康 725000)

凹緣菱紋葉蟬(HishimonussellatusUhler)俗稱綠頭菱紋葉蟬，隸屬于昆蟲綱(Insecta)半翅目(Hemiptera)葉蟬科(Cicadellidae)殃葉蟬亞科(Euscelinae)[1]。凹緣菱紋葉蟬是桑樹黃花型萎縮病、棗瘋病等的重要媒介昆蟲，在中國(guó)分布廣泛[2]。成蟲和若蟲均可造成危害。它們以口針刺吸病樹汁液，亦將病原吸入體內(nèi)，后者移至唾液腺并繁殖；再次取食時(shí)病原菌傳入新寄主，從而導(dǎo)致健康植株發(fā)病。桑樹染病后通常會(huì)出現(xiàn)花變?nèi)~、矮化、黃化、衰退、簇生、叢枝、小葉等癥狀，嚴(yán)重影響蠶桑產(chǎn)業(yè)的健康發(fā)展。凹緣菱紋葉蟬每年4月下旬至10月下旬發(fā)生危害，1 a達(dá)3～4代，繁殖速度快，世代周期短，世代重疊，危害嚴(yán)重，且與養(yǎng)蠶周期重合，難以防治。

秦巴山區(qū)有眾多的小盆地和山間谷地相連接，安康市漢濱、石泉、漢陰、平利、紫陽等縣(區(qū))和漢中市西鄉(xiāng)、勉縣等均位于該地理區(qū)域，形成陜南蠶桑產(chǎn)業(yè)主產(chǎn)區(qū)。各地自然地理?xiàng)l件差異較大，具有很高的環(huán)境異質(zhì)性。桑園種植多分布于坡地和山地，可能阻礙凹緣菱紋葉蟬擴(kuò)散，形成地理隔離，進(jìn)而導(dǎo)致不同種群間受到不同的選擇壓力而出現(xiàn)種群分化。隨著近年來國(guó)家退耕還林、土地流轉(zhuǎn)政策以及蠶桑產(chǎn)業(yè)生產(chǎn)模式的變更，較多的連片桑園和強(qiáng)村大戶出現(xiàn)，這為凹緣菱紋葉蟬擴(kuò)散和傳播提供了條件。為了解陜南桑園凹緣菱紋葉蟬種群的遺傳信息，分析其分布特點(diǎn)，本研究利用 28S rDNA基因D2區(qū)對(duì)其種群遺傳結(jié)構(gòu)和遺傳分化進(jìn)行分析。

28S rDNA存在于編碼細(xì)胞質(zhì)核糖體大亞基中，在生物體內(nèi)含量較高，在進(jìn)化過程中比較保守，堿基替換率低，是研究生物系統(tǒng)發(fā)育較好的分子標(biāo)記[3-4]。同時(shí)在 28S rDNA基因的保守序列中，含有12個(gè)高變區(qū)D1～D12，可解釋從屬到科級(jí)水平上的系統(tǒng)發(fā)生關(guān)系以及屬種間的親緣關(guān)系[5-8]，是非常有用的遺傳進(jìn)化分析工具。Dietrich等[9]用 28S rDNA D2～D10區(qū)序列分析角蟬總科(Membracoidea)與葉蟬科的親緣關(guān)系。戴仁懷等[10]利用 28S rDNA基因 D2區(qū)、 16S rDNA序列和形態(tài)特征對(duì)中國(guó)角頂葉蟬亞科(Deltocephalinae)的系統(tǒng)發(fā)育進(jìn)行研究。目前，未見利用 28S rDNA基因?qū)Π季壛饧y葉蟬進(jìn)行分子系統(tǒng)學(xué)和遺傳多樣性研究的報(bào)道。本試驗(yàn)嘗試這項(xiàng)研究，其結(jié)果將促進(jìn)對(duì)陜南蠶桑產(chǎn)區(qū)凹緣菱紋葉蟬的種群遺傳結(jié)構(gòu)及其分布概況的了解，并為凹緣菱紋葉蟬相近種的分類和快速分子診斷提供依據(jù)。

1 材料與方法

1.1 樣本采集

2017年9月至10月，用黃板誘蟲法采集秦巴山區(qū)陜南各蠶桑區(qū)凹緣菱紋葉蟬成蟲。采集的新鮮標(biāo)本浸泡于φ=75% 乙醇中，于12 h內(nèi)置于-20 ℃的冰箱保存，備用。樣品采集地點(diǎn)及相關(guān)信息見表1。

表1 陜南凹緣菱紋葉蟬種群樣本采集地及其數(shù)量Table 1 Locality and sample size of Hishimones sellatus populations in South Shaanxi

1.2 基因組DNA提取

1.3 PCR擴(kuò)增及測(cè)序

目標(biāo)片段擴(kuò)增用28S D2～D3引物[9]完成，引物序列為：28SD2F：5′-AGTCGKGTTGCTT-GAKAGTGCAG-3′；28SD2R：5′-TTCAATT-TCATTKCGCCTT-3′。反應(yīng)總體系50 μL，包括基因組DNA 2.0 μL，10×PCR擴(kuò)增緩沖液5.0 μL，Mg2+(25 mmol/μL)5.0 μL，dNTPs(10 mmol/μL)1.5 μL，上下游引物(10 pmol/μL)各1.0 μL，高保真TaqDNA聚合酶(5 U/μL)1.0 μL，加雙純水補(bǔ)足至50.0 μL。PCR擴(kuò)增程序?yàn)閇8]：94 ℃預(yù)變性3 min，36×(94 ℃變性1 min，52.2 ℃退火 1 min，72 ℃延伸50 s)，72 ℃終延伸10 min。取 3.0 μL PCR擴(kuò)增產(chǎn)物進(jìn)行10 g/L的瓊脂糖凝膠電泳，檢測(cè)成功后直接送生工生物工程(上海)股份有限公司進(jìn)行純化測(cè)序。

1.4 序列多樣性分析

以Clustal X 1.83軟件[11]基于鄰接法(Neighbor-joining)對(duì)所得序列進(jìn)行比對(duì)，去除兩端不準(zhǔn)確序列后，最終輸出fasta格式的對(duì)齊序列。利用Dnasp 61101軟件[12]統(tǒng)計(jì)各序列間的堿基組成、變異位點(diǎn)，以及基于對(duì)齊序列的種群核苷酸多樣性Pi[13]、單倍型多樣性Hd[Hd=n(1-∑Pi2)/(n-1)，其中n為樣本數(shù)]和種群遺傳分化系數(shù)等，利用Tajima’sD方法對(duì)序列進(jìn)行中性檢驗(yàn)，并根據(jù)數(shù)值服從β分布的特點(diǎn)，利用隨機(jī)抽樣模擬的種群中性進(jìn)化為零假設(shè)，通過與觀測(cè)值的比較進(jìn)行顯著性檢驗(yàn)[14]。

1.5 種群間遺傳差異分析

利用Arlequin 3.5.22軟件[15]，基于Locus-by-Locus模型分析種群變異的分子方差A(yù)MOVA(Analysis of molecular variance)以及種群間的固定系數(shù)Fst(F-statistics)，以成對(duì)單倍型差異比較法(Pairwise difference)計(jì)算種群間遺傳距離。運(yùn)行Dnasp 61101軟件[12]計(jì)算種群間的基因流系數(shù)Nst(Population subdivision at the nucleotide level)和種群間的地理距離(D)，Nst= 1 /(1 + 2Nm)；D= 111.199{(a1-a2)2+ (b1-b2)2cos2[(a1+a2) / 2]}1/2，其中Nm為有效遷移個(gè)體數(shù)，a和b分別代表經(jīng)度和緯度，并進(jìn)一步計(jì)算D的自然對(duì)數(shù)值以及與Nst之間的決定系數(shù)(R2)。

采用Bandelt等[16]的Median-joining模型繪制單倍型間的網(wǎng)絡(luò)進(jìn)化圖，以葉蟬科角頂葉蟬亞科的Dwightlaacutipennis作外群，MEGA 7.0.26軟件[17]運(yùn)行ML極大似然法構(gòu)建各個(gè)單倍型的系統(tǒng)進(jìn)化樹，采用的堿基替換模型為GTR+G，各分支置信度以1 000次自展重復(fù)進(jìn)行檢驗(yàn)。

2 結(jié)果與分析

2.1 凹緣菱紋葉蟬 28S rDNA D2區(qū)序列的特征及多樣性

經(jīng)檢測(cè)，提取的132頭凹緣菱紋葉蟬的總DNA全部合格； 28S rDNA基因D2區(qū)片段的PCR擴(kuò)增和純化測(cè)序，均獲得良好的測(cè)序結(jié)果。測(cè)序序列比對(duì)并去除2端不準(zhǔn)確序列后，得到包括缺失位點(diǎn)的有效片段695 bp(Alignment score=138 469)，序列比對(duì)結(jié)果顯著性檢驗(yàn)P值為 0.001～0.041。綜合所有的多樣性位點(diǎn)共定義8個(gè)單倍型(圖1)，其中存在多態(tài)性位點(diǎn)63個(gè)，突變67個(gè)；其中4個(gè)位點(diǎn)存在3種變體，59個(gè)位點(diǎn)有2種變體。分析各單倍型序列數(shù)和地理分布，H_1有98條序列(表2)，占序列總數(shù)的 74.2%，且都存在于所有7個(gè)地理種群中，為優(yōu)勢(shì)共享單倍型；H_2和H_3僅在1個(gè)種群中出現(xiàn)，為私有單倍型。各種群內(nèi)單倍型的分布從1種～6種不等，其中3號(hào)種群最多6種，且出現(xiàn)29個(gè)個(gè)體共享1個(gè)單倍型，具有較高的保守性。核苷酸位點(diǎn)多樣性Pi和單倍型多樣性指數(shù)Hd分別從核苷酸的突變位點(diǎn)數(shù)和差異基因型兩個(gè)層面展示種群的遺傳多樣性。7個(gè)種群的Pi分布為0～ 0.019 03，Pi最高的是3號(hào)種群，Hd為0～ 0.581 82，7號(hào)種群最高(表2)。Tajima’sD中性檢驗(yàn)是基于等位基因頻率和單倍型數(shù)目變化差異的一種種內(nèi)自然選擇檢驗(yàn)方法，對(duì)每一個(gè)種群和總?cè)后w進(jìn)行檢驗(yàn)(表2)，發(fā)現(xiàn)除4號(hào)種群大于0外，其他種群均為負(fù)數(shù)，其中5號(hào)種群為 -2.222 40，且顯著性檢驗(yàn)P<0.01，說明該種群受到強(qiáng)烈的負(fù)選擇影響[14]。

圖1 28S rDNA基因D2區(qū)序列單倍型及變異位點(diǎn)信息Fig.1 Haplotypes and variation sites of 28S rDNA D2 region sequence

地理種群Geographic population核苷酸多樣性(Pi)Nucleotide diversity單倍型多樣性(Hd)Haplotype diversity單倍型數(shù)量及其分布Haplotype number and distribution中性檢驗(yàn)Tajima’s DP值P value10.015 240.466 67H_1(15), H_7(2), H_8(4), -0.520 080.320200H_1(7)30.019 030.553 91H_1(29), H_4(3), H_5(2), H_6(3), H_7(3), H_8(4)-0.213 670.42940.013 680.485 29H_1(12), H_3(3), H_8(2)0.260 430.67350.006 510.166 67H_1(11), H_5(1)-2.222 40*0.00160.008 420.278 95H_1(17), H_4(1), H_6(2)-0.798 470.27970.013 260.581 82H_1(7), H_2(2), H_5(2)-0.191 360.445合計(jì) Total0.013 780.440 40H_1 ～ H_8-0.742 410.450

注：“( )”中的數(shù)字代表葉蟬數(shù)目。“*”表示差異顯著。

Note: Numbers in “( )” represent the number of leafhoppers.“*” indicate significant difference.

2.2 種群遺傳差異分析

Excel軟件計(jì)算基因流系數(shù)Nst和種群地理距離D的自然對(duì)數(shù)間的決定系數(shù)(R2)，可知R2=0.019 7(P<0.001)，接近于0，說明該基因在本研究范圍內(nèi)與各個(gè)種群之間的基因流動(dòng)和地理距離不相關(guān)，可能存在網(wǎng)狀基因交流。

2.3 單倍型進(jìn)化分析

Median-joining模型是在maximum-parsimony模型[20]和Kruskal’s算法[21]基礎(chǔ)上發(fā)展的，據(jù)此繪制的進(jìn)化圖呈網(wǎng)狀分布。本研究的凹緣菱紋葉蟬 28S rDNA D2區(qū)8種單倍型間的進(jìn)化網(wǎng)絡(luò)圖呈線狀分布(圖2)，在所有的8個(gè)單倍型中，只有H_4和H_1是種內(nèi)網(wǎng)絡(luò)進(jìn)化圖中的中間單倍型[22]，其他單倍型均處于進(jìn)化網(wǎng)絡(luò)邊緣，通過6個(gè)缺失的中間單倍型(mv1～mv6)連接，這與各單倍型在各個(gè)地理種群的分布結(jié)果一致(表2)，說明凹緣菱紋葉蟬的地理聚集性不強(qiáng)，交叉散布于整個(gè)進(jìn)化網(wǎng)絡(luò)圖中。

灰色圓圈代表單倍型 Grey circle represents the haplotype; 灰色圓圈面積代表單倍型相對(duì)頻率 The size of grey circle represents the relative frequency; 實(shí)心黑點(diǎn)代表中間單倍型 Solid black spots represent intermediate haplotype; 文本框內(nèi)數(shù)字代表單倍型間的突變位點(diǎn)數(shù) Boxed numbers represent the mutation sites between haplotypes

圖2 基于 28S rDNA基因的D2區(qū)的8個(gè)單倍型的Median-joining 網(wǎng)絡(luò)關(guān)系圖
Fig.2 Median-joining network of eight haplotypes based on 28S rDNA Gene D2 region sequences

進(jìn)一步以葉蟬科角頂葉蟬亞科的Dwightlaacutipennis的 28S rDNA D2區(qū)(JX845493.1)作外群，構(gòu)建單倍型的ML系統(tǒng)進(jìn)化樹，結(jié)果其 log likelihood值為-1 474.09。自展支持值bootstrap法是對(duì)初始DNA數(shù)據(jù)進(jìn)行放回抽樣生成偽樣本，進(jìn)而估計(jì)節(jié)點(diǎn)可靠度的一種非參數(shù)方法，本研究節(jié)點(diǎn)的自展支持值均大于90，且與網(wǎng)絡(luò)進(jìn)化圖的聚集特征一致，進(jìn)化關(guān)系較為可靠。根據(jù)外群的節(jié)點(diǎn)位置，進(jìn)化樹中單倍型H_1、H_2、H_7、H_8形成一個(gè)分支，H_3～H_6形成一個(gè)分支，其中H_7的枝長(zhǎng)最長(zhǎng)，堿基替換率最高，說明其與其他單倍型的實(shí)際進(jìn)化距離最遠(yuǎn)(圖3)。

枝長(zhǎng)代表堿基替換率 Branch lengths measured in the number of substitutions per site

圖3 基于 28S rDNA D2區(qū)的單倍型的ML系統(tǒng)發(fā)育樹
Fig.3 ML phylogenetic tree of haplotypes based on 28S rDNA D2 region

3 討 論

種群為適應(yīng)不同地域的自然環(huán)境而獨(dú)自演化，可能出現(xiàn)地理種群的遺傳分化和形態(tài)變異[23]，基于 28S rDNA基因D2區(qū)對(duì)陜南桑園凹緣菱紋葉蟬各種群和總?cè)后w進(jìn)行Tajima’sD中性檢驗(yàn)和顯著性分析，結(jié)果表明該種群和總?cè)后w近期沒有受到正選擇的影響，5號(hào)種群受到負(fù)選擇的影響，種群中存在較多的低頻率等位基因，而其他種群等位基因頻率處于穩(wěn)定狀態(tài)，遵循中性進(jìn)化模型[14]。陜南桑園凹緣菱紋葉蟬 28S r DNA基因D2區(qū)的核苷酸多樣性指數(shù)處于較低水平，變異主要來自種群內(nèi)且差異不顯著，表明該群體的大部分差異是隨機(jī)分布，且存在廣泛的基因交流。

選擇恰當(dāng)?shù)耐馊杭从H緣關(guān)系合適的外源物種，是進(jìn)化樹分析的關(guān)鍵，外群選擇不當(dāng)，容易造成進(jìn)化樹的拓?fù)浣Y(jié)構(gòu)發(fā)生變化[24]。葉蟬科被分為16個(gè)亞科，亞科內(nèi)部的親緣關(guān)系知之甚少，且不同亞科間可能存在許多并系[25]。因此本研究首先以蝽科昆蟲為外群，分析整個(gè)葉蟬科的DNA序列，分析獲得角頂葉蟬亞科是分析凹緣菱紋葉蟬的合適外群。以此構(gòu)建的進(jìn)化樹各節(jié)點(diǎn)的自展支持值均大于90%，且與網(wǎng)絡(luò)圖一致，說明該外群的選擇是合適的。

基因流是種群間基因通過個(gè)體或其配子遷移到其他種群而導(dǎo)致等位基因頻率發(fā)生變化的現(xiàn)象，是影響種群間和種群內(nèi)遺傳變異的主要因素[26]。本研究種群內(nèi)的變異遠(yuǎn)大于種群間，說明種群間的基因流增加了種群間的相似度和種群內(nèi)的變異度，也是陜南桑園凹緣菱紋葉蟬沒有形成地理格局的原因。基因流和地理距離之間的相關(guān)性分析顯示，該基因區(qū)域各個(gè)種群之間的基因流動(dòng)和地理距離不相關(guān)，出現(xiàn)這種狀況的原因，一方面可能與凹緣菱紋葉蟬具有一定的飛行能力和遷飛的習(xí)性有關(guān)[27]，成蟲和若蟲能夠適應(yīng)各種桑園生態(tài)環(huán)境，能夠漸次或連續(xù)發(fā)生一些近距離的擴(kuò)散，導(dǎo)致不同地理種群間發(fā)生遺傳物質(zhì)的交流。另一方面可能是人類生產(chǎn)活動(dòng)比如桑苗等宿主植物遠(yuǎn)距離運(yùn)輸?shù)鹊淖饔谩?/p>