穿心蓮全長(zhǎng)轉(zhuǎn)錄組測(cè)序及特性分析

孫銘陽(yáng)，徐世強(qiáng)，顧 艷，梅 瑜，周 芳，李靜宇，王繼華

（廣東省農(nóng)業(yè)科學(xué)院作物研究所/廣東省農(nóng)作物遺傳改良重點(diǎn)實(shí)驗(yàn)室，廣州510640）

0 引言

植物次生代謝產(chǎn)物是生物活性物質(zhì)的主要來(lái)源，它們具有治療多種疾病的臨床應(yīng)用價(jià)值[1]。穿心蓮(Andrographis Paniculata)是一類重要的“寒性”中草藥，其地上部分入藥后可降溫消炎的記載歷史悠久[2]。穿心蓮次生代謝提取物中含有20多種結(jié)構(gòu)類似的二萜類化合物和10多種黃酮類化合物。其中，含量最多的二萜類物質(zhì)為穿心蓮內(nèi)酯，藥理活性豐富[3]。目前，新型冠狀病毒2019(2019 novel coronavirus，2019-nCoV)肺炎疫情已成為全球性的衛(wèi)生危機(jī)。因穿心蓮內(nèi)酯及其熒光衍生物可通過與主要蛋白酶功能位點(diǎn)形成共價(jià)鍵的方式抑制2019-nCoV的活性，現(xiàn)已被建議用于新冠肺炎的治療[4]。

穿心蓮內(nèi)酯合成相關(guān)的分子機(jī)理研究已取得部分進(jìn)展。穿心蓮內(nèi)酯可通過下調(diào)磷脂酰肌醇3激酶/蛋白激酶b信號(hào)通路基因來(lái)抑制腫瘤壞死因子-α誘發(fā)的炎癥[5]。13C核磁共振常用同位素富集光譜顯示，脫氧木酮糖途徑是二萜類前體被修飾為穿心蓮內(nèi)酯所經(jīng)歷的主要途徑[6]。隨高通量測(cè)序技術(shù)的興起，二代Illumina HiSeq短序列測(cè)序(Second Generation Sequencing,SGS)已成為分析基因表達(dá)水平和調(diào)控途徑的常用方法。Tong等[7]利用SGS數(shù)據(jù)揭示了激素對(duì)穿心蓮種子萌發(fā)過程的代謝及信號(hào)轉(zhuǎn)導(dǎo)進(jìn)程的影響。Cherukupalli等[8]將SGS輸出數(shù)據(jù)進(jìn)行二次組裝后，生成了大小為88 Mb的穿心蓮葉片轉(zhuǎn)錄組數(shù)據(jù)庫(kù)，共83800條序列片段，為其次生代謝途徑研究提供了初始平臺(tái)。然而，SGS技術(shù)的讀長(zhǎng)有限；組裝數(shù)據(jù)時(shí)需拼接；無(wú)法精準(zhǔn)評(píng)估可變剪切產(chǎn)生亞型的轉(zhuǎn)錄豐度及較長(zhǎng)轉(zhuǎn)錄本的全長(zhǎng)信息。

如今，以美國(guó)太平洋生物技術(shù)公司的PacBio平臺(tái)為主的第三代超長(zhǎng)單分子實(shí)時(shí)(Single Molecule Real Time,SMRT)測(cè)序技術(shù)已經(jīng)成熟。平均15 kb的超長(zhǎng)讀長(zhǎng)使數(shù)據(jù)無(wú)需拼接，可直接得到高質(zhì)量的轉(zhuǎn)錄本全長(zhǎng)序列和AS亞型信息[9]。為進(jìn)一步明晰穿心蓮次生代謝途徑基因的響應(yīng)情況，本研究利用PacBio平臺(tái)的三代SMRT測(cè)序技術(shù)對(duì)生長(zhǎng)60天的穿心蓮根、莖和倒三葉進(jìn)行全長(zhǎng)轉(zhuǎn)錄組測(cè)序和分析。結(jié)果表明，穿心蓮體內(nèi)功能基因主要參與次生代謝進(jìn)程；預(yù)測(cè)功能為轉(zhuǎn)錄因子(Transcription Factors,TF)的基因中，basic helixloop- helix(bHLH)、 myeloblastosis(MYB) 及WRKYGQK(WRKY)等響應(yīng)植物次生代謝的主要TF家族成員含量排名前10；合成二萜類前體的甲羥戊酸途徑(MVA)和2-C-甲基-D-赤蘚糖醇-4-磷酸途徑(MEP)基因發(fā)生的主要AS事件為內(nèi)含子保留；在搜索到的簡(jiǎn)單重復(fù)序列(Simple Sequence Repeat,SSR)中，AT/AT為優(yōu)勢(shì)重復(fù)單元。以上結(jié)果為穿心蓮次生代謝基因調(diào)控網(wǎng)絡(luò)研究、AS亞型的功能注釋以及穿心蓮遺傳多樣性的評(píng)價(jià)提供分子基礎(chǔ)。

