基于微衛(wèi)星分子標記的浙江省不同滯育型迷卡斗蟋的遺傳分化分析

王佰秋, 沈初澤, 武 雪, 謝慧聰, 宗靖淞, 余哲媛,白 蕓, 李 愷,*, 何祝清,*

(1. 華東師范大學生命科學學院, 上海 200241; 2. 北京師范大學生命科學學院, 北京 100875)

迷卡斗蟋Velarifictorusmicado隸屬于直翅目(Orthoptera)蟋蟀科(Gryllidae)斗蟋屬Velarifictorus，從南至北廣泛分布于東亞地區(qū)(Wangetal., 2020)，也有記錄顯示該物種已經(jīng)成功入侵北美地區(qū)，并不斷擴散(Bowles, 2018)。為適應不同緯度的氣候條件，迷卡斗蟋演化出了兩種截然不同的滯育型:卵滯育型和若蟲滯育型。卵滯育型多分布于東亞的北部地區(qū)，以卵的形式過冬，春末孵化，秋季羽化成蟲，交配產(chǎn)卵后死亡，完成整個世代；若蟲滯育型多分布于東亞的南部地區(qū)，以若蟲形式過冬，夏季羽化成蟲，交配產(chǎn)卵后死亡，卵一般在夏末孵化，在冬季以大齡若蟲形態(tài)過冬，翌年夏季羽化成蟲，完成整個世代(Masaki and Walker, 1987; He and Takeda, 2013)。在中國的浙江省內(nèi)發(fā)現(xiàn)有兩種滯育型種群同域共存的現(xiàn)象(Wangetal., 2020)。線粒體基因分析結(jié)果顯示，卵滯育種群與若蟲滯育種群的線粒體基因已經(jīng)發(fā)生了分化，相似滯育型的種群具有相似的COI基因，而不同滯育型的種群COI基因之間存在一定差異(Wangetal., 2020)。但線粒體基因是母系遺傳，即在不存在突變的情況下后代的線粒體基因完全與母蟲相同。雖然兩種滯育型的迷卡斗蟋羽化成蟲的時間不一致，但可能存在若蟲滯育種群的成蟲存活時間較長，卵滯育種群較早羽化成蟲，兩個種群的個體直接交配從而產(chǎn)生后代的情況。兩個滯育型種群是否有交配，又關(guān)系到如何定義迷卡斗蟋這個物種。因為根據(jù)物種的定義，若兩個種群長期不存在自然交配，且某一特征(包括生活史)存在顯著差異，應將其劃分為兩個物種。因此，通過分子生物學技術(shù)對其進行更為深入的研究十分關(guān)鍵。

在蟋蟀總科中，已有不少基于微衛(wèi)星技術(shù)對種群結(jié)構(gòu)進行研究。例如Vanhala和Cottrell(2008)篩選出針蟋Nemobiussylvestris的25個多態(tài)性微衛(wèi)星位點，用以分析森林破碎化對物種基因交流和遺傳多樣性的影響；Birge等(2007)篩選出針蟋屬Allonemobius的10個多態(tài)性微衛(wèi)星位點，用于研究其屬內(nèi)雜交情況；Bretman等(2008)篩選出可共同用于蟋蟀屬Gryllusbimaculatus和G.campestris的微衛(wèi)星位點；Tinghitella等(2011)對濱海油葫蘆Teleogryllusoceanicus的7個微衛(wèi)星位點進行分析，發(fā)現(xiàn)自其引入至夏威夷后，遺傳多樣性降低; Gupta等(2020)基于家蟋蟀的全基因組開發(fā)了9個微衛(wèi)星位點，可用于種群結(jié)構(gòu)分析。

