陽春砂花絲、花柱轉錄組測序及生物信息學分析

楊博涵， 徐 杰，2， 湯麗云， 蘇 景， 李子翔， 何國振

[1.廣州中醫藥大學中藥學院，廣東廣州 510006； 2.廣東一方制藥有限公司/廣東省中藥配方顆粒企業重點實驗室，廣東佛山 528244；3. 華南農業大學生命科學學院，廣東廣州 510642； 4.陽春市農業試驗場(陽春市春砂仁試驗場)，廣東陽春 529600]

陽春砂(Amomumvillosum)是姜科(Zingiberaceae)豆蔻屬(Amomum)多年生常綠草本植物[1]，其干燥成熟果實為四大南藥之首——砂仁的主流品種，具有理氣安胎、溫脾止瀉、化濕開胃的功效[2]。作為傳統的大宗藥材，陽春砂每年的需求數量大于300萬kg，具有重要的藥用價值及經濟價值[3]。但陽春砂在生產上存在著嚴重的低產問題，其自然結實率僅有1.1%[4]，研究者對此進行了許多研究，認為較低的自然結實率與陽春砂小花的特殊花器結構有著密切的聯系。何國振等將陽春砂特殊的花器結構稱為假合蕊柱，即雌雄蕊貼合在一起，但并未完全合生，在處于開花的狀態下，唇瓣包裹著假合蕊柱，并且呈現出半抱合的狀態，與唇瓣間的距離極近，僅有1.3～2.0 mm，嚴重阻礙了相應蟲媒的授粉[3]。為提高產量，農戶在實際種植中主要采用人工授粉的方式進行勞作，而人工授粉成本高、勞作強度大、根狀莖及花序易被嚴重踩踏損傷等原因是傳統人工授粉方式固有的弊端，嚴重抑制了農戶的生產積極性。

細胞伸長或者細胞分裂所致的植物器官運動與內源激素及相關抑制劑作用相關，而細胞結構及生長速率的不對稱是花器官運動的基礎。有文獻研究報道，臘梅(Chimonanthuspraecox)花絲的運動來源于其兩側表面細胞生長速率的差異[5]。姜科植物馬來良姜(Alpiniamutica)花柱運動部位兩側細胞層數的差異是其花柱卷曲運動的基礎[6]，赤霉素(GA)、茉莉酸(JA)、吲哚-3-乙酸(IAA)等激素則有著調控植物雌雄蕊發育的作用[7-9]。何卓航等研究發現，隨著陽春砂小花的生長，其花絲、花柱遠、近軸側的細胞層數出現了差異，這種不對稱的結構是雌雄蕊運動的基礎，認為陽春砂雌雄蕊的相向運動致使其假合蕊柱的形成[10]；同時對花柱的近端以及遠端軸側的IAA水平進行測定，結果表明，生長時期在保持同樣的水平不同橫切部位中，IAA水平基本上呈現遠端軸側小于近端軸側的趨勢，故花柱運動的原因可能是來自于花柱近、遠軸側IAA水平的差異，該研究從生理層面初步闡述了陽春砂假合蕊柱的形成機制。然而，若要進一步探明陽春砂假合蕊柱形成機制，則須對陽春砂花器官運動及發育的分子機制作深入研究。

因此，本研究對陽春砂不同生長時期的花絲、花柱進行轉錄組測序，通過數據拼接、組裝的方式建立陽春砂花絲、花柱的轉錄組數據庫，將所得的Unigene進行功能注釋、分類以及簡單重復序列(SSR)分子標記。同時，重點關注并篩選可能影響陽春砂雌蕊、雄蕊運動的激素合成與信號轉導途徑關鍵基因。以期進一步解析陽春砂假合蕊柱產生的調控機制及分子機制，并為人為干預陽春砂的花器官結構，提高產量奠定理論基礎。

1 材料與方法

1.1 試驗材料

材料為大田栽培的不同生長時期陽春砂小花的花柱及花絲，陽春砂植株栽培于廣東省陽春市合水鎮那軟村陽春砂種植基地(22°17′N，112°01′E)，于2017年5—7月在該基地開展大田試驗。參照陳紅的方法[11]，按照長度的不同，將陽春砂小花的生長時期進行了劃分(表1)，并于陽春砂花期內，每天上午采摘小花，解剖并分離出不同時期小花的花絲與花柱，迅速置于液氮灌中，隨后轉移至超低溫冰箱(-80 ℃)進行保存。

