薛清 杜虹銳 薛會英 王譯浩 王暄 李紅梅

（1. 南京農業(yè)大學植物保護學院 農作物生物災害綜合治理教育部重點實驗室，南京 210095；2. 西藏農牧學院資源與環(huán)境學院， 林芝 860000）

滑刃線蟲屬（AphelenchoidesFischer，1894）已報道種類約200種，大多數(shù)種類取食真菌，有13個種類為兼性植物寄生線蟲，具有廣泛的寄主范圍［1］。貝西滑刃線蟲（A. besseyi）、草莓滑刃線蟲（A.fragariae）和菊花滑刃線蟲（A. ritzemabosi）因其對作物產(chǎn)量和品質造成嚴重影響而受到廣泛的關注和研究，它們的典型寄主分別為水稻、草莓和菊花，并引起水稻干尖病、草莓夏矮病和菊葉線蟲病［2-4］。針對其他植物寄生性滑刃線蟲種類的研究則相對缺乏。

苜?；芯€蟲A. medicagus是新近描述的一種滑刃線蟲［5］。該線蟲分離自寧波口岸進境的美國干苜蓿草，可在鏈格孢菌（Alternariasp.）、灰葡萄孢菌（Botrytis cinerea）等多種真菌上完成生活史，同時對大豆和苜蓿具有弱寄生性［6］。系統(tǒng)發(fā)育研究顯示，苜蓿滑刃線蟲與貝西滑刃線蟲、福建滑刃線蟲（A. fujianensis）親緣關系較近［5］，但其隸屬的滑刃線蟲屬在墊刃亞目（Tylenchina）中的分類地位尚未完全解析，特別是與真滑刃線蟲屬（Aphelenchus）的高階元系統(tǒng)發(fā)育關系仍然存在較多的爭議［7-8］?，F(xiàn)有的線蟲系統(tǒng)發(fā)育研究大多基于核糖體RNA與線粒體COI基因，這些基因的片段一般較短（常用片段約為450-1 600 bp），所含信息量有限，難以完全解析各物種間的進化關系。動物線粒體基因組一般在10 kb以上，兼具編碼與非編碼區(qū)域，富含豐富的遺傳信息，其中基因結構排列、核苷酸和氨基酸序列，都可成為研究進化關系的良好靶標。因此，基于線粒體基因組的系統(tǒng)發(fā)育已成為近年來分類進化研究方向的熱點之一［9］。

線蟲線粒體基因組的研究相比其他動物類群起步較晚，目前已完成的線粒體基因組測序絕大多數(shù)集中于重要的動物寄生線蟲，而植物寄生類和土壤腐生類的線蟲鮮有報道［10］?；芯€蟲屬種類眾多，然而目前僅有貝西滑刃線蟲完成了線粒體全基因組測序［11］。造成這種現(xiàn)象的重要原因之一是現(xiàn)有的線蟲線粒體基因組測序大多采用操作難度較大的長PCR擴增法，該方法是根據(jù)已有的線粒體基因組信息，設計多對特異性引物進行長片段的擴增與測序［11-13］。然而，線蟲種間的線粒體基因差異較大，且可供參考的線蟲線粒體基因組較少，因此可能會導致引物設計存在較大盲目性，使實踐操作中成功擴增長片段變得困難。近年來，高通量測序技術的發(fā)展為線粒體基因組測序提供了新的策略，已經(jīng)在多種生物中獲得了較好的結果［14］，然而該方法尚未在線蟲的線粒體基因組研究中得到廣泛應用。目前采用高通量測序法僅完成了艾靈頓球孢囊線蟲（Globodera ellingtonae）和哥倫比亞紐帶線蟲（Hoplolaimus columbus）線粒體基因組測序與組 裝［15-16］，尚未見其他線蟲類群的國內外相關報道。

