扁稈荊三棱不同組織混合樣品的轉錄組分析

寧 宇, 鄭彥超, 武高潔, 王義飛,4, 張曼胤

(1.中國林業科學研究院 濕地研究所，北京 100091；2.四川若爾蓋高寒濕地生態系統國家定位觀測研究站，四川 若爾蓋 624500； 3.國家林業局華東林業調查規劃設計院，浙江 杭州 310019；4.濕地生態功能與恢復北京市重點實驗室，北京 100091； 5.河北衡水湖濕地生態系統國家定位觀測研究站，河北 衡水 053000)

扁稈荊三棱(Bolboschoenusplaniculmis)是河漫灘、湖岸、稻田等濕地類環境中的常見雜草。該物種在自然狀態下能通過地下根莖和球莖等結構不斷進行營養繁殖而產生大量無性系分株，這種過程被稱為克隆生長[1-2]。克隆生長賦予該物種生理整合、內部分工等高效利用資源和適應環境的能力，使其經常在自然群落中占據絕對優勢地位[2-3]。已有的研究通常將其視為水稻田中的惡性雜草而研究其防除技術[4-6]。錢希[7]建立了預測扁稈荊三棱再生苗和球莖生長量的線性回歸方程，指出機械和人工移除雖然暫時抑制生長，卻促進了再生苗的蔓延。康學耕和唐恩全等[8]認為其無性繁殖平均10 d一代，并且當年生的種子當年不能萌發，至次年6月中旬萌發出土形成實生苗，形成“即時型”和“延后型”相結合的爆發，對稻田生產極具威脅。該植物雖然在稻田中通常被視為惡性雜草，但其在不同類型的系統中卻可能具有多種應用潛力。在自然濕地中，其球莖由于淀粉含量較高，成為國家一級保護動物白鶴(Grusleucogeranus)在我國東北地區濕地停歇時的主要植物性食物[9]，并且有研究表明它可能是一種歷史悠久的傳統食物[10]，具有食品加工方面的應用潛力。該物種還具有較好的去污染及凈化水體能力，可以作為人工濕地的構建材料[11-12]。但是，目前的工作多屬于經典生物學和實驗生態學范疇，在分子水平的信息比較匱乏，不利于生物學機理的深入研究；因此，亟待補充相關數據，為農林相關工作提供科學依據。

近年來，基于高通量方法的轉錄組測序正成為植物生理生態研究領域的有力工具。轉錄組測序可以獲得某一物種特定組織或器官在某一狀態下的幾乎所有轉錄本序列信息，進而得到的基因功能注釋、蛋白質編碼區序列、基因的表達量、代謝途徑等大量信息，為進一步研究提供基礎數據和重要參考[13-14]。該方法不僅能夠用于有參考基因組序列的物種研究，如水稻[15-16]、葡萄[17]等，也能用于無參考基因組序列的物種[18-20]，并且在探討基因的差異性表達(Differential expression)和共表達(Co-expression)情況等方面具有較強功效。已有的研究顯示，扁桿荊三棱的地下球莖網絡系統在不同環境下具有很強的可塑性，根莖的長度、球莖的數量和營養成分都可能存在差異[3, 21-22]，從而造成不同的種群擴張策略。對該物種進行轉錄組分析，有助于獲取其生長調節相關的基因信息，并為進一步深入研究不同環境條件下的表達模式提供基礎。

針對現有的研究中扁稈荊三棱轉錄組信息缺失的情況，本研究采用高通量測序技術對扁稈荊三棱的莖、葉、根莖和球莖的混合樣品進行轉錄組測序，結合生物信息分析對獲得的Unigene開展功能注釋、代謝通路分析等方面的研究，并探討其潛在的SSR位點情況，以期為扁稈荊三棱遺傳多樣性分析、功能基因挖掘和利用、分子輔助育種等方面的研究提供數據和理論支持。

1 材料和方法

1.1 試驗材料