表1 陽春砂小花生長時期劃分

1.2 試驗方法

1.2.1 RNA提取 分別提取陽春砂小花在各生長時期的花柱和與花絲的總RNA，重復3次。利用NanoDrop 2000和Agilent 2100 Bioanalyzer測定RNA的濃度、純度與完整性，數據檢測合格后的RNA用于構建轉錄組數據庫。此部分試驗與轉錄組的測序委托深圳華大基因科技服務有限公司完成。

1.2.2 cDNA文庫構建及轉錄組序列組裝 使用Oligo(dT)磁珠富集質檢合格且所有樣品的混合RNA。加入打斷試劑，以片段化的mRNA為模板合成1鏈、2鏈cDNA，配制反應體系，使接頭與cDNA連接。PCR反應及產物回收、擴增。PCR產物變性，充分混勻，得到單鏈環形產物，隨后PCR產物變性，即得到文庫。利用Agilent 2100 Bioanalyzer和ABI StepOnePlus real-time PCR System對所得文庫進行檢測，檢測合格后進行轉錄組denovo測序。所得reads通過Trinity軟件進行序列組裝，組裝序列的質量通過BUSCO軟件進行評估。

1.2.3 轉錄組基因功能注釋及數據挖掘 利用生物信息學方法分析獲得的陽春砂Unigene，為獲得全方位的基因功能信息，對組裝所得的Unigene在七大功能數據庫中進行注釋，包括NR、NT、 KOG/COG、GO、KEGG、SwissPro及Interpro。使用MIcroSAtellite (MISA)工具對Unigene進行SSR位點的挖掘。并在KEGG代謝通路中重點關注植物激素生物合成與信號轉導Unigene的注釋情況。

2 結果與分析

2.1 轉錄組數據組裝

使用Illumina HiSeq-2000平臺對不同時期的陽春砂花絲及花柱的cDNA進行測序。測序的結果顯示，從轉錄組中獲得了81.40 Mb的原始讀數，用過濾軟件SOAPnuke去除低質量的reads后共得到68.48 Mb Clean reads，最終獲得10.33 Gb的堿基總數。Q20和Q30的百分比分別為97.20%和93.16%(表2)。說明轉錄組的測序質量較高，可以滿足后續的生物信息學分析。

表2 測序數據質量分析

對轉錄本中Clean reads進行組裝，一共獲得138 590個轉錄本，包含了111 312 938個核苷酸的序列信息，這些片段長度的平均值為803 bp，N50為1 414 bp，GC含量為44.28%。對轉錄本進一步聚類并去除冗余的序列之后得到94 584條Unigene，總長度為 92 501 015 bp。Unigene的N50、N70和GC含量分別為 1 582bp、1 002 bp和44.35%，均大于其對應平均長度，說明組裝的結果良好(表3)。

表3 轉錄本和單基因簇統計分析

2.2 陽春砂轉錄組基因總體注釋情況

將所獲得的結果在七大功能數據庫進行注釋，結果見表4。結果顯示，陽春砂共注釋到94 584條Unigene。其中，共有58 669條Unigene被NR數據庫注釋，占比最多，達到總Unigene的62.03%；NT數據庫有40 205條，占42.51%；SwissProt有40 318條，占總Unigene的42.63%；KEGG數據庫與KOG數據庫分別有44 000條和45 892條Unigene，各占46.52%與48.52%；Interpro有45 793條，占48.42%。GO數據庫注釋到的基因最少，僅有 12 050 條，占總數的12.74%；所得比對結果顯示，在七大數據庫中均能成功注釋的Unigene共有 5 691 條，占總Unigene條數的6.02%。

表4 轉錄組基因注釋情況統計

2.2.1 NR數據庫功能注釋 NR數據庫的注釋結果(圖1)顯示，匹配最多的物種為小果野芭蕉(Musaacuminatasubsp.Malaccensis)，該物種注釋到的基因數量最多，共有45 257條，占比高達77.14%，證明陽春砂與該物種的同源性較高；其他物種依次為油棕(Elaeisguineensis)與海藻(Phoenixdactylifera)，分別有3 338、2 435條Unigene被注釋，占比分別為5.69%、4.15%，陽春砂與這2種植物的同源性相對較低。而剩余的13.02%則分布于其他物種中。

