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

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

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

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

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

1 材料與方法

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

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

1.2 線粒體基因擴增和單倍型網絡關系圖構建

迷卡斗蟋線粒體COI基因引物序列參考Simon等(1994)，LCO: 5′-GGTCAACAAATCATAAAGATAT TGG-3′；HCO: 5′-TAAACTTTCAGGGTGACCAAAAA ATCA-3′，由上海鉑尚生物公司合成，由上海生工公司對擴增產物進行測序，擴增長度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, 日本)檢測測序結果的峰圖并拼接序列，對單峰結果的序列進行拼接編輯，對于具有雙峰的序列的個體重新擴增測序。利用MEGA V.7.0.14(Kumaretal., 2016)中的Clustal W進行序列比對。利用DAMBE V.5.3.8(Xia and Xie, 2001)進行堿基替換飽和性檢測，結果顯示ISS

1.3 線粒體基因單倍型網絡關系圖構建

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

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

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

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

1.5 遺傳多樣性分析

在Genalex v6.5 (Peakall and Smouse, 2012)中計算種群各位點的等位基因數(number of alleles,Na)、有效等位基因數(number of effective alleles,Ne)、觀測雜合度(observed heterozygosity,Ho)、期望雜合度(expected heterozygosity,He)和近交系數(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)中計算遺傳分化指數(Fst)、種群內近交系數(Fis)和基因流(gene flow,Nm)。在Cervus v3.0(Kalinowskietal., 2007)中計算多態(tài)信息含量(polymorphism information content,PIC)。

1.6 遺傳分化分析

將篩選好的微衛(wèi)星位點應用于6個種群共134個個體中。為檢測種群間的分化程度并分析分歧來源，在Arlequin中計算遺傳分化指數(Fst)并進行AMOVA分析。Fst值是利用成對遺傳距離計算的，統計意義為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 遺傳結構分析

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

2 結果

2.1 基于線粒體基因COI的單倍型網絡關系

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

圖1 浙江省3個地區(qū)不同滯育型迷卡斗蟋6個種群基于mtCOI基因的單倍型中介網絡關系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)星位點的遺傳多樣性。結果顯示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)星引物檢測到每個位點的等位基因數(Na)為10.333(位點Vm37和Vm41)～17.833(位點Vm101)，平均13.803；觀測雜合度(Ho)為0.327(位點Vm37)～0.898(位點Vm101)，平均0.743。 位點Vm54(Fis=-0.012)起表現為存在雜合子冗余，Fis值0.006(位點Vm101)～0.631(位點Vm37)表示雜合子缺失(表3)。 等位基因數(Na)的平均值范圍為(6.000±2.966)(BSZ-E)～(18.000±2.449)(TT-E)個。 所有種群等位基因數(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分化程度最高，Fst值為0.21087，與其他種群的之間的Fst值為0.15783～0.19513。BSZ-N與其他種群的Fst值差異較大，Fst值為0.03242～0.05796(P<0.001)(表4)，所有種群之間的Nm值均大于1，為1.87111～36.81399(表4)，表明基因交流頻繁。

表4 浙江省3個地區(qū)不同滯育型迷卡斗蟋6個種群間基于11個微衛(wèi)星位點的基因流Nm值(上三角)和遺傳分化指數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分析結果顯示，總體遺傳變異均來源于個體之間?；诘乩砦恢梅譃?組(group 1-3)，組間差異占比-1.87%，組間種群內差異占比9.35%，組內種群間差異占比12.69%，個體間差異占比79.83%?；跍皖愋头譃槁褱M和若蟲滯育組兩組(group 4和5)，組間差異占比1.27%，組間種群內差異占比7.02%，組內種群間差異占比12.60%，個體間差異占比79.12%(表5)。因此種群間差異主要來源于個體之間，與地理隔離或者滯育型關系無顯著相關性。

表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 遺傳結構

使用Structure對浙江6個種群進行遺傳結構分析，基于基因型頻率分析得最佳分組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個種群的遺傳結構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 討論與結論

3.1 迷卡斗蟋遺傳多樣性

等位基因數(Na)、觀測雜合度(He)、期望雜合度(Ho)和多態(tài)信息含量(PIC)是反映種群遺傳多樣性的重要指標(Maudetretal., 2002)。本研究共成功篩選出迷卡斗蟋11個具有多態(tài)性的微衛(wèi)星位點，每個位點的平均等位基因數(Na)為13.803，平均觀測雜合度(He)為0.743，PIC值范圍為0.826～0.928，篩選的微衛(wèi)星具有較高的多態(tài)性(表3)。利用這11對微衛(wèi)星位點分析了浙江省天童山、大盤山和百山祖3個地區(qū)兩種滯育型共6個種群的遺傳多樣性。 每個種群的等位基因數為(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的自交系數(Fis)最低(表1)。

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

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

3.3 迷卡斗蟋遺傳結構

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

3.4 迷卡斗蟋生活史分化

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