穿心蓮全長轉錄組測序及特性分析

孫銘陽，徐世強，顧 艷，梅 瑜，周 芳，李靜宇，王繼華

（廣東省農業科學院作物研究所/廣東省農作物遺傳改良重點實驗室，廣州510640）

0 引言

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

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

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

1 材料與方法

1.1 實驗材料

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

1.2 全長轉錄組文庫構建

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

1.3 全長轉錄組測序數據篩選

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

1.4 全長轉錄組功能注釋

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

1.5 基因結構分析

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

2 結果與分析

2.1 穿心蓮全長轉錄組測序的數據篩選

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

圖1 FLNC序列長度分布圖

2.2 穿心蓮全長轉錄組的功能注釋

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

圖2 四大數據庫注釋穿心蓮全轉錄組序列韋恩圖

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

KOG注釋結果發現，穿心蓮蛋白除參與維持正常生命活動所必須的途徑（RNA進程、轉錄進程和常規功能）外，主要富集的功能途徑為：翻譯后修飾及蛋白質轉換、信號轉導及胞內運輸、分泌和囊泡運輸。眾所周知，植物次生代謝產物是由公共前體經特異性蛋白修飾而成，如抗生素、激素和色素等。這些產物形成后可貯存在液泡和細胞壁中，部分也可運輸到體外環境。在貯存和運輸的過程中，需要細胞內的信號轉導和囊泡運輸等功能的幫助。因此，該結果體現了穿心蓮細胞生命活動主要為次生代謝的特點。KEGG聚類進一步證明KOG數據庫的注釋結果。如表1。穿心蓮功能基因聚類最多的途徑為代謝途徑、次生代謝產物合成途徑和抗生素合成途徑。

表1 穿心蓮全轉錄組信息的KEGG注釋

2.3 穿心蓮全長轉錄組的基因結構分析

經過CDS和蛋白預測，穿心蓮全長轉錄組信息得到進一步注釋。本研究共識別出64401條mRNA序列和12411條長鏈非編碼RNA(lncRNA)序列。其中，4008條序列的預測翻譯產物含信號肽結構。

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

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

表2 GGPP合成相關基因AS事件統計表

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

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

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

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

圖5 SSR位點的數量占比

3 討論與結論

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

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

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

SSR分子標記技術的核心序列為1～ 6個核苷酸的串聯重復序列，具物種特異性，可用于種質資源鑒定[30]。本研究在穿心蓮全長轉錄組中共找到45550處SSR位點；分布于27236條序列上；包含2個或以上SSR位點的序列有10569條；單元重復次數主要為4～ 7次；包含2 nts的AT/AT和AG/CT為優勢重復單元。該結果有助于穿心蓮的SSR分子標記開發。

本研究利用SMRT超長單分子測序技術來解析穿心蓮幼苗的全長轉錄組信息，發現其體內功能基因主要參與次生代謝物質的合成與貯存相關的生物途徑。活躍的次生代謝使其主要藥效成分穿心蓮內酯及其衍生物得到充分的積累。功能基因啟動子通過結合直接促進次生代謝的轉錄因子以及轉錄后形成多種可變剪切亞型來維持上述次生代謝活性。此外，搜索到的穿心蓮SSR位點可幫助學者對相關基因的遺傳特性進行評估。綜上，本研究為穿心蓮內酯的合成機制研究及遺傳多樣性評價提供了分子基礎。