1 材料與方法

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

選用福建漳州生產(chǎn)用穿心蓮種，種于草炭土(HAWITA)中，保持土壤濕潤(rùn)，于光照培養(yǎng)箱（光照day/night：16 h/8 h；溫度：28℃）中培養(yǎng)60天。將幼苗從土中移出，清洗根部后輕柔擦干，剪取根、莖和倒三葉并迅速混合置于液氮速凍，每個(gè)部位取3次獨(dú)立的生物學(xué)重復(fù)并混合，待RNA提取。穿心蓮生長(zhǎng)于廣東省農(nóng)業(yè)科學(xué)院作物研究所農(nóng)作物遺傳改良重點(diǎn)實(shí)驗(yàn)室，于2020年4—6月進(jìn)行。

1.2 全長(zhǎng)轉(zhuǎn)錄組文庫(kù)構(gòu)建

利用塔克拉生物科技公司的RNAiso Plus試劑提取樣品的總RNA；使用Nanodrop檢測(cè)樣品濃度；利用Agilent 2100檢測(cè)樣品片段完整度。質(zhì)量檢測(cè)合格的RNA樣品用于穿心蓮全長(zhǎng)轉(zhuǎn)錄組文庫(kù)構(gòu)建，建庫(kù)步驟參考Yang等[10]。

1.3 全長(zhǎng)轉(zhuǎn)錄組測(cè)序數(shù)據(jù)篩選

建庫(kù)完成后，上機(jī)（Pacbio三代測(cè)序儀）測(cè)序。利用SMRT Link v6.0軟件對(duì)輸出的原始序列進(jìn)行分析。將原始數(shù)據(jù)中的環(huán)形一致性序列(Circular consensus sequence，CCS)提取后，通過CCS對(duì)序列中是否存在測(cè)序引物序列或者是否存在嵌合序列進(jìn)行篩選、去冗余和校正[11]。最終得到高質(zhì)量全長(zhǎng)轉(zhuǎn)錄本序列，用于后續(xù)功能注釋等分析。

1.4 全長(zhǎng)轉(zhuǎn)錄組功能注釋

將轉(zhuǎn)錄本序列與National Center for Biotechnology Information(NCBI)網(wǎng)站中的blastx作比對(duì)，得到相似性最高的蛋白序列，用于進(jìn)行蛋白功能注釋。利用Eukaryotic ortholgous groups(KOG)進(jìn)行蛋白的直系同源分類；利用swiss prot protein database(SwissProt)數(shù)據(jù)庫(kù)進(jìn)行去冗余后的蛋白功能注釋；利用kyoto encyclopedia of genes and genomes(KEGG)數(shù)據(jù)庫(kù)進(jìn)行基因功能分類和代謝通路聚類；與nonredundant protein database(Nr)數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，取比對(duì)結(jié)果中E值最低（同源性最高）的序列所屬物種進(jìn)行各物種的同源序列數(shù)量統(tǒng)計(jì)[12]。

1.5 基因結(jié)構(gòu)分析

