扁稈荊三棱不同組織混合樣品的轉(zhuǎn)錄組分析

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

(1.中國林業(yè)科學(xué)研究院 濕地研究所，北京 100091；2.四川若爾蓋高寒濕地生態(tài)系統(tǒng)國家定位觀測研究站，四川 若爾蓋 624500； 3.國家林業(yè)局華東林業(yè)調(diào)查規(guī)劃設(shè)計(jì)院，浙江 杭州 310019；4.濕地生態(tài)功能與恢復(fù)北京市重點(diǎn)實(shí)驗(yàn)室，北京 100091； 5.河北衡水湖濕地生態(tài)系統(tǒng)國家定位觀測研究站，河北 衡水 053000)

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

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

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

1 材料和方法

1.1 試驗(yàn)材料

1.1.1 采集地點(diǎn) 樣品采集于野外自然種群，具體位置為河北衡水湖濕地生態(tài)系統(tǒng)國家定位觀測研究站湖岸生態(tài)監(jiān)測區(qū)域內(nèi)(115°38′39.84″，37°39′10.08″)。采樣時間為2018年7月份，該時間內(nèi)扁桿荊三棱多處于克隆生長的起始階段，生長迅速[22]，較為適合采集不同組織部位的樣品。采樣種群所在區(qū)域?qū)倥瘻貛Т箨懠撅L(fēng)氣候區(qū)，整個湖區(qū)的地形為淺碟狀洼地，湖岸為自然平地，土壤以潮濕的黏土為主要類型。保護(hù)區(qū)內(nèi)全年平均降水量為506.3 mm，主要集中于7-8月份，全年日照平均總時數(shù)為2 571.7 h[23]。

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

1.2 試驗(yàn)方法

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

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

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

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

2 結(jié)果與分析

2.1 測序與Denovo組裝結(jié)果

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

表1 測序數(shù)據(jù)質(zhì)控前后對比Tab.1 Comparison of sequence data after quality control

注：表中reads數(shù)量為雙端測序的結(jié)果。

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

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

2.2 Unigene的注釋

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

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

2.3 扁稈荊三棱Unigene的KOG分類

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

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

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

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

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

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

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

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

2.6 扁稈荊三棱轉(zhuǎn)錄組中潛在SSR位點(diǎn)的分布

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

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

3 結(jié)論與討論

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

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