綜上所述，通過微衛(wèi)星能夠較好地研究屬內(nèi)以及種下的種群遺傳結(jié)構(gòu)。因此，本研究新開發(fā)獲得11個微衛(wèi)星位點，并據(jù)此對浙江省3個地區(qū)(天童山、大盤山和百山祖)迷卡斗蟋的卵滯育型和若蟲滯育型種群進行測序分析，研究其遺傳結(jié)構(gòu)和基因交流情況，探究浙江省同域分布不同滯育型迷卡斗蟋種群之間的遺傳分化的情況。

1 材料與方法

1.1 迷卡斗蟋樣本的獲取與基因組DNA的提取

本研究所用的迷卡斗蟋共6個種群134個標本分別采集于浙江省天童山、大盤山和百山祖國家級自然保護區(qū)，采集時間地點及種群名稱具體見表1，根據(jù)采集時間及若蟲成蟲形態(tài)判別其滯育型，即：9和10月份采集獲得成蟲及6月份采集獲得若蟲為卵滯育型；反之，9和10月份采集獲得若蟲及6月份采集獲得成蟲為若蟲滯育型。所有材料采集后首先保存在無水乙醇中，后于實驗室更換為新的無水乙醇并保存在-20℃冰箱中。利用AxyPrep Genomic DNA Miniprep Kit(AXYGEN)試劑盒對蟋蟀后腿肌肉組織提取基因組DNA，利用Thermo Scientific NanoDrop One超微量分光光度計測定總DNA濃度。

1.2 線粒體基因擴增和單倍型網(wǎng)絡關(guān)系圖構(gòu)建

迷卡斗蟋線粒體COI基因引物序列參考Simon等(1994)，LCO: 5′-GGTCAACAAATCATAAAGATAT TGG-3′；HCO: 5′-TAAACTTTCAGGGTGACCAAAAA ATCA-3′，由上海鉑尚生物公司合成，由上海生工公司對擴增產(chǎn)物進行測序，擴增長度658 bp。PCR反應體系(50 μL): 2×Taq PCR Mix 25 μL, ddH2O 19 μL, 正反向引物(10 μmol/L)各2 μL, 基因組DNA 2 μL。擴增程序: 94℃ 4 min; 94℃ 30 s, 45℃ 30 s, 72℃ 30 s, 共35個循環(huán)；最后72℃ 5 min。利用ATGC V.7.0.2(Genetyx Corporation, Tokyo, 日本)檢測測序結(jié)果的峰圖并拼接序列，對單峰結(jié)果的序列進行拼接編輯，對于具有雙峰的序列的個體重新擴增測序。利用MEGA V.7.0.14(Kumaretal., 2016)中的Clustal W進行序列比對。利用DAMBE V.5.3.8(Xia and Xie, 2001)進行堿基替換飽和性檢測，結(jié)果顯示ISS

1.3 線粒體基因單倍型網(wǎng)絡關(guān)系圖構(gòu)建

使用DnaSP V.5.10.01(Librado and Rozas, 2009)軟件將COI基因的序列fasta文件按地理種群分組，生成單倍型文件.nex和.arp文件。在ARLEQUIN中導入.arp文件，計算各種群的單倍型頻率。將單倍型文件統(tǒng)計為.txt文件。在PopART(Leigh and Bryant, 2015)中導入單倍型文件.nex和種群的單倍型頻率.txt文件。利用中介(median-joining, MJ)網(wǎng)絡方法構(gòu)建單倍型網(wǎng)絡關(guān)系圖以檢測不同組之間的系統(tǒng)發(fā)育關(guān)系。

1.4 微衛(wèi)星分子標記開發(fā)及篩選

