薏苡轉錄組測序及基因功能注釋

歐陽雨晴，李玲玲，石好宇，趙藝，朱甫臻，尹艷燕，王靈芝（北京中醫藥大學生命科學學院，北京100029）

轉錄組（transcriptome）是指細胞在特定發育階段或生理條件下，所表達的全部RNA 的總和，可反映相同基因在不同環境下表達的差異，揭示不同基因間的相互作用及其各自的功能[1-2]。轉錄組測序，又稱為RNA 測序（RNA-seq），將樣本總RNA 反轉錄成cDNA 后，進行高通量測序來確定樣本中整體轉錄組的表達情況。自2005年以來，以Solexa 技術、SOLiD 技術和454 技術為代表的第二代測序技術，推動了RNA-seq 技術的迅速發展。RNA-seq 可在限定時間內全面而快速地獲取特定樣品的幾乎所有mRNA 序列信息，該技術在基因功能分析、特異表達基因挖掘、信號途徑構建、代謝途徑鑒定、基因家族確定等方面發揮重要作用[3]。RNA-seq 技術已經成功應用于許多藥用植物轉錄組信息的收集，如黃連[4]、人參[5]、黃三七[6]、銀杏[7]、紅豆杉[8]和阿育魏實[9]等，在植物有效成分生物合成的功能基因挖掘和分子標記育種方面取得了較大進展[10]。

薏苡（Coix lacryma-jobiL.var.mayuenStapf），是禾本科（Gramineae）玉蜀黍族（Maydeae），薏苡屬CoixL.植物[11]，為藥食同源作物，具有利水滲濕、健脾止瀉之功能，中醫臨床常用于治療水腫、排尿困難和腹瀉等[12]。薏苡仁含有多種化學成分如多糖、蛋白質、脂質、多酚、類胡蘿卜素、植物甾醇、螺甾內酯和內酰胺[13]，具有多種藥理功能，如抗腫瘤[14]、免疫調節[15]、抗氧化[16]、降血壓[17]等功能。其中，薏苡仁油因具有顯著的抗腫瘤活性，已被研制成康萊特注射液，廣泛應用于臨床，并被美國食品和藥物管理局（FDA）批準在美國進行臨床試驗[18]。薏苡仁油的主要活性成分為三酰甘油類，占薏苡仁油含量的86.96%～94.39%[19]，其中甘油三油酸酯為《中國藥典》（2020年版）薏苡仁質量標準的指標性成分。研究表明1，3-二油酸-2-亞油酸甘油酯[20]和1-棕櫚酸-2-亞油酸-3-油酸甘油酯[21]具有顯著抗腫瘤活性。根據溶解度不同，種子儲藏蛋白分為清蛋白、球蛋白、醇溶蛋白和谷蛋白，谷蛋白約占總蛋白含量的40%[22]，具有多種潛在的生物活性。美國山核桃谷蛋白源多肽有降壓，抗氧化和抗腫瘤活性[23]。莧菜谷蛋白水解產物有降血糖作用[24]。胡桃楸谷蛋白水解產物具有免疫調節作用[25]。本課題組前期研究表明，薏苡仁谷蛋白源小分子肽具有免疫調節活性[26]。

薏苡仁的生物活性主要取決于其品種。薏苡在中國除青海和寧夏以外的大多數省份中均有分布，藥材的質量與生產區域的環境密切相關[27]，由于地理環境和栽培條件的差異，使得我國薏苡仁種質資源極為豐富多樣。但薏苡基因序列信息缺乏，阻礙其遺傳育種與藥用價值的開發。本研究采用Illumina Hiseq 測序平臺技術，構建薏苡轉錄組數據庫，并進行功能注釋，為谷蛋白和三酰甘油等成分合成關鍵基因提供序列信息。

1 材料與方法

1.1 藥材

薏苡種植于北京中醫藥大學的藥材種植園，采集不同發育時期的根、葉、穗、幼苗及籽粒，液氮冷凍后，－80℃保存，用于后續RNA 提取。

1.2 儀器與試藥

EASY spin 植物RNA 快速提取試劑盒（北京博邁德生物技術有限公司）；1%瓊脂糖凝膠電泳檢驗RNA 完整性；Qubit Fluorometer 熒光定量儀（美國BioTek）；NanoDrop 超微量分光光度計（美國Thermo）；Agilent 2100 生物分析儀（美國安捷倫科技有限公司）。

