宋發(fā)軍，楊瑞霜，呂昕芮，劉祖懿，王紅瑩，孟艷艷

（中南民族大學(xué) 生命科學(xué)學(xué)院& 生物技術(shù)國(guó)家民委重點(diǎn)實(shí)驗(yàn)室& 武陵山區(qū)特色資源植物種質(zhì)保護(hù)與利用湖北省重點(diǎn)實(shí)驗(yàn)室，武漢 430074）

重樓是我國(guó)名貴中藥材，重樓皂苷是其主要活性成分［1］. 因其藥用價(jià)值和經(jīng)濟(jì)價(jià)值高，且重樓藥材市場(chǎng)的供需矛盾突出，當(dāng)前重樓屬多個(gè)物種資源已瀕臨枯竭. 同時(shí)，由于重樓的藥用部位生長(zhǎng)周期長(zhǎng)，種子深度休眠，加之重樓皂苷的生物合成途徑仍不甚清楚，完全依賴于植物中提取，因此加劇了重樓屬植物的資源危機(jī)［2］.為此，研究者們從分子角度對(duì)重樓屬植物展開了相關(guān)研究，以期為該屬植物資源的保護(hù)和可持續(xù)利用提供科學(xué)參考［3］. 然而現(xiàn)階段重樓屬植物中關(guān)于重樓皂苷合成、種子休眠與萌發(fā)等的基因信息仍較少.這嚴(yán)重阻礙了人們對(duì)上述問(wèn)題的解析. 故而，快速獲得重樓屬植物的大量基因信息并研究其功能特征，從而解析重樓皂苷的合成途徑、種子萌發(fā)等機(jī)制，為后續(xù)人工調(diào)控重樓皂苷的生產(chǎn)與合成奠定基礎(chǔ)，顯得至關(guān)重要.

近年來(lái)，隨著高通量測(cè)序技術(shù)的快速發(fā)展，轉(zhuǎn)錄組測(cè)序技術(shù)在獲取大量基因信息、新基因挖掘、基因家族鑒定、代謝途徑及系統(tǒng)進(jìn)化關(guān)系的分析等方面發(fā)揮著重要作用［4］.前人也利用轉(zhuǎn)錄組測(cè)序技術(shù)分別以重樓屬植物的藥用組織、不同萌發(fā)狀態(tài)的種子等為材料，展開了相關(guān)研究［5-8］.這些研究均為二代測(cè)序技術(shù)，且主要以滇重樓為對(duì)象，而對(duì)另一個(gè)重樓正品藥源植物——七葉一枝花的研究較少.

在諸多測(cè)序技術(shù)中，第三代PacBio 技術(shù)在轉(zhuǎn)錄組測(cè)序方面有著更好的應(yīng)用前景，它在植物優(yōu)良品種選育、分析次生代謝產(chǎn)物生物合成途徑等研究中發(fā)揮著重要作用，尤其適用于基因信息資源匱乏的藥用植物領(lǐng)域［9-11］. 這是因?yàn)樽鳛樽钚乱淮臏y(cè)序技術(shù)，其最大讀長(zhǎng)可達(dá)40 kb，且結(jié)果不受序列特異性影響，也不會(huì)產(chǎn)生由PCR 過(guò)程中堿基突變等原因造成的測(cè)序錯(cuò)誤［12］. 另外，利用PacBio 測(cè)序技術(shù)還可以檢測(cè)到多種DNA 修飾如5mC、5hmC、m6A、硫代磷酸等［12-14］. 上述諸多優(yōu)點(diǎn)，使得PacBio 技術(shù)在對(duì)無(wú)基因組信息物種的研究上有著極大的優(yōu)勢(shì)，該技術(shù)已經(jīng)被應(yīng)用于多個(gè)藥用植物的相關(guān)研究［10，15］.