提取得到迷卡斗蟋的基因組DNA后，由上海派森諾公司進行測序，采用Illumina MiSeq系統(tǒng)收集數(shù)據(jù)并產(chǎn)生標準的數(shù)據(jù)文件格式，數(shù)據(jù)經(jīng)過過濾后，共獲得1.31 Gb數(shù)據(jù)量。根據(jù)引物重復單位和多態(tài)性的類型選擇101對引物，其中包含16對二核甘酸重復模體的引物，36對三核甘酸重復模體的引物，44對四核甘酸重復模體的引物，2對五核甘酸重復模體的引物，3對六核甘酸重復模體的引物。選擇4頭迷卡斗蟋(包括2頭卵滯育個體，2頭若蟲滯育個體)的基因組DNA進行微衛(wèi)星位點的擴增，初步篩選出能夠擴增出清晰條帶并具有多態(tài)性的位點。

通過種群個體再次篩選初篩的引物以檢測位點的遺傳多樣性參數(shù)。采用熒光PCR擴增微衛(wèi)星基因，通過毛細管電泳檢測擴增產(chǎn)物的條帶大小。選擇浙江天童山若蟲滯育種群，對20頭個體進行微衛(wèi)星位點的再次篩選。由上海生工公司合成3種熒光標記的引物，分別為FAM, HEX和NED。引物信息見表2。采用96孔板擴增，不同熒光標記且大小相差50 bp以上的擴增產(chǎn)物混合送樣以節(jié)約成本，PCR產(chǎn)物在ABI 3730XL測序儀(Applied Biosystems)上進行基因分型，在GeneMapper v4.0(Applied Biosystems)上進行分析。

1.5 遺傳多樣性分析

在Genalex v6.5 (Peakall and Smouse, 2012)中計算種群各位點的等位基因數(shù)(number of alleles,Na)、有效等位基因數(shù)(number of effective alleles,Ne)、觀測雜合度(observed heterozygosity,Ho)、期望雜合度(expected heterozygosity,He)和近交系數(shù)(inbreeding coefficient,Fis)；在Arlequin(Excoffier and Lischer, 2010)中進行哈德-溫伯格平衡(Hardy-Weinberg equilibrium, HWE)和種群的連鎖不平衡(linkage disequilibrium, LD)檢測。 在Microchecker里檢測是否存在無效等位基因，在Arlequin v3.5.2.2(Excoffier and Lischer, 2010)中計算遺傳分化指數(shù)(Fst)、種群內(nèi)近交系數(shù)(Fis)和基因流(gene flow,Nm)。在Cervus v3.0(Kalinowskietal., 2007)中計算多態(tài)信息含量(polymorphism information content,PIC)。

1.6 遺傳分化分析

將篩選好的微衛(wèi)星位點應用于6個種群共134個個體中。為檢測種群間的分化程度并分析分歧來源，在Arlequin中計算遺傳分化指數(shù)(Fst)并進行AMOVA分析。Fst值是利用成對遺傳距離計算的，統(tǒng)計意義為10 000個置換。AMOVA分析時采用兩種分析方法對種群進行分組并分析遺傳分化來源，(1)按照地理位置劃分：百山祖種群BSZ-E和BSZ-N為一組(group 1)，大盤山種群DP-E和DP-N為一組(group 2)，天童山種群TT-E和TT-N為一組(group 3)；(2)按照滯育型劃分：BSZ-E, DP-E和TT-E為卵滯育組(group 4)，BSZ-N, DP-N和TT-N為若蟲滯育組(group 5)。

1.7 遺傳結(jié)構(gòu)分析

利用Structure v2.3.4(Pritchardetal., 2000)軟件進行種群結(jié)構(gòu)分析，在線 http:∥taylor0.biology.ucla.edu/struct_harvest/確定最佳的種群組數(shù)，即最佳K值；利用Clumpp v1.1(Jakobsson and Rosenberg, 2007)進行Clusters重復抽樣分析以減少結(jié)果的偶然性；利用Distruct v1.1圖形化顯示種群遺傳結(jié)構(gòu)。

2 結(jié)果

2.1 基于線粒體基因COI的單倍型網(wǎng)絡關(guān)系

