基于線粒體基因組的云南脊扁蝽遺傳分化與譜系地理分析

摘要：[目的]分析中國云南脊扁蝽遺傳多樣性及系統發生關系，推測其起源及分鼓時間，并探討其歷史生物地理分布格局的成因。[方法]以2009-2020年采集于中國西南10個地區的云南脊扁蝽，采用高通量測序測定線粒體基因組，采用MEGA11，DnaSP6和Arlequin3.5.2.2等軟件分析其遺傳多樣性；以其他8個扁蝽科昆蟲作為外群，采用PhyIoSuite vl.2.3等軟件重建的云南脊扁蝽系統發育樹，并利用BEAST vl，10.4軟件采用寬松分子鐘和化石校正點推測估算云南脊扁蝽的分歧時間。[結果]本研究獲得30條長度為10 890 bp線粒體蛋白質編碼序列。共檢測到26個單倍型，單倍型多樣性Hd為0 989；核苷酸多態指數丌值為0.026 21，表現為高Hd、低丌的現象。AMOVA分析表明，云南脊扁蝽遺傳分化主要來自種群間。系統發育分析顯示中國西南地區云南脊扁蝽種群明顯分為4個支系，單倍型網絡圖與系統發育樹結果一致。分歧時間結果表明云南脊扁蝽在距今13.83 Ma的中新世中期開始分化，綠春種群（哀牢山脈南部）為中國地區的原始種群，后沿橫斷山脈分化出3個支系。[結論]第4紀冰期一間冰期輪回事件使中國西南地區云南脊扁蝽形成了現在的地理分化格局。

關鍵詞：云南脊扁蝽；線粒體基因組；單倍型；遺傳分化；譜系地理

中圖分類號：Q963 文獻標識碼：A 文章編號：1001-1498（2024）06-0045-09

云南脊扁蝽（Neuroctenusyunnanensis Hsiao）為半翅目（Hemiptera）扁蝽科（Aradidae）脊扁蝽屬（Neurovtenus）昆蟲，是蝽類中較為原始的種類。該物種身體扁平，常棲息于枯木及腐朽樹枝的樹皮下，若蟲長期在固定區域生活，成蟲對環境耐受性較強但壽命較短，飛行能力一般，無法遠距離遷徙，是理想的譜系地理研究物種。目前，對扁蝽的研究主要集中于傳統分類，系統發育等方面，關于扁蝽科昆蟲的譜系地理學研究較少，僅對蕭喙扁蝽（Brachyrhynchus hsiaoi Hsiao）、臺灣脊扁蝽（Neuroctenus taiwanicus Kromilev）、海南脊扁蝽（Neuroctenus hainanensis Liu）進行了相關研究，對其起源及分歧時間進行了分析和估計。而云南脊扁蝽目前的研究僅涉及傳統分類方面，但傳統的形態學分類不能準確的反映種群的譜系進化關系，需進行譜系地理學研究為云南脊扁蝽的種群遺傳學奠定基礎并幫助了解扁蝽科昆蟲的進化關系及演化歷史。

昆蟲的線粒體基因組通常是一類14～20 kb的雙鏈環狀DNA分子，因其序列較短、拷貝數高、母系遺傳、進化速率快、突變率高等特點，所以被廣泛應用于種群遺傳分化、系統發育和譜系地理學等方面的研究。且線粒體DNA在推斷種群或單倍型歷史方面比核基因更有優勢。但目前多數研究都是基于線粒體DNA單個基因片段，而單個線粒體基因片段具有一定的局限性，如在蝎蛉科Panorpidae研究中，coxl用于物種界定時區分度過低。與單個基因片段相比，線粒體全基因組在種群遺傳結構、系統發育和譜系地理學研究等方面更具優勢，得出的結果更加穩定可靠。

本研究對采集于中國西南地區10個地理種群的云南脊扁蝽線粒體基因組進行序列測定與對比分析?；跍y定的數據，對云南脊扁蝽的種群遺傳多樣性、系統發育、種群地理歷史遷移進行分析，并依據化石時間點和寬松分子鐘來估算種群分歧時間，再結合地質事件及氣候變化歷史資料，探討云南脊扁蝽在我國的起源及種群譜系地理模式。