本研究利用二代測序技術Illumina平臺對苜?；芯€蟲進行低覆蓋全基因組測序，從獲得的原始數(shù)據(jù)中提取線粒體信息，通過多種方法組裝線粒體基因組后，對線粒體基因組的基因結構和排列、堿基組成，蛋白編碼基因和密碼子使用情況，以及RNA基因和非編碼區(qū)的結構等進行預測分析，結合已發(fā)表的線蟲線粒體基因組序列對苜?；芯€蟲的系統(tǒng)發(fā)育關系進行探討，旨在為滑刃線蟲的分類與線粒體基因組研究提供新的科學方法和依據(jù)。

1 材料與方法

1.1 材料

供試苜?；芯€蟲來自本實驗室在灰葡萄孢平板上保存的純培養(yǎng)，為該種描述時使用的模式群體［5］。挑取300條線蟲用無菌蒸餾水清洗后，放入裝有8 μL蒸餾水的200 μL PCR管中。采用液氮凍融法進行線蟲DNA的粗提，將PCR管在液氮中冷凍1 min，然后迅速放入65℃水浴鍋中溶解2-3 min，重復凍融步驟3次后，向PCR管中加入100 mg/mL蛋白酶K 1 μL和10×PCR緩沖液1 μL，混勻后于65℃保溫1.5 h，然后于95℃保持15 min。使用Qubit 1x ds DNA HS檢測試劑盒及 Qubit?3.0型熒光定量儀測定 DNA溶液的濃度。

1.2 方法

1.2.1 組裝種子序列的獲取 本研究采用線粒體COI基因部分序列作為種子，用于線粒體基因組組裝的起始序列。采用引物對COI-F1（5'-CCT ACT ATG ATT GGT GGT TTT GGT AAT TG-3'）與COI-V2（5'-GTA GCA GCA GTA AAA TAA GCA CG-3'）擴增COI基因片段，預測片段大小約為720 bp［17］。擴增體系為25 μL，含有DNA模板1 μL、濃度10 μmol/L的上游和下游引物各2 μL、100 μmol/L dNTPs 2 μL、10×ExTaq PCR緩沖液2.4 μL、25 mmol/L MgCl22 μL以 及5 U/μL ExTaq酶0.1 μL，加ddH2O補 足 至25 μL。PCR反應條件為：94℃ 5 min；94℃ 30 s，51℃ 30 s，72℃ 2 min，循 環(huán)42次；72℃延 長10 min。PCR產(chǎn)物經(jīng)凝膠電泳檢測后，由北京擎科生物科技有限公司采用Sanger法進行測序。

1.2.2 高通量測序與線粒體基因組組裝 采用標準的 Illumina TruSeq Nano DNA LT文庫制備實驗流程構建苜?；芯€蟲DNA樣本所需的基因組上機文庫，由南京派森諾基因科技有限公司使用 Hiseq X ten PE150平臺進行高通量測序，測序產(chǎn)生2 Gb數(shù)據(jù) 量，得到低覆蓋全基因組測序數(shù)據(jù)。采用FSATP［18］對下機數(shù)據(jù)進行質控，去除接頭污染，過濾掉長度低于50 bp、Q值低于20或n數(shù)量大于3的reads。

從GenBank中下載貝西滑刃線蟲（A. besseyi）、傷殘短體線蟲（Pratylenchus vulnus）、香蕉穿孔線蟲（Radopholus similis）、燕麥真滑刃線蟲（Aphelenchus avenae）、秀麗隱桿線蟲（Caenorhabditis elegans）、大豆孢囊線蟲（Heterodera glycines）以及花生根結線蟲（Meloidogyne arenaria）等7種線蟲的線粒體基因組序列作為參考基因組。在Linux系統(tǒng)下使用NextGenMap［19］根據(jù)線粒體基因組參考序列快速比對全基因組測序數(shù)據(jù)，比對閾值設定為0.3。通過SAMtools［20］提取比對結果中至少在單向符合要求的候選線粒體序列。利用上述獲得的COI基因擴增序列作為“種子”，結合自定義循環(huán)腳本與NOVOPlasty的種子序列擴展（seed-extend）算法，對已提取的線粒體基因組進行組裝。使用MitoZ［21］軟件并結合已發(fā)表的貝西滑刃線蟲［11］線粒體基因組信息，對獲得的苜蓿滑刃線蟲線粒體基因組進行注釋。根據(jù)基因組初步組裝結果，在nad4與cox1兩個基因上分別設計正向引物AM_F（5'- TTA TAC TCT ATT CGT TTT AGT GTC AGG TTT-3'）與反向引物AM_R（5'- TCA AAC AGC AAA ACC ACC TTG ATA CT-3'），以擴增非編碼片段。PCR采用 Prime STAR MAX酶，反應條件為：94℃ 5 min；94℃ 30 s，48℃ 30 s，62℃ 2 min，循 環(huán)38次；62℃延 長10 min。最后，將基因組組裝結果與擴增獲得的非編碼區(qū)進行拼接，將拼接后的基因組序列上傳至GenBank，登錄號為MZ424209。