2.2.2 KOG 數據庫功能注釋 在KOG數據庫中，共有45 892條Unigene被注釋到(圖2)，共分為25個大類功能區，包括一般功能預測、信號轉導機制、轉錄等功能。在為數眾多且不同的功能分類中，注釋到的基因的數量差異較為顯著，一般功能預測類基因數量最多，共有11 177條Unigene，占比為24.36%，其次是信號轉導機制，有8 271條Unigene被注釋，占比為18.02%；除此之外，負責轉錄功能的Unigene有5 470個，占比為11.92%；有1 172條Unigene注釋到負責次生代謝產物生物合成、運輸和代謝功能區中，占比為2.55%。

2.2.3 GO數據庫功能注釋 GO數據庫中的功能分類注釋結果見圖3，結果表明，陽春砂花絲及花柱中共有12 050條Unigene注釋到不同的功能節點上，共涉及到生物學過程、細胞組分和分子功能3個大類，55個亞類。歸入到細胞組分的17個亞類中，以細胞、細胞部分、膜和細胞器功能的Unigene數量最多，分別有5 335、5 229、4 052、3 821條。在分子功能的14個亞類中，催化活性和結合注釋數量最多，分別有5 594、5 569條。而涉及生物學過程的24個亞類中，以代謝過程(6 077條)和細胞進程 (6 057 條)為主。結果表明，同一個Unigene可以注釋到多個功能結點上， 因此同一個功能分支的總注釋數大于注釋到該功能分支的總Unigene數。

2.2.4 KEGG數據庫功能注釋 對轉錄組測序獲得的Unigene進行KEGG代謝通路富集分析后發現，共有44 000條 Unigene得到注釋，涉及到137個通路。對富集顯著的前20條通路進行分析(表5)，其中，共有8 660條Unigene被代謝途徑通路所注釋，占比最多，達到注釋總數的19.68%。而次生代謝產物的生物合成通路共有4 191條Unigene被注釋到，占注釋總數的9.53%，另外，本研究主要關注的植物激素轉導通路一共有1 656條Unigene被注釋到，占比為3.76%。

在所有的代謝通路中，共有23個與次生代謝相關的代謝通路，對其所有的Pathway進行分析后(表6)發現，一共有4 191個Unigene被次生代謝產物的生物合成途徑所注釋；而在本研究主要關注的植物激素合成途徑中，62個Unigene被油菜素內酯生物合成途徑所注釋，另有59個Unigene被玉米素生物合成途徑所注釋。

2.2.5 SSR特征分析 轉錄組測序所得Unigene的SSR的檢測最終結果顯示，共有18 895個SSR被檢測到(圖4)。總共有14種多堿基重復SSR，其中， 三核苷酸重復的SSR 數目最多，共有6 313個，占SSR位點總數的33.41%；三核苷酸重復中AGG/CCT重復基元最多，共有1 434個；單核苷酸和二核苷酸重復分別有6 238個(33.01%)和5 128個(27.14%)，二核苷酸中重復基元最多的是AG/CT，共有2 909個；六核苷酸重復數目共有546個(2.89%)，四核苷酸重復330個(1.75%)和五核苷酸重復340個(1.80%)。

表5 次生代謝相關代謝通路的注釋情況

表6 次生代謝相關代謝通路的注釋情況

2.3 重點關注通路的Unigene 注釋情況

2.3.1 植物激素信號轉導途徑中Unigene注釋情況 在KEGG代謝通路中，共有1 760條Unigene注釋到植物激素的信號轉導途徑中(表7)。赤霉素的信號轉導途徑中注釋到的5個基因，對應的Unigene最多，共有438個；其中，植物光敏色素互作因子(PIF)對應的Unigene最多，達到318個；而水楊酸信號轉導途徑注釋到3個基因，包括非表達病程相關蛋白(NPR1)、TGA轉錄因子及病程相關蛋白1(PR1)基因，但對應的Unigene最少，僅有74個；茉莉酸的信號轉導途徑中注釋到的茉莉酸ZIM結構域蛋白(JAZ)、茉莉酸氨基酸合成酶(JAR1)、MYC2轉錄因子和茉莉酸受體蛋白(COI1)一共對應了176條Unigene；乙烯和油菜素內酯信號轉導途徑注釋到的基因最多，均為8個，分別注釋到139條及191條Unigene；其他植物激素中，生長素、細胞分裂素以及脫落酸的信號轉導途徑分別注釋到329、245、168條Unigene；統計結果中，各激素信號轉導途徑關鍵基因均有注釋，如細胞分裂素中的磷酸轉移蛋白(AHP)、生長素中的生長素載體(AUX1)、生長素/吲哚乙酸蛋白(Aux/IAA)轉錄家族以及脫落酸信號轉導途徑中的ABF轉錄因子等。綜上，研究較為深入的植物激素基本均在KEGG通路中有所注釋，且其信號轉導途徑的關鍵基因都有一定數量的Unigene被注釋到。