1材料和方法

1.1材料

本研究于2009-2020年期間共收集中國西南地區10個地理種群云南脊扁蝽標本，具體信息見表1。為確保DNA的完整性，采集后需立即放人無水乙醇中，回到實驗室后置于-80℃冰箱中保存。并在每個地理種群隨機挑選3頭個體，編號為1～3，用于實驗分析。另外選取扁蝽科8個種的作為外群進行數據分析。

1.2DNA提取、測序

選取云南脊扁蝽的頭胸部和腿提取全基因組，使用天根生化科技有限公司的試劑盒提取DNA，參照使用手冊提取基因組DNA。將處理好后的樣品送至北京貝瑞和康生物技術有限公司進行測序，采用第2代重測序技術，測序模式為Novaseq6 000-S4-150PE。

1.3序列處理

得到測序數據后，以等脊扁蝽Neuroctenusparus （GenBank登錄號為NC_012 459）為參考基因組并利用Geneious 2023.2.1軟件進行拼接，得到云南脊扁蝽線粒體基因組。將拼接好的序列導出，用MEGA11軟件進行比對后將序列上傳至NCBI官網，利用BLAST功能進行近緣種多序列比對，保證其序列的準確性。再采用在線軟件MITOS2 （http：//mitos2.bioinf.uni-leipzig.de/index.Py）注釋線粒體DNA。

1.4遺傳多樣性分析

使用MEGA 11軟件計算基于線粒體蛋白質編碼基因的自裔位點（Singleton sites）、變異位點（Variable sites）、簡約位點（Parsimony information sites）和保守位點（Conserved sites）；使用DnaSP6軟件分析單倍型個數、單倍型多樣性指數、核苷酸多平均差異數等；用Arlequin veision 3.5.2.2軟件的AMOVA模塊分析估算種群間的遺傳分化指數（ F-statistics，Fst）及基因流參數（Nm），顯著性的檢驗通過1 000次的重復取樣實現；基于獲得的單倍型序列，使用PopART軟件用來構建種群的單倍型網絡圖。

1.5系統發育分析

基于獲得的10個地理種群的云南脊扁蝽的單倍型數據與NCBI上下載的8個扁蝽科昆蟲作為外群構建種群系統發育樹（表2）。使用Phylosuitev1.2.3軟件構建最大似然法（maximumlikelihood，ML）和貝葉斯法（Bayesian inference，Bl）系統發育樹。首先使用軟件中Modle Finder功能模塊查找構建引系統發育樹和ML系統發育樹的最佳堿基替換模型（GTR+I+G），基于該模型使用Mrbayes模塊構建引系統發育樹，手動選擇外群，MCMC世代數為10 000 000，采樣頻率1 000，兩次運行，4條MCMC鏈，BuminFraction值為0.25，當ASDSF值低于0.01時，引運行收斂，得到引系統發育樹。再使用IQ-tree模塊構建最大似然ML樹，手動選擇外群，“Bootstrap”選擇Ultrafast，“Num of Bootstrap”使用默認值5 000，運行得到ML系統發育樹。引和ML系統發育樹使用在線網站TVBOT（https：//www.chiplot.online/Nbot.html）可視化并進行美化。

1.6分歧時間估計

基于蛋白質編碼基因利用BEAST vl.10.4計算云南脊扁蝽中國西南地區種群單倍型分歧時間?；贕TR模型和寬松分子鐘模型，并采用3個化石標本作為校正點，MCMC鏈總代數設為20 000 000代，每1 000代取樣1次，舍棄樣本比例設為25％。運行結果使用Tracer vl.7.2查看，檢驗各參數是否收斂（ESSgt;200），并使用TreeAnnotato vl.10.4獲取MCC樹（最大可信進化分枝樹），最后使用FigTree vl.4.4查看MCC樹并對其美化。

2結果與分析

2.1遺傳多樣性分析