結(jié)果顯示迷卡斗蟋134個個體共分為18個單倍型，6個種群顯著分為兩個世系，且同一滯育型的種群聚為一個支系，即卵滯育型種群BSZ-E, DP-E和TT-E聚為一支，若蟲滯育型種群BSZ-N, DP-N和TT-N聚為另一支(圖1)。

圖1 浙江省3個地區(qū)不同滯育型迷卡斗蟋6個種群基于mtCOI基因的單倍型中介網(wǎng)絡關(guān)系Fig. 1 Median-joining haplotype network based on mtCOIgene of six populations of Velarifictorus micado withdifferent diapause types from three areas in Zhejiang Province種群信息見表1。For population information, see Table 1. 下同The same below.

2.2 微衛(wèi)星位點篩選

在101對微衛(wèi)星引物中，共篩選出16對能夠擴增出具有多態(tài)性的清晰條帶，進一步檢測這些微衛(wèi)星位點的遺傳多樣性。結(jié)果顯示16個微衛(wèi)星位點中有4個位點(Vm11, Vm27, Vm92和Vm95)偏離哈德-溫伯格平衡，Microchecker檢測顯示位點Vm11和Vm92存在無效等位基因，這很可能是造成其偏離的原因。位點Vm59的遺傳多態(tài)性低，其觀測雜合度為0.250，期望雜合度為0.219。其余11個微衛(wèi)星位點遺傳多態(tài)性較高，可用于后續(xù)分析(表2)。

2.3 遺傳多樣性

在6個種群共134個個體中，利用11對微衛(wèi)星引物檢測到每個位點的等位基因數(shù)(Na)為10.333(位點Vm37和Vm41)～17.833(位點Vm101)，平均13.803；觀測雜合度(Ho)為0.327(位點Vm37)～0.898(位點Vm101)，平均0.743。 位點Vm54(Fis=-0.012)起表現(xiàn)為存在雜合子冗余，F(xiàn)is值0.006(位點Vm101)～0.631(位點Vm37)表示雜合子缺失(表3)。 等位基因數(shù)(Na)的平均值范圍為(6.000±2.966)(BSZ-E)～(18.000±2.449)(TT-E)個。 所有種群等位基因數(shù)(Na)平均為13.803。BSZ-E的觀測雜合度(Ho)最低，為0.514±0.292，TT-N的Ho值最高，為0.850±0.100，平均0.743。Fis值為0.06990(TT-N)～0.18849(TT-E)，表示雜合子缺失(表1)。

表2 浙江省天童山迷卡斗蟋若蟲滯育種群的16個微衛(wèi)星位點信息及遺傳多樣性Table 2 Information of sixteen microsatellite loci in nymphal diapause populations of Velarifictorus micadofrom Tiantongshan, Zhejiang Province and their genetic diversities

續(xù)表2 Table 2 continued

表3 浙江省3個地區(qū)不同滯育型迷卡斗蟋6個種群基于11個微衛(wèi)星位點的遺傳多樣性信息Table 3 Genetic diversity information of six populationsof Velarifictorus micado with different diapause typesfrom three areas in Zhejiang Province based on11 microsatellite loci

2.4 遺傳分化和基因流

迷卡斗蟋種群間的Fst值范圍為0.01340(DP-N與TT-N)～0.21087(BSZ-E與BSZ-N)，其中BSZ-E與其他種群Fst差異大(P<0.001)。BSZ-E與BSZ-N分化程度最高，F(xiàn)st值為0.21087，與其他種群的之間的Fst值為0.15783～0.19513。BSZ-N與其他種群的Fst值差異較大，F(xiàn)st值為0.03242～0.05796(P<0.001)(表4)，所有種群之間的Nm值均大于1，為1.87111～36.81399(表4)，表明基因交流頻繁。