轉(zhuǎn)錄本在以上4個(gè)數(shù)據(jù)庫(kù)的蛋白比對(duì)結(jié)果中選取E值小于e-5，且rank最高的多肽序列對(duì)應(yīng)的堿基區(qū)域定為編碼區(qū)(Coding Sequence,CDS)。利用ANGEL軟件預(yù)測(cè)與以上4個(gè)數(shù)據(jù)庫(kù)比對(duì)不上的轉(zhuǎn)錄本序列的CDS區(qū)和蛋白序列[13]。對(duì)沒有注釋到以上4個(gè)數(shù)據(jù)庫(kù)的轉(zhuǎn)錄本序列進(jìn)行l(wèi)ncRNA分析，利用cnci軟件和CPC軟件進(jìn)行編碼能力的預(yù)測(cè)，取2個(gè)軟件均預(yù)測(cè)為“非編碼”的序列作為最終的lncRNA結(jié)果[14-15]。利用Signal P網(wǎng)站進(jìn)行蛋白信號(hào)肽預(yù)測(cè)。利用MISA軟件搜索SSR，配置參數(shù)信息為：1、2 nts至少重復(fù)6次；3 nts至少重復(fù)5次；4 nts至少重復(fù)4次；5 nts至少重復(fù)4次和6 nts至少重復(fù)4次才會(huì)被認(rèn)為是SSR。2、如果2個(gè)SSR的距離短于100 bp，就會(huì)被合并作為一個(gè)SSR。對(duì)于AS事件，本研究利用Cogent軟件組裝出經(jīng)聚類后轉(zhuǎn)錄本的CDS序列，并以組裝CDS序列作為參考，利用SUPPA軟件進(jìn)行可變剪切分析[16]。

2 結(jié)果與分析

2.1 穿心蓮全長(zhǎng)轉(zhuǎn)錄組測(cè)序的數(shù)據(jù)篩選

濃度合格的樣品總RNA經(jīng)雙鏈cDNA合成及加接頭，形成啞鈴型單分子雙鏈cDNA文庫(kù)后上機(jī)測(cè)序。通過提取下機(jī)序列中的CCS，共得到全長(zhǎng)序列487309條，其中全長(zhǎng)非嵌合(Full-Length Non-Chimeric,FLNC)序列354407條。FLNC序列在去掉5’引物、3’引物和polyA序列后即為物種原始全長(zhǎng)轉(zhuǎn)錄本序列。由圖1可知，穿心蓮幼苗的FLNC序列長(zhǎng)度主要集中在200～ 3200 nts之間。因基因組時(shí)刻處于轉(zhuǎn)錄mRNA的動(dòng)態(tài)中，所以測(cè)序得到的單分子FLNC序列存在大量冗余。將FLNC序列進(jìn)行聚類、矯正和去冗余后，最終得到76812條高質(zhì)量轉(zhuǎn)錄本信息，用于后續(xù)功能注釋和結(jié)構(gòu)分析。

圖1 FLNC序列長(zhǎng)度分布圖

2.2 穿心蓮全長(zhǎng)轉(zhuǎn)錄組的功能注釋

利用四大功能注釋數(shù)據(jù)庫(kù)KOG、Nr、Swiss Prot及KEGG對(duì)76812條穿心蓮全轉(zhuǎn)錄組信息進(jìn)行注釋，各數(shù)據(jù)庫(kù)比對(duì)情況見圖2。共有63855條序列得到了注釋，其中Nr數(shù)據(jù)庫(kù)注釋到的信息最多，共63822條序列；Swissport數(shù)據(jù)庫(kù)注釋到54430條序列；KOG數(shù)據(jù)庫(kù)注釋到43656條序列；KEGG數(shù)據(jù)庫(kù)注釋到30412條序列。由4個(gè)數(shù)據(jù)庫(kù)共同注釋到的序列有25010條。

圖2 四大數(shù)據(jù)庫(kù)注釋穿心蓮全轉(zhuǎn)錄組序列韋恩圖

Nr注釋結(jié)果顯示，與穿心蓮最相近的物種中的前三位是唇形目的芝麻(Sesamum indicum)（比對(duì)到34684條序列）、管狀花目的猴面花(Erythranthe guttata（)比對(duì)到7101條序列）以及管狀花目的旋蒴苣苔(Dorcoceras hygrometricum（)比對(duì)到1790條序列）。