本研究的得到30個線粒體全基因組長度介于15 035-15 180 bp之間。注釋后利用MEGA11將13個蛋白質編碼基因進行手動對比，拼接為1個PCGs數據集，最終獲得30條長度為10 890bp的序列。共檢測到9 952個保守位點、938個變異位點，而變異位點中有826個簡約信息位點及112個自裔位點。T、C、A、G的含量分別為38.3％、16.3％、30.5％、14.9％。

本研究共檢測到26個單倍型，單倍型多樣性Hd為0.989，核苷酸多態性指數π值為0.026 21。其中易貢（1）單倍型數量最少，其次為通麥（2）和黃連山（2）。就核苷酸多態性而言，達德（0.027 43）和百花嶺（0.008 94）π值較大，其余地理種群π值均較小（表2）。而Hd和π值越大，群體的遺傳多樣性就越高，兩者通常具有一致性。但除達德、百花嶺、易貢外其余地區均呈現高Hd、低丌的現象，對于該現象有以下幾個可能原因：核苷酸多樣性的提高需要較長的時間，新的單倍型卻只需單個堿基變異；環境劇烈變化導致瓶頸效應，后又出現奠基者效應（少數個體迅速擴張成為新的種群，但核苷酸多樣性還停留在較低水平）。

AMOVA分析結果顯示，云南脊扁蝽各種群間Fst值在-0.244～0.999之間。云南地理種群中姚家坪與其他地理種群的Fst值均較大；黃連山僅與達德的Fst值為0.216，與其他地理種群的Fst值均在0.8以上；昔馬、那邦、銅壁關之間的Fst值均小于零，Nm值均為inf，表明這3個地理種群間基因交流密切，而百花嶺與昔馬、那邦、銅壁關的Fst值分別為0.138、0.213、0.134，對應Nm值分別3.107、1.837、3.227，表明其存在一定的基因交流；西藏地理種群中易貢與通麥Fst值為0.578，Nm值為0.363，而嘎美與易貢、通麥Fst值分別為0.000和0.125，其對應的Nm值為3.500和inf，這些數據表明嘎美可能為西藏地理種群的原始種群（表3）。對整個群體而言，變異來源主要是地理種群間，達到了81.14％（表4），表明各地理種群間有著較大的遺傳差異；種群內部變異差異較低，僅為18.86％，結合前面的分析，可能是經歷瓶頸效應群體擴張后長期的內部基因交流而導致的。

2.2系統發育分析

基于線粒體蛋白質編碼基因構建的引和ML系統發育樹拓撲結構相似，兩棵樹都將云南脊扁蝽分為4個支系，支系I包括了云南昔馬、那邦、銅壁關、百花嶺的全部個體以及達德的少量個體；支系Ⅱ包括西藏的通麥、嘎美、易貢3個地理種群的全部個體；支系Ⅲ僅為云南姚家坪；支系Ⅳ包括達德大部分個體和黃連山全部個體。單倍型網絡圖的分析結果與引和ML系統發育樹分析結果相似，都能夠明顯區分出4個支系（圖1，2，3）。

2.3分歧時間估計

分歧時間結果如圖4所示，赫氏刺緣扁蝽和同扁蝽與其他扁蝽科昆蟲在始新世晚期（37.54 Ma，95％置信區間：33.67～41.23 Ma）形成分化；赫氏尖胸霜扁蝽和無脈扁蝽屬在漸新世早期（33.48Ma，95％置信區間：30.12-37.11 Ma）形成分化；短喙扁蝽屬在漸新世晚期（24.35 Ma，95％置信區間：21.84～27.06 Ma）形成分化；本研究的云南脊扁蝽與近緣種等脊扁蝽在中新世中期（13.83Ma，95％置信區間：12.33～15.51 Ma）形成分化，其中支系Ⅳ于更新世中期（1.94 Ma，95％置信區間：1.68～2.2 Ma）形成分化，隨后支系Ⅲ在更新世中后期（1.41 Ma，95％置信區間：1.21～1.6 Ma）分化，最后支系Ⅱ與支系I也于更新世末期（1.16 Ma，95％置信區間：1-1.33 Ma）形成分化。

3討論