表7 轉錄組中植物激素信號轉導途徑Unigene注釋情況

2.3.2 植物激素合成途徑中Unigene注釋情況 轉錄組中激素合成途徑酶基因Unigene注釋情況(表8)顯示，赤霉素注釋的酶基因最多，包括內根-貝殼杉烯合成酶(KS)、古巴焦磷酸合成酶(CPS)、內根-貝殼杉烯氧化酶(KO)及GA3氧化酶(GA3ox)等8個合成通路中的關鍵酶基因，一共對應67條Unigene；油菜素內酯與水楊酸合成途徑中均有1個酶基因被注釋到，為異分支酸合酶(ICS)與油菜素內酯-6-氧化酶2(CYP85A2)，分別對應9條和24條Unigene；茉莉酸合成途徑2個關鍵酶基因丙二烯氧化物合成酶(AOS)以及丙二烯氧化物環化酶(AOC)分別注釋到6條和7條Unigene，生長素一共注釋到42條Unigene，色氨酸轉氨酶(TAA1)及吲哚-3-丙酮酸單加氧酶(YUC)分別注釋到了其中22條和20條；乙烯合成途徑中的乙烯合成前體合成酶(ACS)與乙烯合成前體氧化酶(ACO)分別對應9和14條Unigene；脫落酸及玉米素的生物合成途徑中均注釋到3個關鍵酶基因，分別有29條和45條Unigene。

表8 轉錄組中植物激素生物合成途徑Unigene 注釋情況

3 討論與結論

近年來，轉錄組測序技術發展迅速，同時伴隨著多個生物信息學分析平臺的加入，大量研究即便在沒有基因組數據支撐的前提下，通過轉錄組測序所得的結果也可在挖掘植物未知基因、明確相關生理功能的代謝途徑及基因調控機制等方面提供海量相關信息。無需參考基因組的數據便可對目標物種進行分子生物學方面的研究，已成為中草藥材基因信息挖掘的重要研究手段[12]。利用該技術，前人已從多種植物中挖掘出有關花器官發育的基因，并分析了相關基因的調控機制。李梅等利用雌蕊缺失茶樹(Camelliasinensis)花的花、花蕾、花芽作為材料，基于轉錄組測序技術，共發掘出了34個花器官發育差異表達的相關基因，包括KNOX家族基因、WUS類基因、花器官發育 ABCDE 模型的相關基因等，證明茶樹花即使是在雌蕊缺失的情況下，有關雌蕊發育的基因也是同樣存在[13]。而在長瓣兜蘭(Paphiopedilumdianthum)花蕾和花朵的轉錄組測序中發現，相較于花蕾期，花朵時期中調控花朵發育相關的基因，如AG、C2H2-ZEP等基因的表達量明顯下調，大量基因表達下調的原因可能是在花朵時期，花的各部分器官已完成分化，調控器官分化基因的相關任務已完成[14]所致。位明明等通過該技術初步明確了多個參與編碼橡膠樹(Heveabrasiliensis)花器官發育的基因家族，如Mads-box、MYB及AP2等基因家族[15]。