表4 浙江省3個地區(qū)不同滯育型迷卡斗蟋6個種群間基于11個微衛(wèi)星位點的基因流Nm值(上三角)和遺傳分化指數(shù)Fst值(下三角)Table 4 Gene flow Nm values (above the diagonal) and pairwise genetic differentiation coefficient Fst values(below the diagonal) between six populations of Velarifictorus micado with different diapause typesfrom three areas in Zhejiang Province based on 11 microsatellite loci

基于地理位置和滯育型兩種分組方法的AMOVA分析結(jié)果顯示，總體遺傳變異均來源于個體之間。基于地理位置分為3組(group 1-3)，組間差異占比-1.87%，組間種群內(nèi)差異占比9.35%，組內(nèi)種群間差異占比12.69%，個體間差異占比79.83%?；跍皖愋头譃槁褱M和若蟲滯育組兩組(group 4和5)，組間差異占比1.27%，組間種群內(nèi)差異占比7.02%，組內(nèi)種群間差異占比12.60%，個體間差異占比79.12%(表5)。因此種群間差異主要來源于個體之間，與地理隔離或者滯育型關(guān)系無顯著相關(guān)性。

表5 浙江省3個地區(qū)不同滯育型迷卡斗蟋6個種群基于11個微衛(wèi)星位點的AMOVA分析Table 5 AMOVA of six populations of Velarifictorus micado with different diapause typesfrom three areas in Zhejiang Province based on 11 microsatellite loci

2.5 遺傳結(jié)構(gòu)

使用Structure對浙江6個種群進行遺傳結(jié)構(gòu)分析，基于基因型頻率分析得最佳分組K值為3，當K=3時，DeltaK值最大，即6個種群分為3組，主要劃分為BSZ-E一組，BSZ-N中有一半個體與DP-E, DP-N, TT-E和TT-N為一組，其余個體為一組(圖2: A)。當K=2時，主要劃分為BSZ-E一組，其余個體為一組(圖2: B)。

圖2 基于Structure v2.3.4分析當K=3(A)和K=2(B)時浙江省3個地區(qū)不同滯育型迷卡斗蟋6個種群的遺傳結(jié)構(gòu)Fig. 2 Genetic structure of six populations of Velarifictorus micado from three areas in Zhejiang Provinceanalyzed by Structure v2.3.4 with K=3 (A) and K=2 (B)

3 討論與結(jié)論

3.1 迷卡斗蟋遺傳多樣性

等位基因數(shù)(Na)、觀測雜合度(He)、期望雜合度(Ho)和多態(tài)信息含量(PIC)是反映種群遺傳多樣性的重要指標(Maudetretal., 2002)。本研究共成功篩選出迷卡斗蟋11個具有多態(tài)性的微衛(wèi)星位點，每個位點的平均等位基因數(shù)(Na)為13.803，平均觀測雜合度(He)為0.743，PIC值范圍為0.826～0.928，篩選的微衛(wèi)星具有較高的多態(tài)性(表3)。利用這11對微衛(wèi)星位點分析了浙江省天童山、大盤山和百山祖3個地區(qū)兩種滯育型共6個種群的遺傳多樣性。 每個種群的等位基因數(shù)為(6.000±2.966)～(18.000±2.449)個，其中BSZ-E的最低(表1)。所有種群的觀測雜合度(Ho)平均為0.743(表3)，其中BSZ-E的觀測雜合度最低，為0.514±0.292，由此可見BSZ-E的遺傳多樣性不高，可能是由于采集時僅局限于一個區(qū)域；TT-N的觀測雜合度(He)最高，為0.850±0.100，與此對應的TT-N的自交系數(shù)(Fis)最低(表1)。

3.2 迷卡斗蟋遺傳分化和基因交流