KOG注釋結(jié)果發(fā)現(xiàn)，穿心蓮蛋白除參與維持正常生命活動(dòng)所必須的途徑（RNA進(jìn)程、轉(zhuǎn)錄進(jìn)程和常規(guī)功能）外，主要富集的功能途徑為：翻譯后修飾及蛋白質(zhì)轉(zhuǎn)換、信號(hào)轉(zhuǎn)導(dǎo)及胞內(nèi)運(yùn)輸、分泌和囊泡運(yùn)輸。眾所周知，植物次生代謝產(chǎn)物是由公共前體經(jīng)特異性蛋白修飾而成，如抗生素、激素和色素等。這些產(chǎn)物形成后可貯存在液泡和細(xì)胞壁中，部分也可運(yùn)輸?shù)襟w外環(huán)境。在貯存和運(yùn)輸?shù)倪^程中，需要細(xì)胞內(nèi)的信號(hào)轉(zhuǎn)導(dǎo)和囊泡運(yùn)輸?shù)裙δ艿膸椭Ｒ虼?，該結(jié)果體現(xiàn)了穿心蓮細(xì)胞生命活動(dòng)主要為次生代謝的特點(diǎn)。KEGG聚類進(jìn)一步證明KOG數(shù)據(jù)庫(kù)的注釋結(jié)果。如表1。穿心蓮功能基因聚類最多的途徑為代謝途徑、次生代謝產(chǎn)物合成途徑和抗生素合成途徑。

表1 穿心蓮全轉(zhuǎn)錄組信息的KEGG注釋

2.3 穿心蓮全長(zhǎng)轉(zhuǎn)錄組的基因結(jié)構(gòu)分析

經(jīng)過CDS和蛋白預(yù)測(cè)，穿心蓮全長(zhǎng)轉(zhuǎn)錄組信息得到進(jìn)一步注釋。本研究共識(shí)別出64401條mRNA序列和12411條長(zhǎng)鏈非編碼RNA(lncRNA)序列。其中，4008條序列的預(yù)測(cè)翻譯產(chǎn)物含信號(hào)肽結(jié)構(gòu)。

AS事件可在轉(zhuǎn)錄后水平豐富基因的轉(zhuǎn)錄本種類。利用生物信息學(xué)軟件將全長(zhǎng)轉(zhuǎn)錄組數(shù)據(jù)中的序列進(jìn)行聚類，可比較同一聚類下的任意兩條轉(zhuǎn)錄本發(fā)生AS事件的位置。經(jīng)統(tǒng)計(jì)，共有11273個(gè)穿心蓮基因發(fā)生了AS事件，體現(xiàn)出其轉(zhuǎn)錄本的多樣性。其中，產(chǎn)生2個(gè)AS亞型的基因最多，超過總數(shù)的1/4（2845個(gè)，25.24%）。數(shù)量次之的為產(chǎn)生3個(gè)和4個(gè)AS亞型的基因，分別有1935個(gè)(17.16%)和1389個(gè)(12.32%)。本研究中出現(xiàn)次數(shù)最多的AS事件為內(nèi)含子保留事件，其中啟動(dòng)子可變型內(nèi)含子保留事件發(fā)生頻率最高。其次為5’端可變事件。

穿心蓮內(nèi)酯的二萜類前體骨架(E,E,E)-香葉基香葉基二磷酸酯[(E,E,E)-geranylgeranyl pyrophosphate，GGPP]的合成原料為質(zhì)體內(nèi)的MEP途徑和胞質(zhì)中的MVA途徑產(chǎn)生的異戊烯基二磷酸(IPP)和二甲基烯丙基二磷酸(DMAPP)[17]。本研究對(duì)MEP和MVA途徑基因進(jìn)行AS事件分析。由圖3可知，兩途徑的基因均發(fā)生了AS事件。GGPP合成前期的基因主要以內(nèi)含子保留和5’端可變?yōu)橹鳌Ｉ蒊PP和DMAPP后，香葉基香葉基焦磷酸合成酶(GGPPS)基因和法呢基二磷酸酯合成酶(FPPS)基因發(fā)生了內(nèi)含子保留、3’端可變和外顯子跳躍AS事件。由表2可知，在發(fā)生內(nèi)含子保留的 GGPPS 基因中，Isoform0050046、Isoform0042454、Isoform0041790、Isoform0053892、Isoform0054050 和Isoform0054809為啟動(dòng)子可變式內(nèi)含子保留AS事件產(chǎn)生的轉(zhuǎn)錄本亞型。