1.3 轉錄組文庫的建立和測序分析

樣品檢測合格后，等量混合EASY spin 植物RNA 快速提取試劑盒提取的各薏苡組織的總RNA，經帶有Oligo（dT）的磁珠純化mRNA 以構建cDNA 文庫。 采用Illumina HiSeq 4000 以PE 100 的測序讀長對cDNA 進行測序，得到原始數據（Raw reads）。使用SOAPnuke 處理Raw reads，去除低質量，包含接頭以及未知堿基N含量大于5%的reads，得到Clean reads。使用Trinity 對Clean reads 進行De novo組裝，使用Tgicl 對組裝的轉錄本進行聚類，去除冗余獲得unigene。N50 值用于評估組裝過程中的連續性。

1.4 轉錄本功能注釋

將組裝得到的薏苡轉錄本數據，通過BLAST 軟件分別與NT、NR、COG、KEGG 以及SwissProt 數據庫進行比對；基于NR 注釋，使用Blast2GO 軟件獲得GO 注釋，使用InterProScan5軟件獲得InterPro 注釋。

1.5 生物信息學分析

使用Blast 軟件將unigene 序列按照 NR、SwissProt、KEGG、COG 數據庫的優先順序進行比對，得到unigene 編碼區5'-3'方向的核酸序列，即CDS。使用ESTScan 軟件對未比對上的unigene 序列，進行CDS 預測。使用MISA 軟件對unigene 進行SSR 分析，并用Primer 3.0 軟件設計SSR 引物。使用HISAT 軟件將Clean reads 比對到unigene，然后通過變異檢測軟件GATK 識別SNP 位點。

1.6 基因表達量計算

基因表達量用 FPKM（Fragments Per Kilobase of transcript per Million mapped reads）表示。使用Bowtie2 軟件將Clean reads 比對到unigene，然后利用RSEM 計算各樣品的基因表達水平。FPKM計算公式如下：

FPKM＝cDNA Fragments/[Mapped Fragments（Millions）×Transcript Length（kd）]

式中，cDNA Fragments 表示比對到某一轉錄本上片段數目，即雙端Reads 數目；Mapped Fragments（Millions）表示比對到轉錄本上的片段總數，以106為單位；Transcript Length（kd）表示轉錄本的長度，以1×103個堿基為單位。

2 結果

2.1 測序結果與序列組裝

電泳圖譜（見圖1）表明樣品RNA 完整性較好。根據HiSeq 測序樣品質量標準，各樣品均為A級，滿足后續cDNA 文庫構建要求。使用Illumina HiSeq 4000 平臺共得到了65.60 Mb Raw reads，過濾后得到99.95%的Clean reads（65.57 Mb），Q20為96.09%。讀本中堿基分布均勻沒有明顯的AT或GC 分離現象（見圖2A），表明該RNA-seq 質量較高。使用Trinity 軟件對Clean reads 進行De novo組裝，共獲得102 092 個transcripts。經Tgicl軟件將transcrips 聚類并除去冗余之后，獲得了76 164 個unigene，平均長度為1023 bp，GC 含量為49.36%，N50 為1700 bp，表明組裝效果較好（見表1），其中序列長度為300 bp 數目最多，有18 321 條（見圖2B）。

表1 組裝過程質量指標Tab 1 Quality index of assembly process

圖1 薏苡各組織RNA 電泳圖Fig 1 RNA electrophoretogram of Coix tissues

圖2 Clean reads 堿基含量（A）和unigene 長度（B）分布圖Fig 2 Clean reads base content（A）and unigene length（B）distribution

2.2 基因功能注釋

將unigene 序列與NT、NR、COG、KEGG、SwissProt、GO 和InterPro 七大功能數據庫進行比對，共有59 460 個unigene 獲得注釋，占總數的78.07%。其中，在NT 數據庫中，獲得注釋的unigene 數量最多，有56 121 個，占73.68%。而在COG 數據庫中，獲得注釋的unigene 數量最少，僅有20 670 個，占27.14%（見表2）。unigene 功能注釋Venn 圖（見圖3A）顯示，共有14 177 個unigene在NR、Interpro、Swissprot、COG 和KEGG 5 個數據庫同時得到注釋，占unigene 總數的23.84%。