核苷酸多樣性（π）和單倍型多樣性（Hd）不僅是評價遺傳多樣性的重要指標還是分析生物種群發展歷史及其結構特征的主要依據。若Hd≥0.5，π≥0.005，表明種群較為穩定，歷史較長。本研究單倍型多樣性Hd為0.989，核苷酸多態性指數丌值為0. 026 21（表2），表明中國西南地區云南脊扁蝽具有較為漫長的演化歷史且較為穩定。依據Wright提出的理論： Fstlt;0.05，種群分化程度較低或無分化；0.050.25，高度分化。以及Kumar等認為：Nmlt;1時，種群分化程度高；14時，相當于同一個種群。支系I中的通壁關、昔馬、那邦種群之間Fst均小于0.05，Nm值均為inf，表明這3個地理種群可以認為是同一個種群，支系Ⅱ中的通麥和嘎美與上述情況一樣，也可以看作同一個種群，在另外兩個支系中未發現基因交流如此密切的種群。

系統發育關系及分歧時間估計結果表明云南脊扁蝽分化形成于13.83 Ma，中國地區云南脊扁蝽原始種群為綠春種群（支系Ⅳ）。綠春位于哀牢山脈南部，研究證明由于亞洲和印度板塊的碰撞擠壓，在17～26 Ma時期形成哀牢山脈。據此推測有部分綠春種群可沿哀牢山向西北方向及越南方向擴張。本研究數據未包含哀牢山脈其他區域的云南脊扁蝽，但白曉拴在景東（哀牢山脈中段）及范明瀾在越南老街?。òЮ紊矫}南端）有采集記錄，表明云南脊扁蝽在哀牢山脈存在。

分歧時間結果顯示瀘水姚家坪（橫斷山脈高黎貢山中下段）種群為橫斷山脈原始種群，從哀牢山脈北端至瀘水姚家坪一帶地勢都較為平坦，且白曉拴于永平（橫斷山脈怒山南端）有相關采集記錄，推測云南脊扁蝽可朝西北方向通過怒山南端遷移擴散至瀘水姚家坪。其后姚家坪種群沿高黎貢山向南、向北擴散至百花嶺（高黎貢山南端）及嘎美（伯舒拉嶺南端）地區。而第四紀發生了（-2.6Ma）冰期-間冰期輪回事件，導致該地區出現物種不斷地擴散、隔離以及種群間相互混雜的現象，并形成一種特殊的冰期演化方式，即物種通過垂直遷移來適應冰期。所以推測哀牢山脈的云南脊扁蝽不再向西北方向遷移，而橫斷山脈的種群通過垂直遷移的方式來避免冰期影響，兩個種群便逐漸形成分化，最后于1.94 Ma分化為兩個地理種群。隨著第四紀冰期進入希夏邦馬冰期（0.8～1.1Ma），該地區物種進人海拔較低的山谷作為生物避難所，橫斷山脈物種種群間的基因交流被進一步阻斷，推測姚家坪云南脊扁蝽種群因該原因在1.16 Ma分化形成支系I、Ⅱ。這與冰期發生時間相吻合。進入間冰期后，冰期避難所變為了間冰期物種種群擴張的物種泵（species pump），使種群再一次擴散。但由于環境劇烈變化導致的瓶頸效應，新發生種群由原始種群的部分個體大量繁衍而成（奠基者效應），導致遺傳多樣性下降，所以大部分子群體均呈現高Hd、低π的現象。在后續冰期-間冰期輪回事件中，盈江云南脊扁蝽各種群（支系I）可向低緯度遷移，與哀牢山的達德種群產生一定程度的基因交流，所以在盈江種群支系I中出現了達德種群的單倍型，具體遷移路徑見圖5。本研究中部分數據只有采集記錄而沒有分子數據，故云南脊扁蝽種群的地理歷史遷移還需進一步補充研究。

4結論

地理氣候變化是影響云南脊扁蝽譜系分布格局的主要因素；不同地理種群遷移擴張路徑與山脈走向以及地理隔離等存在明顯的相關性，加上第四紀氣候劇烈變化的影響，最終形成了云南脊扁蝽如今的譜系地理分布格局。