本研究以中藥重樓正品藥源之一的七葉一枝花為研究材料，進(jìn)行了PacBio全長(zhǎng)轉(zhuǎn)錄組的測(cè)序工作，并對(duì)所得的轉(zhuǎn)錄本進(jìn)行分析了和注釋.本研究將為重樓屬植物的分子研究提供數(shù)據(jù)支撐和基因資源，同時(shí)為重樓皂苷合成相關(guān)基因的篩選及其種子休眠機(jī)制等方面的研究奠定基礎(chǔ).

1 材料與方法

1.1 實(shí)驗(yàn)材料

2016 年4 月在湖北省巴東縣采集的4 年生七葉一枝花的5 個(gè)組織（根狀莖、莖、葉、花、果莢）及2016 年秋季采收的種子.每種樣品采集3份，其中新鮮根狀莖洗凈后切成小塊立即進(jìn)行生物反應(yīng)滅活處理（液氮冷凍），其他部位直接進(jìn)行滅活處理，于-80 ℃凍存.

1.2 RNA提取、反轉(zhuǎn)錄及測(cè)序

分別提取根狀莖、莖、葉、花、果莢和種子6個(gè)部位的RNA［RNeasy Plus Mini Kit（#74134），Qiagen公司］，等量混合后，使用Clontech SMARTer PCR cDNA 試劑盒合成第一條鏈cDNA，再進(jìn)行PCR 優(yōu)化，并利用BluePippinTM將測(cè)序數(shù)據(jù)按照長(zhǎng)度分為1～2 kb、2～3 kb、3～6 kb、5～10 kb 四個(gè)庫(kù)，進(jìn)行擴(kuò)大PCR 反應(yīng)（若轉(zhuǎn)錄本大于3～6 kb 則需要利用BluePippinTM進(jìn)行大小篩選）后即可進(jìn)行轉(zhuǎn)錄組測(cè)序.轉(zhuǎn)錄組測(cè)序工作委托深圳華大基因公司完成.

1.3 Reads Of Insert序列及全長(zhǎng)轉(zhuǎn)錄本的獲取

Sequel 測(cè)序 平 臺(tái)共 包 含16 個(gè)SMRT cell，每 個(gè)SMRT cell包含100萬(wàn)個(gè)ZMW（zero-mode waveguides），每個(gè)ZMW 為一個(gè)測(cè)序單元.同一ZMW 中的所有subreads 來(lái)自同一個(gè)轉(zhuǎn)錄本，其中堿基出錯(cuò)率隨機(jī)，通過(guò)subreads間比對(duì)以提高堿基質(zhì)量，并利用SMRT軟件獲取Reads Of Insert.使用pbtranscript.py腳本識(shí)別全長(zhǎng)轉(zhuǎn)錄本序列，檢測(cè)這些序列是否包含5′端引物、3′末端引物和polyA 尾巴，區(qū)分全長(zhǎng)序列和非全長(zhǎng)序列.

1.4 高質(zhì)量一致性序列（isoform）的獲取和校正

同一孔中所有subreads 先進(jìn)行矯正得到環(huán)形一致性序列，然后孔與孔之間矯正得到高質(zhì)量一致性序列，即isoform.利用SMRT 進(jìn)行isoform 水平的聚類，使用非全長(zhǎng)序列對(duì)聚類的isoform 進(jìn)行Quiver 質(zhì)量校正.使用cd-hit分別對(duì)高質(zhì)量（QV＞99%）和低質(zhì)量（QV＜99%）isoform 進(jìn)行去冗余分析，計(jì)算覆蓋度并對(duì)其序列長(zhǎng)度及分布進(jìn)行統(tǒng)計(jì).

1.5 轉(zhuǎn)錄本序列的功能注釋、CDS預(yù)測(cè)及SSR檢測(cè)

利 用NR、NT、GO、COG、KEGG、SwissProt、InterPro七個(gè)數(shù)據(jù)庫(kù)對(duì)得到的轉(zhuǎn)錄本序列進(jìn)行功能注釋，并根據(jù)聚類到相同isoform 的序列數(shù)量，計(jì)算isoform 覆蓋度.隨后選取與數(shù)據(jù)庫(kù)匹配度較高的片段，進(jìn)行Blast注釋以預(yù)測(cè)CDS序列，而未能進(jìn)行功能注釋的片段則利用EST Scan來(lái)預(yù)測(cè)CDS序列.然后設(shè)計(jì)特定引物，通過(guò)PCR技術(shù)檢測(cè)七葉一枝花簡(jiǎn)單序列重復(fù)（Simple sequence repeat，SSR）的物種特異性.