1.2.3 基因注釋與系統(tǒng)發(fā)育研究 采用MITOS在線服 務（http：//mitos.bioinf.uni-leipzig.de/index.py）對組裝的線粒體基因組進行注釋［22］，并結合貝西滑刃線蟲線粒體基因組已有的注釋對其進行修正。利用MITOS軟件的MiTFi 程序對tRNA二級結構進行預 測［23］，并 用FORNA（http：//rna.tbi.univie.ac.at/forna）［24］軟件繪制tRNA二級結構。

為研究苜?；芯€蟲在線蟲門的系統(tǒng)發(fā)育地位，使用線蟲門31個不同種類的線粒體基因組作為參考序列，以科斯格羅夫異索蟲（Thaumamermis cosgrovei）作為外群構建基于極大似然法的系統(tǒng)發(fā)育樹。使用Geneious7.1.3軟件將12個蛋白編碼基因（PCGs）按照無脊椎動物線粒體密碼子分別轉換為氨基酸序列，并在MAFFT軟件中進行多重比對［25］。采用SMS軟件［26］對12個基因進行評估，選擇的模型為MtZoa+G+I+F。按照nad6-nad5-nad4-nad4Lnad3-nad2-nad1-cytb-cox3-cox2-cox1-atp6的順序構建序列矩陣，利用CIPRES運算平臺［27］的RAxML v8.2.12［28］軟件進行1 000次bootstrap運算，獲得極大似然自舉值（BS），構建極大似然系統(tǒng)發(fā)育樹。

2 結果

2.1 苜?；芯€蟲線粒體基因組的基因結構與 排列

苜?；芯€蟲DNA樣本經(jīng)Hiseq X ten PE150平臺測序共獲得213 159 776條原始序列，通過質控處理后得到211 293 116條高質量序列，占總原始序列的99.12%，測序總體質量較高。與7種線蟲的線粒體基因組序列進行比對和過濾后，共獲得線粒體序列1 430 492條，占獲得高質量序列的0.68%。拼接共獲得13 890 bp序列，為典型的閉合雙鏈DNA分子（圖1）。

圖1 苜?；芯€蟲線粒體基因組結構與基因排列順序Fig.1 Genome structure and gene arrangement in A. medicagus

通過線粒體基因組基因排列順序研究發(fā)現(xiàn)，苜?；芯€蟲與貝西滑刃線蟲、松材線蟲和擬松材線蟲的基因排列完全一致，而與燕麥真滑刃線蟲及根結線蟲、孢囊線蟲、短體線蟲等墊刃類線蟲存在較大的差異（圖2）。經(jīng)過注釋分別屬于36個基因，其中PCGs 12個，分別是atp6、cob、cox1-3、nad1-6和nad4L，無atp8基因；編碼線粒體核糖體RNA的2個rRNA（rrnS和rrnL）和22個tRNA基因。非編碼區(qū)位于nad4與cox1基因間，該部分未能通過測序片段完成組裝。采用結合非編碼區(qū)兩側區(qū)域的引物AM_F與AM_R進行擴增，并對獲得的片段進行測序，結果顯示非編碼區(qū)長度約1 600 bp，顯著少于貝西滑刃線蟲的3 142 bp。通過Sanger測序獲得其中1 393 bp序列，其余部分由于存在大量AT串聯(lián)重復，測序未能成功。對獲得的非編碼區(qū)進行分析顯示，除AT短重復外，該區(qū)域至少存在9個103 bp的串聯(lián)重復。將Sanger測序獲得的非編碼區(qū)與高通量測序的組裝結果進行拼接，共獲得14 411 bp的苜?；芯€蟲基因組，堿基A、T、C、G的含量分別為23.7%、55.9%、10.2%、10.2%，存在明顯的AT 偏好性。