表2 Unigenes 功能注釋Tab 2 Function annotation of unigenes

2.2.1 Unigene 的NR 數據庫比對分析 將unigene與NR 數據庫進行相似序列匹配，結果如圖3B 所示，共有49 092 個unigene 在NR 數據庫中注釋成功。其中薏苡與高粱（Sorghum bicolor）同源的序列最多，占NR 總注釋的51.73%，其次匹配的物種有玉米（Zea mays，26.29%）、小米（Setaria italica，10.16%）和水稻（Oryza sativaJaponica Group，3.21%）。

2.2.2 Unigene 的GO 功能分類分析 如圖3C 所示，共有27 407 個unigene 在GO 數據庫中注釋成功，得到156 122 條對應關系。三大功能中生物學過程最多，為67 477 條（43.22%）；細胞組分數目為55 713 條（35.69%）；分子功能最少為32 932條（21.09%）。這三大功能可詳細劃分為56個分支，生物學過程22個、細胞組分17個、分子功能17個。在生物學過程類型中，與新陳代謝過程相關的unigene 為15 060 條，所占比例最高；在細胞組分功能分類中，細胞與細胞成分所涉及的unigene 數量一致且最多，有14 043 條；在分子功能類型中，注釋上的unigene 涉及結合、抗氧化、催化、翻譯調節及營養儲存等多方面功能活性，其中，涉及到結合功能的unigene 有14 728 條，所占比例最高，對應蛋白質標簽功能的unigene 最少（1 條）。

2.2.3 Unigene 的COG 注釋結果分析 如圖3D所示，共有20 670 個unigene 在COG 數據庫中注釋成功，根據其功能的不同可分為25 類。一般功能基因類別中unigene 數量最多，占unigene總量的27.33%；核結構類unigene 有4 條，所占比例最小。此外，翻譯、核糖體結構與生物合成（19.86%），轉錄（17.50%），信號轉導機制（12.13%），復制、重組與修復（16.72%），翻譯后修飾、蛋白質轉運及分子伴侶（12.76%），細胞膜、細胞壁生物合成（10.42%），細胞周期調控、細胞分裂與染色體重排（14.61%），碳水化合物轉運與代謝（11.55%），氨基酸轉運與代謝（7.51%）等功能類別的unigene 表達豐度也相對較高。

2.2.4 Unigene 的KEGG 代謝途徑分析共有30 941 個unigene 在KEGG 數據庫中注釋成功，涉及126 個代謝途徑（見圖3E）。與新陳代謝相關的unigene 數量最多，參與影響氨基酸代謝（1363 個）、次級代謝產物代謝（1085 個）、碳水化合物代謝（2418 個）、能量代謝（704 個）、脂質代謝（3308 個）及核酸代謝（1020 個）等多個方面，其中，參與聚糖合成與代謝的基因最少，僅有496 個，占unigene 總量的1.60%。

圖3 Unigene 功能注釋圖Fig 3 Functional annotation of unigenes

經統計分析，注釋到三酰甘油代謝路徑（ko00561）的unigene 有388 個，結合三酰甘油目前已知的生物合成途徑，篩選出257 個與三酰甘油生物合成途徑相關的unigene，并對其表達量進行分析（見表3）。二酰基甘油酰基轉移酶（diacylglycerolO-acyltransferase，DGAT） 是依賴酰基輔酶A 的三酰甘油從頭合成途徑的關鍵酶和限速酶。分別有19、5、1 個unigene 編碼DGAT1、DGAT2、DGAT3，其中編碼DGAT3 的unigene 表達量最高，FPKM 值為30.20。磷脂：二酰基甘油酰基轉移酶（phospholipid：diacylglycerol acyltransferase，PDAT）是不依賴酰基輔酶A 的合成途徑的關鍵酶，共有10 個unigene 編碼PDAT，FPKM 值最高為35.49。除此之外與三酰甘油生物合成的途徑相關的酶中表達量較高的還有醛脫氫酶（NAD＋）[aldehyde dehydrogenase（NAD＋），ALDH]、1-酰基-sn-甘油-3-磷酸酰基轉移酶（1-acyl-sn-glycerol-3-phosphate acyltransferase，plsC）和酰基甘油脂肪酶（acylglycerol lipase，MGLL），其FPKM 分別為993.54、113.86、70.14。