2 結(jié)果與分析

2.1 測(cè)序結(jié)果及數(shù)據(jù)質(zhì)量控制

基于Pacific Biosience RS Ⅱ平臺(tái)，按照測(cè)序片段長(zhǎng)度構(gòu)建4 個(gè)庫(kù)進(jìn)行測(cè)序.如表1 所示，整合4 個(gè)庫(kù)的數(shù)據(jù)，共獲得的1219115492 個(gè)堿基分布于357362 個(gè)reads 中，QC（Quality control）值在86%～93%之間.由于PacBio 測(cè)序錯(cuò)誤分布隨機(jī)，一般情況下高質(zhì)量的文庫(kù)讀段更長(zhǎng)，本研究構(gòu)建的文庫(kù)中有較多3 kb 以上的片段，且短讀段的QC 值較高，由此可以看出此次轉(zhuǎn)錄組測(cè)序結(jié)果較好，可為后續(xù)的數(shù)據(jù)組裝提供很好的原始數(shù)據(jù).

表1 測(cè)序數(shù)據(jù)概覽Tab.1 Sequencing data overview

2.2 高質(zhì)量isoform的獲取及分析

整合4個(gè)庫(kù)共獲得58763個(gè)isoforms，去除冗余后進(jìn)行質(zhì)量分析，篩選出QV＞99%的高質(zhì)量isoforms共52537個(gè)（表2），所得序列的相關(guān)信息均已提交在NCBI Bioproject 數(shù)據(jù)庫(kù)中（https：//www.ncbi.nlm.nih.gov/bioproject/，BioSample accession：SAMN09762366）.

表2 七葉一枝花中isoform的數(shù)量Tab.2 Number of isoform in Paris polyphylla Smith var. chinensis

對(duì)52537 個(gè)高質(zhì)量isoforms 進(jìn)行分析（表3），平均長(zhǎng)度為2607 bp，這遠(yuǎn)大于一代測(cè)序（600 bp）和二代測(cè)序（200 bp）長(zhǎng)度.其N50 達(dá)到2998 bp，表明本次測(cè)序結(jié)果中序列拼接質(zhì)量較好.對(duì)高質(zhì)量isoforms的序列長(zhǎng)度進(jìn)行統(tǒng)計(jì)，可見這些isoforms 主要聚集于1800 bp和3000 bp（圖1）.

圖1 七葉一枝花中高質(zhì)量isoform長(zhǎng)度分布直方圖Fig.1 Histogram for high quality isoform length in Paris polyphylla Smith var. chinensis

表3 七葉一枝花中高質(zhì)量isoform數(shù)Tab.3 Number of high quality isoform in Paris polyphylla Smith var. chinensis

綜上所述，本研究獲得的高質(zhì)量isoform 普遍為較長(zhǎng)的序列，說(shuō)明文庫(kù)質(zhì)量較好，而短序列質(zhì)量相較長(zhǎng)序列高，進(jìn)一步從側(cè)面反映出本次構(gòu)建isoform文庫(kù)整體質(zhì)量?jī)?yōu)良，可為后續(xù)基因篩選及相關(guān)分析提供更好的參考.

2.3 轉(zhuǎn)錄本的功能注釋及覆蓋度分析