表2 GGPP合成相關(guān)基因AS事件統(tǒng)計(jì)表

圖3 穿心蓮MEP和MVA途徑的AS事件示意圖

有2803條序列的預(yù)測(cè)蛋白為TF家族成員，含量最多的前10類如圖4所示。有研究表明，bHLH、basic leucine zipper(bZIP)、MYB及WRKY家族成員可直接調(diào)節(jié)次生代謝途徑基因的表達(dá)水平[18]。

圖4 含量排名前十位的TF家族

SSR位點(diǎn)廣泛存在于基因組中，具有重復(fù)性好、共顯性高、穩(wěn)定性強(qiáng)且多態(tài)性豐富等特點(diǎn)，可作為分子標(biāo)記來(lái)輔助育種[19]。本研究在穿心蓮全長(zhǎng)轉(zhuǎn)錄組中共找到45550處SSR位點(diǎn)，分布在27236條序列上。包含2個(gè)或以上SSR位點(diǎn)的序列有10569條。含2 nts重復(fù)單元的位點(diǎn)有20235個(gè)；含3 nts重復(fù)單元的位點(diǎn)有14403個(gè)；含4 nts重復(fù)單元的位點(diǎn)有6638個(gè)；含5 nts重復(fù)單元的位點(diǎn)有2411個(gè)；含6 nts重復(fù)單元的位點(diǎn)有1863個(gè)。單元重復(fù)次數(shù)主要集中在4～ 7次。包含不同堿基數(shù)重復(fù)單元的復(fù)合式SSR位點(diǎn)有7076個(gè)。由圖5可知，包含2 nts的AT/AT是穿心蓮轉(zhuǎn)錄組SSR中占比最高(27.7%)的重復(fù)單元，其次為AG/CT重復(fù)單元(12.8%)。

圖5 SSR位點(diǎn)的數(shù)量占比

3 討論與結(jié)論

穿心蓮內(nèi)酯作為穿心蓮的主要二萜類次級(jí)代謝物質(zhì)，具抗炎、抗感染及抗過敏等多重功效[20]。從轉(zhuǎn)錄水平分析穿心蓮次生代謝途徑的分子機(jī)制可幫助解析其藥用價(jià)值。基于Pacbio平臺(tái)的SMRT轉(zhuǎn)錄組測(cè)序無(wú)需組裝即可捕獲單個(gè)RNA分子的全長(zhǎng)序列。相較于SGS技術(shù)，SMRT測(cè)序能夠更準(zhǔn)確地反應(yīng)被測(cè)序物種轉(zhuǎn)錄組的實(shí)時(shí)信息，如發(fā)現(xiàn)更多的AS事件和剪切位點(diǎn)、發(fā)現(xiàn)新基因以及豐富基因組注釋等[10]。本研究利用三代SMRT測(cè)序技術(shù)對(duì)穿心蓮不同組織進(jìn)行全長(zhǎng)轉(zhuǎn)錄組測(cè)序分析，共得到高質(zhì)量序列76812條。本次測(cè)序N50值為1414 bp，大于1000 bp，說(shuō)明測(cè)序結(jié)果組裝完整性較高[11]。功能注釋結(jié)果顯示，穿心蓮蛋白主要參與翻譯后修飾、信號(hào)轉(zhuǎn)導(dǎo)及運(yùn)輸過程，體現(xiàn)次生代謝產(chǎn)物的前體在形成后的修飾和貯存過程。

TF是基因表達(dá)的開關(guān)，可通過與啟動(dòng)子結(jié)合來(lái)加速或緩和次生代謝進(jìn)程[21]。在中藥鼠尾草(Salvia Sclarea)中，SsWRKY18、SsWRKY40及SsMYC（bHLH家族成員）基因促進(jìn)二萜類物質(zhì)合成，過表達(dá)株系具細(xì)菌和真菌抗性[22]。番茄(Solanum lycopersicum)bZIP型TF可調(diào)控花青素積累[23]。中藥黃芩(Scutellaria Baicalensis)的SbMYB8基因可促進(jìn)黃酮類物質(zhì)合成[24]。同時(shí)，腺毛體具有合成和分泌次生代謝產(chǎn)物的能力。煙草(Nicotiana Benthamiana)Cys2-His2(C2H2)鋅指TF可正調(diào)控腺毛體的發(fā)育起始進(jìn)程[25]。本研究對(duì)穿心蓮全轉(zhuǎn)錄本信息進(jìn)行TF家族聚類。結(jié)果表明，以上參與次生代謝途徑的重要TF家族成員基因在穿心蓮轉(zhuǎn)錄組中的含量占比位居前10。該結(jié)果為后續(xù)穿心蓮次生代謝調(diào)控途徑的TF挖掘提供了參考依據(jù)。