本研究通過RNA-Seq技術構建了陽春砂不同時期花絲以及花柱總樣品的轉錄組數據庫，共獲得94 584條Unigene，Unigene的平均長度為842 bp，N50為1 582 bp，GC含量為44.35%。其中，共有 62 174 條Unigene被注釋到GO、KEGG、NR和KOG等七大數據庫，占總Unigene的 65.73%。注釋率較高，測序所得結果與拼接質量較好。這些注釋的Unigene為陽春砂的代謝途徑、基因功能分類、花器官發育及雌雄蕊運動分析等方面提供參考依據。此外，一共有32 410條Unigene未被注釋到，占總Unigene的34.27%，推測這些Unigene可能為陽春砂中的非編碼 RNA 序列或現有基因數據庫尚未完善所致。在NR數據庫中共注釋到58 669條Unigene，其中77.14%的Unigene被注釋到姜科植物小果野芭蕉中，兩者比對所得同源基因數目較多的原因可能是親緣關系較近。轉錄組所得Unigene在GO功能分類中，注釋到生物學過程、細胞組分以及分子功能三大類，分別有 24、17、14個亞類，主要集中在代謝過程、催化活性、細胞部分和細胞過程等功能。而在KOG數據庫的分析結果中，共有 45 892 條Unigeine得到了注釋，在25個功能大類中，一般功能預測的Unigene最多，共有11 117條，信號轉導機制、轉錄以及翻譯后的修飾、蛋白質轉換與代謝功能次之，核結構與細胞運動最少，分別僅有346條和75條Unigene。通過KEGG數據庫分析發現，被注釋到44 000條Unigene參與了23類共計137個KEGG代謝通路，這些通路主要集中在代謝途徑、次生代謝產物的生物合成、RNA轉運以及植物激素信號轉導等，可以通過所注釋的通路全面了解陽春砂花絲及花柱的代謝途徑信息。同時，在陽春砂花絲、花柱轉錄組中一共挖掘到18 895個SSR位點，在6種不同的核苷酸重復類型中均有分布，證明本研究轉錄組中位點的類型豐富，分布密度較大，具有良好的多態性潛能；在二核苷酸和三核苷酸中，重復最多的基序是分別是AG/CT和AGG/CTT，該結果與王煥的研究結果[16]一致。

植物激素是植物通過自身代謝所產生的一些有機信號分子，可在低濃度下產生明顯的生理效應，并在合成部位發揮功能[17-18]。而花器官的運動受到激素調節，Luo等研究發現，外源施加IAA后，會明顯影響姜科植物花柱的彎曲程度[19]。有些植物激素也會調控姜科植物藍豬耳二長雄蕊的翻轉與伸長[20]。此外，其他植物激素同樣也在雌雄蕊發育中起著重要的調控作用。JA可調控花藥的開裂，JA生物合成通路的相關基因DAD1與OPR3在生長素感知缺陷的擬南芥突變體中的表達量提高，可提前使得花藥產生開裂，說明此過程是一種負向調節，即生長素可以通過介導JA的生物合成來調節花藥的發育[21-22]。在模式植物擬南芥的其他研究中，茉莉酸與赤霉素信號轉導相關基因也在擬南芥雄蕊的發育過程中起到了重要的調控作用[8，23-24]。在本研究重點關注的植物激素合成與信號轉導途徑中一共有2012條Unigene被注釋，囊括了赤霉素、生長素、茉莉酸及脫落酸等9種重要的植物激素，為進一步從分子層面闡述陽春砂花器官的發育及運動提供了大量候選基因。

筆者所在課題組于2012年提出通過生理手段改變陽春砂的花器結構，以適應更多種類的昆蟲傳粉而提高產量的研究思路[25]。近10年來，已經在陽春砂的花芽分化規律[26]、花器結構特征[3]、假合蕊柱的形成過程[10]、某些生殖生物學特性[27]、落果規律[28]等多方面進行了研究，明晰了陽春砂生產過程中由于特殊花器官結構導致傳粉昆蟲種類受限和低產的原因。因此，下一步的研究方向需對本研究中所獲有關陽春砂花器官的發育及運動的關鍵候選基因進行更深層次的功能分析與驗證，并對陽春砂假合蕊柱形成與相關基因表達模式的關系進行分析與探討，不僅可為深入開展陽春砂花器官發育與運動相關基因的克隆及相互作用的模式分析提供數據基礎，也為進一步深入揭示陽春砂假合蕊柱的分子調控機制提供了理論依據，以期形成人為干預的效果進而改變陽春砂花器官的形態，拓寬陽春砂唇瓣與假合蕊柱之間的距離，為陽春砂的授粉蟲媒打造一個相對有利的授粉環境，從而提高自然結實率，達到增加農民收入的目的。