使 用NR、NT、GO、COG、KEGG、SwissProt、InterPro 共7 個(gè)數(shù)據(jù)庫(kù)，對(duì)所得轉(zhuǎn)錄本進(jìn)行功能注釋，在52537 個(gè)高質(zhì)量isoforms 中總計(jì)有40725 個(gè)被成功注釋，涵蓋了七葉一枝花所有高質(zhì)量isoforms的77.52%.由圖2 可知，在5 個(gè)數(shù)據(jù)庫(kù)中共同檢索到的isoforms 有11901 個(gè)，在NR 中單獨(dú)檢索到3878 個(gè)isoforms，在KEGG 中單獨(dú)檢索到52 個(gè)isoforms，在COG中單獨(dú)檢索到的isoforms有5個(gè)，在SwissProt中單獨(dú)檢索到的為38 個(gè)，在InterPro 中單獨(dú)檢索到isoforms139個(gè).

圖2 七葉一枝花中isoform注釋結(jié)果的韋恩圖Fig.2 Venn diagrams of isoform annotation results in Paris polyphylla Smith var. chinensis

其中，COG 功能注釋的結(jié)果顯示（表4），參與復(fù)制、重組和修復(fù)的isoforms 最多，達(dá)5094 個(gè)；屬于一般功能預(yù)測(cè)的有3815 個(gè).此外，還有少部分isoforms分布于防御機(jī)制、胞外機(jī)構(gòu)等途徑.與此同時(shí)，還發(fā)現(xiàn)了653 個(gè)分布于次級(jí)代謝產(chǎn)物的生物合成、轉(zhuǎn)運(yùn)與代謝等途徑的isoforms，這些isoforms 極有可能與重樓皂苷合成途徑相關(guān).

表4 七葉一枝花中isoform的COG注釋分布Tab.4 COG annotation distribution of isoform in Paris polyphylla Smith var. chinensis

GO 注釋中的isoforms 共分為生物學(xué)過(guò)程（藍(lán)色）、細(xì)胞組成（綠色）和分子功能（紅色）3 大類（圖3），49個(gè)小類.在生物學(xué)過(guò)程這一大分類分布于代謝過(guò)程的isoforms達(dá)2724個(gè)，居首位；在細(xì)胞組成中分布于細(xì)胞和細(xì)胞部分的isoforms 最多，均為2367 個(gè)；而在分子功能這一大分類中，大多數(shù)isoforms被歸為催化活性一類，達(dá)2606個(gè).

圖3 七葉一枝花中isoform的GO注釋分布圖Fig.3 GO annotation distribution map of isoform in Paris polyphylla Smith var. chinensis

KEGG 注釋中，52537 個(gè)isoforms 分布于碳水化合物的代謝、翻譯、折疊、分類和降解、脂質(zhì)代謝等分類中.為進(jìn)一步了解基因的生物學(xué)功能，對(duì)注釋到的isoforms 進(jìn)行代謝路徑分析，結(jié)果見表5，共有31934 個(gè)isoforms 分屬于134 個(gè)代謝通路.其中，前3個(gè)代謝通路分別為代謝通路、次生代謝物的生物合成、剪接體，其注釋到的isoforms 數(shù)分別為7412、4631、2235 個(gè).此外，ABC 轉(zhuǎn)運(yùn)體、RNA 轉(zhuǎn)運(yùn)等代謝途徑中，也有大量的isoforms富集.

表5 七葉一枝花KEGG分析中部分代謝通路Tab. 5 Some metabolic pathways in KEGG analysis of Paris polyphylla Smith var. chinensis

2.4 CDS預(yù)測(cè)

為進(jìn)一步了解所得isoform 的基因信息，利用Blast 和EST Scan 注釋預(yù)測(cè)其CDS 序列，結(jié)果如表6所示.本研究最終獲得39343條CDS序列，平均長(zhǎng)度759 bp.對(duì)預(yù)測(cè)的CDS 序列長(zhǎng)度進(jìn)行整合，結(jié)果如圖4. 這些CDS 序列長(zhǎng)度分布在200 nt 到3000 nt 之間，其中長(zhǎng)度≥500 nt 及≥1000 nt 的序列占總序列數(shù)的比例分別為70.45%和30.82%，這表明本研究的測(cè)序結(jié)果較為準(zhǔn)確可靠，并且所獲得的長(zhǎng)片段isoform 可直接用于基因的全長(zhǎng)克隆及功能分析等研究.