表3 薏苡三酰甘油生物合成中編碼關鍵酶的unigene Tab 3 Unigenes encoding the key enzymes involved in the synthesis of triacylglycerol in Coix

2.3 CDS 預測

CDS 預測可為薏苡基因組圖譜繪制及基因功能研究提供關鍵依據。使用BLAST 軟件比對得到48 449 個CDS 序列（63.61%），通過ESTScan 軟件預測得到2711 個CDS 序列（3.56%）。總共獲得51 160 個CDS 序列，占unigene 總數的67.17%。CDS 長度分布如圖4A 所示，長度在200～800 bp 內的CDS 序列有33 550 個（69.25%）；長度為300 bp 的數量最多，9793 個（20.21%）；其中與谷蛋白相關的CDS 序列有9 個（見表4），分別來自Swissprot、NR 數據庫和NT 數據庫，可用于后續谷蛋白基因克隆和定點突變研究。

表4 薏苡谷蛋白相關CDS 序列Tab 4 Glutelin related CDS of Coix

2.4 SSR 檢測結果

使用MISA 軟件對unigene 進行SSR 檢測，共搜索得到分布于11 101 個unigene 中的3324 個SSR，占unigene 總序列的29.94%。SSR 長度特征見表5及圖4B，SSR 類型豐富，從單核苷酸到六核苷酸均有出現，其中，三核苷酸重復最多（7899，占71.2%），出現頻率最高的基序為CCG/CGG（3690），最低的是ACT/AGT（84）。在二核苷酸重復類型中，AG/CT 基序出現的頻率最高（1407），而CG/CG 出現的頻率最低（242）。

圖4 薏苡CDS（A）和SSR（B）分布圖Fig 4 CDS（A）and SSR（B）distribution in Coix

表5 SSR 長度統計結果Tab 5 Statistics of SSR length distribution

2.5 SNP 位點分析結果

使用GATK 軟件對樣本基因進行SNP 檢測，共鑒定出22 764 個SNP，其中轉換15 141個（66.51%），發生A-G 替換的有7704 個，發生C-T 替換的有7437 個；變異類型為顛換的共7623 個（33.49%），發生A-C、A-T、C-G 及G-T替換的數目分別為1831、1699、2322 及1771 個（見表6）。

表6 SNP 變異類型統計Tab 6 Statistics of SNP variation types

2.6 Unigene 表達量計算

根據組裝結果，使用 Bowtie2 軟件將Clean reads 比對到unigene，后使用RSEM 軟件計算樣品的基因表達水平，計算獲得72 188 條unigene的表達量。

3 討論