Wright(1978)指出Fst<0.05表示種群間分化程度較低。BSZ-E與其他種群的遺傳分化系數(shù)(Fst)的范圍為0.15783～0.21087，其次為BSZ-N與其他種群的，這表明百山祖的種群與其他地區(qū)的迷卡斗蟋種群遺傳距離大(表4)。Wright(1978)認為當Nm>1時，種群間存在基因交流，減緩了種群間的分化?；蛄鹘Y(jié)果表示，各種群之間存在著頻繁的基因交流，雖然其中BSZ-E與其他種群的基因交流程度較低，為1.87111～2.66805，但是其他種群之間的基因交流十分頻繁，其中DP-N與TT-N的基因交流最頻繁，基因流值為36.81399(表4)?？梢姷乩砦恢酶覍儆谕N滯育型之間的迷卡斗蟋基因交流最大。

3.3 迷卡斗蟋遺傳結(jié)構(gòu)

對迷卡斗蟋的遺傳結(jié)構(gòu)進行分析，當分為兩組時(K=2)，BSZ-E獨立為一組，其他種群聚為一組。當進一步分為3組時(K=3)，BSZ-N中的一半與其他種群聚為一組，剩下的一半單獨為一組(圖2)，可見，百山祖地區(qū)與其他地區(qū)的迷卡斗蟋之間也許由于地理距離緯度的原因基因交流減少。而AMOVA方差分析結(jié)果顯示，基于地理位置分為3組時，個體間差異占比最高(表5)，說明組與組之間、種群與種群之間并沒有顯著的遺傳分化，迷卡斗蟋間的遺傳差異主要來源于種群內(nèi)部個體之間的差異，而不是地理上的差異或者不同的生活周期。綜上研究結(jié)果表明，不同滯育型之間的種群之間仍存在著頻繁的基因交流。

3.4 迷卡斗蟋生活史分化

基于線粒體基因的種群地理學分析顯示，迷卡斗蟋的mtCOI基因共18個單倍型，分化為兩個單倍型組(圖1)，分別對應于卵滯育種群和若蟲滯育種群，這與Wang等(2020)的研究結(jié)果一致。然而利用微衛(wèi)星數(shù)據(jù)的種群遺傳學分析顯示，具有兩種不同滯育型的種群間并未呈現(xiàn)顯著的遺傳分化，同域分布的卵滯育和若蟲滯育迷卡斗蟋存在廣泛的基因交流(圖2)，這說明在同一地區(qū)兩種不同滯育型的種群的成蟲階段可能仍具有相當程度的重合。有一種可能性為若蟲滯育種群在夏季成蟲后，一直存活至秋季，與秋季剛成蟲的卵滯育種群成蟲交配。綜合前人線粒體基因分析的結(jié)果，我們推測調(diào)控迷卡斗蟋生活史分化的主效基因可能位于線粒體基因組上，或很大程度伴隨母系遺傳。類似的現(xiàn)象在部分其他的昆蟲類群中也有報道，如大猿葉蟲甲Colaphellusbowringi不同地理種群雜交后代的滯育率主要取決于其母本(Laietal., 2008)，紅頭毛蛾Zygaenatrifolii父本對后代滯育率完全沒有貢獻(Wipking and Kurtz, 2000)。所以不同滯育型的兩者“雜交”后，卵由母蟲在其發(fā)生季節(jié)產(chǎn)出，只有與母蟲保持一致的滯育型，才能繼續(xù)存活，因此mtDNA清晰地分化為兩支。假設“雜交”后代保持雄蟲的滯育型且能夠繼續(xù)存活并產(chǎn)生后代，在mtDNA結(jié)果中應該能發(fā)現(xiàn)交錯的現(xiàn)象，而非如今清晰地劃分為兩組。綜上所述，兩個種群存在頻繁的基因交流，但“雜交”后代只保持了母蟲的滯育型且繼續(xù)繁衍。這既解釋了為何同域不同滯育型種群之間存在頻繁的基因交流，而mtCOI基因卻能夠清晰地分為兩支。將來，要理清不同滯育型迷卡斗蟋的遺傳分化的調(diào)控機制，還要對同域分布種群進行進一步的GWAS分析。