圖4 七葉一枝花CDS序列長(zhǎng)度分布Fig.4 CDS sequence length distribution in Paris polyphylla Smith var. chinensis

表6 七葉一枝花中CDS預(yù)測(cè)結(jié)果Tab.6 CDS forecast results in Paris polyphylla Smith var. chinensis

2.5 SSR檢測(cè)

轉(zhuǎn)錄本聚集后，繼續(xù)對(duì)所得序列進(jìn)行SSR 分析（圖5）.SSR 檢測(cè)中雙核苷酸重復(fù)的含量最高，達(dá)到9925 個(gè)，占全部SSR 位點(diǎn)的46.30%；其次是單核苷酸重復(fù)，為8360 個(gè).此外三核苷酸重復(fù)為2498 個(gè)，占總SSR 位點(diǎn)的37.54%.該項(xiàng)結(jié)果為后續(xù)七葉一枝花及重樓屬植物分子標(biāo)記的開發(fā)提供了依據(jù).

圖5 七葉一枝花SSR檢測(cè)結(jié)果Fig.5 SSR test results in Paris polyphylla Smith var. chinensis

3 討論

通過(guò)轉(zhuǎn)錄組測(cè)序技術(shù)可獲得大量的轉(zhuǎn)錄本信息，這對(duì)于七葉一枝花等沒有基因組信息的物種而言，是獲取其基因序列和功能信息的重要手段. 因此，近年來(lái)人們對(duì)重樓屬植物的轉(zhuǎn)錄組研究逐漸增多.

但是對(duì)比相關(guān)研究發(fā)現(xiàn)，大多數(shù)重樓屬植物的測(cè)序研究主要以云南重樓為對(duì)象，對(duì)七葉一枝花的研究較少，并且由于測(cè)序目的不同，選取的組織和采用的測(cè)序技術(shù)也有所差異.前人研究中選取的大多為根狀莖或種子等單一組織，重樓屬植物中果莢和花等組織的測(cè)序研究尚未見報(bào)道.本研究利用七葉一枝花的根狀莖、莖、葉、花、果莢及種子6個(gè)組織進(jìn)行混合測(cè)序，所得的基因信息系統(tǒng)和全面.并且，已有的研究均采用二代測(cè)序技術(shù)，雖然獲得的unigenes 數(shù)量較多，但大多為200～300 bp 左右的短片段，難以拼接和有效使用.如LING 等人為研究重樓種子的休眠機(jī)制，對(duì)云南重樓種皮和種子進(jìn)行Illumina 測(cè)序，共鑒定得到146671 個(gè)平均長(zhǎng)度為923 bp 的unigenes［16］.LIU 等對(duì)云南重樓4 年生根和8 年生根進(jìn)行NGS 測(cè)序，總共獲得87577 個(gè)平均長(zhǎng)度為614 bp 的unigenes［17］. 相比二代測(cè)序所得的短片段序列信息，本研究使用PacBio 測(cè)序獲取的isoform 大多數(shù)為1800～3000 bp 的片段，平均長(zhǎng)度為2607 bp.因此，本研究利用PacBio 測(cè)序技術(shù)所獲得的長(zhǎng)片段的isoforms，可直接用于基因全長(zhǎng)預(yù)測(cè)與克隆及功能鑒定等實(shí)驗(yàn)，有效減少了擴(kuò)增目標(biāo)基因全長(zhǎng)等工作量，并可為二代測(cè)序數(shù)據(jù)的拼接提供參考模板，為后續(xù)研究提供了極大的便利.