1.1.1 采集地點 樣品采集于野外自然種群，具體位置為河北衡水湖濕地生態系統國家定位觀測研究站湖岸生態監測區域內(115°38′39.84″，37°39′10.08″)。采樣時間為2018年7月份，該時間內扁桿荊三棱多處于克隆生長的起始階段，生長迅速[22]，較為適合采集不同組織部位的樣品。采樣種群所在區域屬暖溫帶大陸季風氣候區，整個湖區的地形為淺碟狀洼地，湖岸為自然平地，土壤以潮濕的黏土為主要類型。保護區內全年平均降水量為506.3 mm，主要集中于7-8月份，全年日照平均總時數為2 571.7 h[23]。

1.1.2 取樣方法 選擇長勢均勻的扁稈荊三棱種群，在其范圍內隨機設置5個1 m×1 m的小樣方，每個小樣方內隨機選擇3個單株，分別采集其莖、葉、球莖和根莖部分，分開盛裝并做好標記。總共采集15個扁稈荊三棱的單株，每個單株的取樣量為5 g左右。各部分采樣后立即使用去離子水沖洗去除雜質，然后封裝放入足量的干冰環境中暫時保存，并在24 h內存于-80 ℃低溫條件下，以避免RNA降解。

1.2 試驗方法

1.2.1 RNA的提取、建庫與測序 分別取莖、葉、根莖和球莖1 g的材料，采用RNeasy Plant Mini Kits(Qiagen, Inc., Valencia, CA, USA)試劑盒提取總RNA，對總RNA樣品使用DNase(TaKaRa, Japan)純化后，從每類組織的產物中取10 μg等量混合組成RNA池，然后使用Oligo(dT)磁珠分離和富集poly(A)mRNA，向富集后的mRNA中加入Fragmentation Buffer 并進行打斷，再以片段化的mRNA為模板，用隨機引物(Random hexamers)合成第一鏈cDNA，加入緩沖液、dNTPs、RNaseH 和DNA polymerase 合成第二鏈cDNA。經過試劑盒純化回收、黏性末端修復后，在cDNA 的3′ 末端加上堿基“A”并連接接頭，用瓊脂糖凝膠電泳進行片段大小選擇，然后進行PCR 擴增得到測序文庫。構建好的文庫用Agilent 2100 Bioanalyzer質檢合格后，使用Illumina Solexa HiSeq 2000進行雙端測序(Paired-end)。測序所獲原始數據開放存儲于figshare(doi：10.6084/m9.figshare.8150792)

1.2.2 測序數據的處理與拼裝 測序得到的原始reads并不都是有效reads，在后續組裝和分析前，首先對Raw reads進行預處理，移除含poly-N、接頭和低質量的片段，獲得Clean reads。然后使用Trinity，采用Paired-end的拼接方法，將具有一定長度overlap的reads 連成更長的片段，得到不含未知堿基的組裝片段，之后利用TGICL軟件聚類去冗余延伸得到一套Unigene，并最終選取每個Loci下最長的轉錄本作為確定的Unigene。

1.2.3 Unigene序列的同源比對及功能分類 把拼接獲得的所有Unigene通過Blast比對(Evalue<1e-5)到主要蛋白數據庫中，以獲得其通路和功能等方面的注釋。采用的數據庫包括NR(非冗余蛋白數據庫)、Swiss-Prot(蛋白質序列數據庫)、GO(基因本體數據庫)、KEGG(京都基因與基因組百科全書)和KOG(蛋白質真核直系同源數據庫)等，得到跟給定Unigene具有最高序列相似性的蛋白(如有并列，取第1條)，從而得到試驗中所獲得的Unigene序列的功能注釋信息。