RNA-seq 技術具有靈敏度高、通量高、成本低、操作簡單、重復性好、檢測范圍廣、分辨率高等諸多優勢，不局限于已知的基因組序列信息，還可用于發現未知基因及非編碼RNAs，精確地識別可變剪切位點以及SNP，提供更加全面的樣本轉錄組信息[3，28-29]。目前，RNA-seq 技術已成功應用于多種藥用植物轉錄組信息的采集。李依民等[6]應用Illumina HiSeq 2000 測序平臺，以PE150 的測序讀長，對黃三七根莖進行RNA-seq 技術，得到63 322 086 條Clean reads 和52 575 個unigene，為今后黃三七功能基因的鑒定、代謝途徑的分析及奠定基礎。Majeed 等[8]構建了紅豆杉參考轉錄組，經注釋得到129 869 個unigene，同時產生了7041 個SSR 標記，為破譯紅豆杉遺傳多樣性及其保護提供了有價值的信息。Amiripour 等[9]對阿育魏實的花序進行RNA-seq 技術，De novo組裝后產生68 051 個unigene，其中43 156 個獲得注釋，生成了8751 個SSR 標記，是目前第一個可用于轉錄組表征、基因發現、開發SSR 分子標記和阿育魏實基因組研究的數據庫。鄧楠等[30]應用Illumina HiSeq 2000 測序平臺對膜果麻黃不同發育期種子進行RNA-seq 技術，得到12 999 122 條初始序列，組裝后獲得49 449 條編碼基因，并根據KEGG 通路找到了有關姜醇烷、二芳基庚及芪類合成途徑的編碼基因片段230 個。張紹鵬等[31]獲得了珠子參根莖的轉錄組信息數據，為該物種資源的有效開發和利用提供了重要的理論基礎。李鐵柱等[32]利用高通量測序技術構建杜仲葉片和果實的轉錄組數據庫，共獲得54 471 338 條Raw reads，拼接組裝后得到49 610 條unigene，總長度為37 616 729 bp。本研究使用Illumina HiSeq 4000平臺對薏苡不同發育時期的各種組織部位RNA進行RNA-seq 技術，構建薏苡轉錄組數據庫。總共獲得76 164 個unigene，總長度為77 917 702 bp，數據量大且組裝效果較好，可為薏苡分子生物學研究及相關藥用價值的開發提供重要依據。將unigene 與七大功能數據庫比對后分別有49 092（NR：64.46%）、56 121（NT：73.68%）、32 748（Swissprot：43.00%）、20 670（COG：27.14%）、30 941（KEGG：40.62%），27 407（GO：35.98%）以及 29 789（Interpro：39.11%）個unigene 得到注釋。未得到注釋的unigene 占21.93%，推測原因可能是由于數據庫基因信息缺乏、測序技術的局限、序列長度較短未包含蛋白質功能域，或者是由于序列本身是非編碼序列或不完整序列或薏苡特有的新基因[33]。NR 數據庫比對結果顯示，薏苡與高粱具有最近的親緣關系（51.73%），這與Ottoboni等[34]研究結果一致。GO 功能分類中涉及代謝過程的unigene 數量最多（15 060）與KEGG 數據庫的注釋結果相一致。由KEGG 注釋結果可知，代謝分支中，除涉及全局和概覽圖途徑外，涉及脂質代謝途徑的unigene 數量最多（3308），提示脂類化合物在薏苡生長發育過程中有著重要作用。本研究分析了三酰甘油生物合成的相關代謝通路的257 個unigene，為今后提高三酰甘油含量奠定了理論基礎。

質譜鑒定中的數據庫搜索是肽/蛋白質結構解析的重要步驟，通過擴展蛋白數據庫可顯著提高肽/蛋白質鑒定的準確性和靈敏度。Goecks 等[35]采用RNA-seq 技術建立了果蠅雌性寄生蜂腹部轉錄本序列文庫，液相色譜-串聯質譜法鑒定寄生蜂毒液腺腔分離純化的肽序列，將肽譜數據與轉錄組數據庫信息進行比對，最終從寄生蜂L.boulardi毒液中鑒定獲得了129 種蛋白，從L.heterotoma毒液中鑒定得到176 種蛋白。本研究共獲得51 160 個CDS 序列，成功構建了薏苡蛋白數據庫，解決了以往薏苡蛋白序列信息量較少的問題，豐富了禾本科植物蛋白數據庫，為薏苡及其近緣物種蛋白及多肽序列的分析、鑒定提供重要依據。薏苡仁谷蛋白水解物具有顯著的降壓活性，本研究中獲得的CDS 序列中有9 個與谷蛋白相關，借助分子標記輔助育種，選擇高表達谷蛋白相關基因親本材料，可為薏苡抗高血壓親本的選育，提高薏苡仁的藥用價值奠定了理論基礎。推測到的3324 個SSR中CCG/CGG 基序出現最多（36 090），與水稻、高粱、玉米的SSR 檢測結果相一致[36]。基于這3324個SSR 位點，開發SSR 分子標記，將為薏苡的遺傳多樣性研究和進化研究奠定基礎。得到的22 764個SNP 位點中，轉換位點以A-G 居多（7704），顛換位點最多的是C-G（2322），與張得芳等[37]對花葉海棠的SNP 位點研究相一致。薏苡生長發育的過程中需要大量基因的調控，而SNP 位點的改變可能會導致所在基因序列的改變，從而對基因的轉錄、翻譯過程造成影響，得到結構和功能可能發生變化的蛋白質，進而造成薏苡植株性狀的改變[38]。因此，通過對SNP 位點和其改變所造成的植株性狀變化，可為薏苡的分子育種提供理論參考。