轉(zhuǎn)錄本注釋中，本研究所得的isoforms 有77.52%在NR，NT等數(shù)據(jù)庫(kù)中得到注釋，在LI等［18］的滇重樓基因組測(cè)序分析結(jié)果中，轉(zhuǎn)錄本長(zhǎng)度多為500～900 nt，有94.36%的序列得到注釋；兩者都在NR 數(shù)據(jù)庫(kù)中得到了最大占比的注釋. 另外，KEGG分析中，本研究所得的isoforms更多的注釋到代謝通路、次生代謝物的生物完成、剪接體等通路，而LI 的研究中，更多的isoforms被注釋到代謝通路、RNA 轉(zhuǎn)運(yùn)、內(nèi)質(zhì)網(wǎng)中蛋白質(zhì)加工等通路［18］.這可能是因?yàn)?，兩個(gè)研究工作中所選用的測(cè)序方法不同，因此獲得轉(zhuǎn)錄本的數(shù)量和類型都產(chǎn)生了差異，也可能跟實(shí)驗(yàn)材料不同有一定的關(guān)系. 此外，與LIU 等的研究結(jié)果［19-20］對(duì)比，雖然本研究獲取并注釋到的isoform 總量 略 少（40725/65535），但 注 釋到NR、SwissProt、KEGG 等數(shù)據(jù)庫(kù)的isoforms 數(shù)量更多，且KEGG 和GO 的注釋功能分布中各類別占比基本一致.綜上，本研究所得的isoforms 數(shù)量和注釋的質(zhì)量都較為可靠，為后期相關(guān)研究中的基因挖掘、篩選和功能驗(yàn)證提供了有效資源.

植物種子休眠的解除伴隨著胚乳中大量貯藏蛋白質(zhì)的降解，同時(shí)與DNA 復(fù)制、糖酵解、信號(hào)轉(zhuǎn)導(dǎo)、能量代謝等途徑有關(guān)［19-20］. 此外，重樓皂苷的合成途徑主要分為異戊烯焦磷酸的合成、甾體碳骨架的形成、甾體皂苷元的形成與修飾以及甾體皂苷元糖基化重樓皂苷等4部分［21-22］.本研究不僅獲得了大量與糖類代謝、信號(hào)轉(zhuǎn)導(dǎo)、能量轉(zhuǎn)換等與種子萌發(fā)有關(guān)的轉(zhuǎn)錄本，同時(shí)也富集到了大量參與次生代謝物的生物合成及降解相關(guān)的isoform，對(duì)上述isoforms進(jìn)行更加詳細(xì)的篩選和分析，將為深入研究重樓皂苷合成相關(guān)基因以及重樓種子萌發(fā)等問(wèn)題奠定良好的基礎(chǔ).

由圖5 可知，SSR 檢測(cè)結(jié)果中雙核苷酸重復(fù)占據(jù)優(yōu)勢(shì)，其次為單核苷酸重復(fù)，而三苷酸重復(fù)較少，這與張成才、李俊仁等人的研究結(jié)果有所差異［8］，這可能是由于物種差異性和實(shí)驗(yàn)樣本的不同而造成的. 這些差異的SSR 位點(diǎn)，將進(jìn)一步為重樓屬植物的遺傳多樣性分析、目標(biāo)基因標(biāo)定、遺傳圖譜構(gòu)建和分子標(biāo)記輔助育種等提供基礎(chǔ). 另外，由于PacBio 測(cè)序技術(shù)獲得的轉(zhuǎn)錄本長(zhǎng)度較長(zhǎng)，更接近真實(shí)的轉(zhuǎn)錄本分布.本研究所預(yù)測(cè)到的CDS 長(zhǎng)度分布能較為真實(shí)地顯示七葉一枝花中CDS 的長(zhǎng)度分布情況，這為后續(xù)重樓屬植物的基因研究提供了重要參考.

綜上，本研究利用PacBio 測(cè)序技術(shù)獲取的七葉一枝花全長(zhǎng)轉(zhuǎn)錄組測(cè)序的數(shù)據(jù)質(zhì)量較為可靠和全面，該結(jié)果進(jìn)一步豐富了七葉一枝花乃至整個(gè)重樓屬植物的基因資源，為重樓屬植物的二代測(cè)序數(shù)據(jù)拼接提供了良好的參考模板，同時(shí)也為解析重樓皂苷合成途徑、種子萌發(fā)機(jī)制等問(wèn)題提供了良好的依據(jù).