1.2.4 轉錄組中的SSR分析 使用MISA軟件(http://pgrc.ipk-gatersleben.de/misa/)搜索所獲得的Unigene上的潛在SSR位點。搜索條件為2個堿基的motif(基序)，需要至少6個重復，3個堿基的motif需要至少5個重復，4～6個堿基的motif至少需要4次重復。如果2個SSR序列的距離短于100 bp，就會被認為是同一個潛在SSR標記。此外，為排除polyA結構的干擾，剔除了單核苷酸重復的情況。

2 結果與分析

2.1 測序與Denovo組裝結果

基于Paired-end測序獲得了大量的樣本數據。為排除測序錯誤率的影響，使用滑動窗口法去除低質量片段，質量閾值20(錯誤率=1%)，窗口大小5 bp，長度閾值35 bp，并切除reads中含N的部分序列。最終得到了有效RNA-seq數據共0.32億條，GC平均含量為44.76%，共32 417 032條有效reads，平均長度95.77 bp。對比質控前數據，有效率為87.11%。測序數據質量控制前后對比情況見表1。

表1 測序數據質控前后對比Tab.1 Comparison of sequence data after quality control

注：表中reads數量為雙端測序的結果。

Note: reads amount in the table is from pair-end sequencing.

將樣本的有效reads合并進行de novo拼接，對拼接序列去重復，最終得到了79 098個長度大于200 bp的轉錄本，大小73 Mb。取每個Loci下最長的轉錄本作為Unigene，得到了59 788個Unigene，大小50 Mb，其最大長度、平均長度和N50分別為10 923，842,1 402 bp，長度大于500 bp和1 000 bp的Unigene數目分別為30 171和18 542。圖1展示了所有Unigene的長度分布情況。

2.2 Unigene的注釋

將所獲得的Unigene序列比對到公共數據庫(NR、 Swiss-Prot、 KEGG、Go、KOG)，進行序列相似性分析(Evalue

圖1 所獲得的Unigenes長度分布情況Fig.1 The length distribution of acquired Unigenes

2.3 扁稈荊三棱Unigene的KOG分類

將獲得的Unigene與KOG功能分類數據庫進行比對，結果表明，有13 334個Unigene 被注釋在25種功能分類中(圖2)。其中，歸類到信號傳導機制(Signal transduction mechanisms)的Unigene數量最多(2 836個，4.74%)，其次是蛋白質翻譯后修飾與轉運、分子伴侶(Posttranslational modification, protein turnover, chaperones)(2 465個，4.12%)及一般功能(General function prediction only)(2 142個，3.58%)。最小的類別是胞外結構(Extracellular structures)(33，<1‰)和細胞運動(Cell motility)(6，<1‰)。

A.RNA 加工與修飾；B.染色質結構與變化；C.能量產生與轉化；D.細胞周期調控與分裂，染色體重排；E.氨基酸運輸與代謝；F.核苷酸運輸與代謝；G.碳水化合物運輸與代謝；H.輔酶運輸與代謝；I.脂類運輸與代謝；J.翻譯，核糖體結構與生物合成；K.轉錄；L.復制、重組與修復；M.胞壁/膜生物發生；O.蛋白質翻譯后修飾與轉運、分子伴侶；P.無機離子運輸與代謝；Q.次生產物合成，運輸及代謝；R.一般功能；S.功能未知；T.信號傳導機制；U.胞內分泌與膜泡運輸；Z.細胞構架。

A.Nucleotide transport and metabolism; B.Chromatin structure and dynamics; C.Energy production and conversion; D.Cell cycle control, cell division, chromosome partitioning; E.Amino acid transport and metabolism; F.Nucleotide transport and metabolism; G. Carbohydrate transport and metabolism; H.Coenzyme transport and metabolism; I. Lipid transport and metabolism; J.Translation, ribosomal structure and biogenesis; K.Transcription; L. Replication, recombination and repair; M. Cell wall/membrane biogenesis; O.Posttranslational modification, protein turnover, chaperones; P.Inorganic ion transport and metabolism; Q.Secondary metabolites biosynthesis, transport and catabolism; R.General function prediction only; S.Function unknown; T.Signal transduction mechanism; U.Intracellular trafficking, secretion, and vesicular transport; Z. Cytoskeleton.

圖2 Uingenes的KOG分類情況
Fig.2 KOG classification of Unigenes

2.4 扁稈荊三棱Unigene的GO分類

對Unigene進行GO功能分類預測，并有24 670個Unigene被注釋上GO分類，并按照基因的分子功能，參與的生物過程和所處的細胞位置進行分類。如圖3所示，在生物過程(Biological process)分類中，主要聚集于代謝過程(Metabiolic process)(28.66%)和細胞過程(Cellular process)(23.97%)；在細胞組分(Cellular component)分類中，前三位依次是細胞(Cell)(16.00%)、細胞部分(Cell part)(16.00%)和細胞器(Organelle)(10.32%)；在分子功能(Molecular function)分類中，主要聚集于結合(Binding)(25.85%)和催化活性(Catalytic activity)(23.12%)。

圖3 Uingenes的GO分類Fig.3 GO classification of Unigenes

2.5 扁稈荊三棱Unigene的KEGG通路分析

扁稈荊三棱共有13 141個Unigene映射到KEGG數據庫中并得到注釋，獲得的參考代謝通路(Pathway)數目為291。其中包含Unigene最多的是代謝通路是核糖體(ko03010)，共有405條Unigene，其次是剪切體(ko03040)，共有360條Unigene。Pathway富集性分析的前10個Pathway數據列于表2。參與淀粉與糖代謝通路(ko00500)的Unigene有326條，共統計篩選出60條參與淀粉合成的Unigene，編碼9個關鍵酶(表3)，其中5個Unigene編碼α-糖苷酶；11個Unigene編碼β-呋喃果糖苷酶；13個Unigene編碼己糖激酶；6個Unigene編碼果糖激酶；4個Unigene編碼葡萄糖-6-磷酸異構酶；5個Unigene編碼葡萄糖磷酸變位酶；6個Unigene編碼葡萄糖-1-磷酸腺苷酰基轉移酶；6個Unigene編碼淀粉合成酶及4個Unigene編碼1，4-α-葡聚糖分支酶，另外參與植物激素信號轉導的Unigene為279條，參與根莖伸長生長的赤霉素受體GID1的Unigene有4條，與赤霉素(GA)生物合成過程中的關鍵酶赤霉素3-β-雙氧化酶(GA3ox)同源的Unigene有5條，編碼細胞分裂素合酶的Unigene有5條，共計14條Unigenes。

表2 扁稈荊三棱中包含Unigenes數目最多的10個代謝通路Tab.2 Top ten metabolic pathways involving B.planiculmis Unigenes

表3 扁稈荊三棱中地下莖伸長及球莖淀粉合成相關基因Tab.3 The genes and transcription factors related to sprouting and starch synthesis in B.planiculmis

2.6 扁稈荊三棱轉錄組中潛在SSR位點的分布

利用MISA軟件對扁稈荊三棱的59 788條Unigene 序列進行掃描，掃描的堿基容量為50 341 770。結果表明，12 823個Unigene 中含有潛在SSR位點，包含多于1個SSR位點的Unigene有2 531個，發生頻率(含有SSR的Unigenes 數量與總Unigenes 數量之比)為21.45%。在剔除單核苷酸重復類型后，總共有9 698個SSR位點，平均0.193個SSR/kbp，SSR的分布頻率(SSR的個數與總Unigenes 的數量比)為16.22%，其中二核苷酸motif(6 971個)和三核苷酸motif(2 332個)是主要的串聯序列類型，分別占SSR總數的71.88%和24.05%,兩者合計占所有SSR 位點的95.93%。所有潛在SSR位點的重復單元詳細分布情況見表4。

表4 扁稈荊三棱轉錄組中潛在SSR位點重復單元的分布Tab.4 Occurrence frequency of different microsatellites motifs in B.planiculmis transcriptome

3 結論與討論

扁稈荊三棱為稻田和濕地環境中的常見雜草，隸屬于莎草科(Cyperaceae)三棱草屬(Bolboshoenus)，屬于無參考基因組序列的物種，目前尚未發現扁桿荊三棱或其同屬物種的轉錄組研究先例。具體到本研究，在研究策略上，采用莖葉、球莖、根莖混合樣品進行轉錄組測序，有助于在節約試驗成本的基礎上發掘到更多的轉錄本[24-25]。在數據質量上，通過高通量測序的技術手段，首次對扁稈荊三棱的轉錄組信息進行獲取和分析。獲得了59 788個Unigene，其最大長度為10 923 bp、平均長度為842 bp，N50為1 402 bp(即50%的數據分布在長度不小于1 402 bp的Unigene上)，注釋比率達到58.91%。以KEGG的代謝途徑數據庫為依據，可鑒別出291個代謝通路，并篩選到與地下莖伸長生長相關的Unigene 9條，與地下球莖中淀粉合成相關的Unigene 60條。此外，還獲得了潛在的SSR位點9 698個，分布在12 823個Unigene中，發生頻率為21.45%，其中二堿基和三堿基的motif為主要類型。類比已有的克隆植物轉錄組研究，如毛竹(Phyllostachysedulis) (均長612 bp，注釋比率69.75%)[26]、巨龍竹(Dendrocalamussinicus)(均長723 bp，注釋比率69.22%)[27]，本研究的Unigene長度較為進步，有助于獲得較完整的基因表達信息。對于未獲得注釋的Unigene，可能是由于長度較短無法匹配，也有可能是非編碼序列或者是新的基因[28]，需要進一步的試驗予以驗證。在分析思路上，GO和KOG的功能分類對初步了解基因的功能起著重要作用，而KEGG數據庫中的參考pathway不僅可以推測基因的功能，還可以研究基因在不同代謝通路中所在位置及作用，三者相輔相成，成為新物種中發掘功能基因的重要手段[29-30]。本研究通過KEGG 數據庫中的pathway分析篩選與淀粉合成相關的基因，同時結合所篩選到的Unigene在相關數據庫中予以功能注釋，進一步確保了所獲得基因的可靠性。

本研究結果有助于理解扁桿荊三棱的地下根(球)莖系統的發展過程。對克隆植物矢竹(Pseudosasajaponica)的研究表明，地下根莖的生長發育受到激素的調節，其中赤霉素在啟動地下莖節間快速伸長上具有重要的作用，而細胞壁生長及細胞骨架組織等相關功能基因則在維持節間快速伸長生長上具有重要作用[31]，扁稈荊三棱中也檢測到了調控地下莖伸長生長的赤霉素代謝途徑相關基因，但尚未發現細胞生長的相關基因。這可能是由于該物種在7月份正處于克隆生長的起始階段，地下莖系統尚在形成過程中有關。該發現為防治該植物種群的過度擴張提供了新思路。在通路分析中，核糖體、剪接體和內質網上的蛋白加工居于前三位，表明可能在進行較大量的蛋白質合成工作，這與其地下根莖系統的快速擴張需求是相符合的[9, 22]。嘌呤代謝通路也位于前列，作為核苷酸合成、能量物質及一些初級產物如蔗糖、多糖及磷脂類物質的前體物質，嘌呤及其代謝途徑在許多生物學過程中具有重要的作用[32]。這也從另一方面印證了該植物地下莖系統的能量與物質需求。此外，本研究的數據也檢測到較多淀粉合成相關的基因及通路，與木薯(Manihotesculenta)[33]、馬鈴薯(Solanumtuberosum)[34]等主要以地下塊莖儲存淀粉的植物相一致。已有的研究認為扁稈荊三棱球莖中的淀粉是一種“儲備”行為，可支持其度過環境脅迫時期，并為球莖萌苗的快速生長提供能量[3, 22]，結果進一步表明這種儲備從地下系統擴張的早期就已經開始，為進一步挖掘該植物球莖的飼養價值或野生動物保護價值提供了支持。在研究展望方面，由于本研究采用的是不同組織材料等量混合后測序和組裝，雖然有助于得到更多的轉錄組數據，成本控制也較優秀，但是由于可變剪切等特異性機制[35]，本研究結果在解釋組織差異化表達方面具有局限性。在未來的研究中，應注重對該植物不同生長發育階段，或不同組織的轉錄組情況進行考察，具體探討其對水淹、鹽漬和種間競爭等生境中主要限制條件的響應機制。