圖2 墊刃亞目不同線蟲種群的線粒體基因排列順序Fig.2 Mitochondrial gene arrangement among different species in suborder Tylenchina

2.3 苜?；芯€蟲線粒體基因組的rRNA基因和tRNA基因結構

苜?；芯€蟲線粒體基因組中rrnS基因長度為910 bp，位于trnH和nad3基因間；rrnL基因長度為667 bp，位于trnE與trnS2之間，兩者AT含量分別為84.9%和83.2%。MUSCEL序列比對顯示，苜?；芯€蟲的rRNA與貝西滑刃線蟲的rrnL和rrnS的相似度最高，分別為72.1%和77.2%。對苜?；芯€蟲基因組進行同源性比對分析，預測得到了22個tRNA基因的序列，有關tRNA基因的位置見表1。22個tRNA基因的長度在53-57 bp 范圍之間，其中最長的為trnK，最短的為trnS1。

表1 苜?；芯€蟲的線粒體基因位置分布與起始終止密碼子Table 1 Placements and start/stop codons of mitochondrial genes in A. medicagus

對tRNA基因二級結構的預測顯示，得到的tRNA基因大部分均無法折疊形成典型的后生動物的三葉草結構，只有trnI基因序列的二級結構能折疊成典型的三葉草結構，其余21個tRNA 的二級結構均為非典型（圖3）。這些非典型的二級結構均屬于可變莖（variable arm）及配對臂缺失，呈TV-loop 結構。

圖3 苜?；芯€蟲線粒體基因組tRNA二級結構預測Fig.3 Predicted secondary structure of tRNA in the mitochondrial genome of A. medicagus

2.4 苜?；芯€蟲系統(tǒng)發(fā)育分析

將12個PCGs（atp6、cob、cox1-3、nad1-6和nad4）基因分別轉換為氨基酸序列，聯(lián)合構建了包含3 948個信息位點的序列矩陣。采用極大似然法構建了包括線蟲門兩個亞目20個科共30種線蟲的系統(tǒng)發(fā)育樹（圖4），從圖中可以看出，苜蓿滑刃線蟲與貝西滑刃線蟲互為姐妹群，且獲得高度支持（BS=100），且與松材線蟲、擬松材線蟲處在高度支持（BS=100）的滑刃線蟲科（Aphelenchoididae）單系中。燕麥真滑刃線蟲與滑刃科線蟲的親緣關系雖然較遠，但還是以較高的支持率與小桿亞目（Rhabditina）、旋尾亞目（Spirurina）、滑刃線蟲科、斯氏線蟲科（Steinernematidae）以及全凹線蟲科（Panagrolaimidae）的部分種類聚類在一個大分支中（BS=97）。以短體線蟲科（Pratylenchidae）、異皮線蟲科（Heteroderidae）與根結線蟲科（Meloidogynidae）為代表的墊刃類線蟲與寄生脊椎動物的盤尾絲蟲科（Onchocercidae）互為姐妹群，但拓撲結構支持率較低（BS=53），且兩者構成的分支在系統(tǒng)發(fā)育樹的地位尚無法解析（BS=29）。

圖4 基于12個蛋白編碼基因的氨基酸序列構建的極大似然系統(tǒng)發(fā)育樹Fig. 4 Phylogenetic tree based on amino acid sequences from concatenated 12 protein coding genes by maximum likelihood method

3 討論