含內(nèi)含子的前體mRNA可通過高度受控的AS過程從單個(gè)基因上產(chǎn)生多種轉(zhuǎn)錄本，從而增加胞內(nèi)或胞間轉(zhuǎn)錄組及蛋白組的復(fù)雜性[26]。AS事件主要分為：外顯子跳躍型、3’端可變型、5’端可變型及內(nèi)含子保留型四大類。其中包括外顯子互斥型（不同亞型跳躍保留不同的外顯子）、可變啟動(dòng)子型（5’端內(nèi)含子保留導(dǎo)致首個(gè)外顯子不同的多種亞型）及可變終止子型（3’端內(nèi)含子保留導(dǎo)致最后一個(gè)外顯子不同的多種亞型）[27]。內(nèi)含子保留型是植物中較為常見的AS事件，擬南芥(Arabidopsis Thaliana)和水稻(Oryza sativa)的內(nèi)含子保留事件發(fā)生概率高達(dá)64%和55%[28]。Xu等[29]利用三代SMRT測(cè)序技術(shù)從丹參的根中鑒定出參與迷迭香酸和丹參酮生物合成的酶編碼基因中的AS事件。本研究顯示，穿心蓮內(nèi)酯二萜類前體GGPP的合成途徑基因發(fā)生的主要AS類型為內(nèi)含子保留，與前人研究一致。其中，1個(gè)GGPPS基因出現(xiàn)了可變啟動(dòng)子式內(nèi)含子保留AS事件，產(chǎn)生6個(gè)首位外顯子不同的AS亞型。該結(jié)果體現(xiàn)穿心蓮全長(zhǎng)轉(zhuǎn)錄組的蛋白編碼多樣性。與本研究類似，類胰島素生長(zhǎng)因子1基因的轉(zhuǎn)錄由兩個(gè)啟動(dòng)子控制。啟動(dòng)子的交替使用配合其他位置的AS事件，使該基因產(chǎn)生若干第一外顯子不同的mRNA亞型[27]。

SSR分子標(biāo)記技術(shù)的核心序列為1～ 6個(gè)核苷酸的串聯(lián)重復(fù)序列，具物種特異性，可用于種質(zhì)資源鑒定[30]。本研究在穿心蓮全長(zhǎng)轉(zhuǎn)錄組中共找到45550處SSR位點(diǎn)；分布于27236條序列上；包含2個(gè)或以上SSR位點(diǎn)的序列有10569條；單元重復(fù)次數(shù)主要為4～ 7次；包含2 nts的AT/AT和AG/CT為優(yōu)勢(shì)重復(fù)單元。該結(jié)果有助于穿心蓮的SSR分子標(biāo)記開發(fā)。

本研究利用SMRT超長(zhǎng)單分子測(cè)序技術(shù)來(lái)解析穿心蓮幼苗的全長(zhǎng)轉(zhuǎn)錄組信息，發(fā)現(xiàn)其體內(nèi)功能基因主要參與次生代謝物質(zhì)的合成與貯存相關(guān)的生物途徑?；钴S的次生代謝使其主要藥效成分穿心蓮內(nèi)酯及其衍生物得到充分的積累。功能基因啟動(dòng)子通過結(jié)合直接促進(jìn)次生代謝的轉(zhuǎn)錄因子以及轉(zhuǎn)錄后形成多種可變剪切亞型來(lái)維持上述次生代謝活性。此外，搜索到的穿心蓮SSR位點(diǎn)可幫助學(xué)者對(duì)相關(guān)基因的遺傳特性進(jìn)行評(píng)估。綜上，本研究為穿心蓮內(nèi)酯的合成機(jī)制研究及遺傳多樣性評(píng)價(jià)提供了分子基礎(chǔ)。