線粒體基因組信息在物種多樣性與系統(tǒng)發(fā)育研究上有著廣泛的應用前景，但因傳統(tǒng)PCR擴增測序方法效率低，成本較高，目前僅有105個物種完成了線粒體基因組測序。相比超過25 000種已描述的線蟲［29］，大多數(shù)線蟲的線粒體基因信息仍然未知。已有研究顯示，線蟲的線粒體基因結構為一個雙鏈閉合的單環(huán)或雙環(huán)［10］，大小在12-22 kb，一般包括36個基因，即編碼氧化磷酸化所需酶的12個蛋白編碼基因（PCGs），分別是atp6、cob、cox1-3、nad1-6和nad4L，編碼線粒體核糖體RNA組分的2個rRNA（rrnS和rrnL），以及翻譯不同線粒體蛋白的22個tRNA。目前已知只有寄生脊椎動物的旋毛屬（Trichinellaspp.）和鞭蟲屬（Trichurisspp.）包含atp8基因，而大多數(shù)線蟲缺少atp8基因［12］。本研究共獲得14 411 bp的苜?；芯€粒體基因組序列，另有約1 600 bp序列因存在大量AT串聯(lián)重復，測序未能成功。對基因組成分析顯示，苜蓿滑刃線蟲的線粒體基因組與其他大多數(shù)線蟲較為相似，具有除atp8外的36個基因；在基因排列與組成上，與貝西滑刃線蟲、松材線蟲和擬松材線蟲這3種滑刃科線蟲的基因排序結構完全相同且序列相似度最高，而與燕麥真滑刃線蟲以及傷殘短體線蟲、香蕉穿孔線蟲、大豆孢囊線蟲和花生根結線蟲等墊刃類線蟲存在較大差異，顯示滑刃線蟲科與真滑刃線蟲科以及其他墊刃類線蟲均存在較遠的親緣關系，這與前期線粒體基因組系統(tǒng)發(fā)育研究［8，13］與核糖體RNA系統(tǒng)發(fā)育結果一致，支持滑刃科線蟲與墊刃類線蟲分別有著獨立的進化起源［30］。

對tRNA的比較分析顯示，苜蓿滑刃線蟲與貝西滑刃線蟲［11］最為相似，兩者具有相似的基因長度（分別為53-57 bp與51-59 bp），比哥倫比亞紐帶線蟲的基因長度變化范圍（52-72 bp）更窄。苜?；芯€蟲tRNA二級結構中有一個為典型的三葉草結構（rrnI），其余均為可變莖及配對臂缺失，貝西滑刃線蟲tRNA也均為可變莖及配對臂全部缺失，兩者較為相近，而哥倫比亞紐帶線蟲中有7個tRNA具三葉草結構［16］，明顯多于苜?；信c貝西滑刃線蟲線蟲，因此二級結構同樣支持滑刃科線蟲具有更為相近的親緣關系，而與墊刃類線蟲存在一定差異。

本研究中測序獲得的原始數(shù)據(jù)大部分為細胞核基因，僅有0.68%屬于線粒體基因，這與昆蟲中0.5%-1.4%的比例類似［31-32］，但是明顯高于哥倫比亞紐帶線蟲的0.2%線粒體序列含量［16］。截至到 2021年初，部分測序公司的建庫費用已降低至200-300元/文庫，單個樣品二代測序費用亦低至50元/Gb。因此，隨著測序成本的降低，通過增加測序通量彌補線粒體數(shù)據(jù)含量低，進而獲得充足的測序深度和組裝質量的策略在實踐中成為可能。

本研究首次完成了苜?；芯€蟲的線粒體基因測序，結果顯示，低覆蓋全基因組測序法可以成功組裝出線粒體基因組中大部分序列，證明了利用該方法獲取線蟲線粒體全基因組序列的可行性。同時，線粒體基因組與傳統(tǒng)的核糖體RNA相比擁有更多的信息位點，構建的系統(tǒng)發(fā)育樹也具有更高的支持率，因此低覆蓋全基因組測序法在線蟲線粒體基因組研究上將有廣闊的應用前景。

4 結論

通過低覆蓋全基因組測序法獲得了苜?；芯€蟲線粒體基因組，研究表明其基因的構成排列與貝西滑刃線蟲、松材線蟲及擬松材線蟲相同，系統(tǒng)發(fā)育地位相近。本研究證明了低覆蓋全基因組測序法獲取線蟲線粒體全基因組序列的可行性。