鴉膽子轉錄組測序及基因功能注釋

石好宇，王慧蕓，趙藝，朱甫臻，尹艷燕，繆劍華，王靈芝*（.北京中醫藥大學生命科學學院，北京 0009；.廣西藥用植物園，廣西藥用植物資源保護與遺傳改良重點實驗室，南寧 53003）

中藥鴉膽子是苦木科植物鴉膽子[

Brucea javanica

（L.）Merr.]的干燥成熟果實，又名鴉蛋子、老鴉膽、苦榛子、苦參子，主產于我國廣東、廣西、云南、福建、海南和臺灣等地，始記于清代《本草綱目拾遺》。鴉膽子味苦、性寒，有小毒，歸大腸、肝經，有清熱解毒、截瘧止痢等功效。鴉膽子含有多種化學成分，主要包括苦木素類、類固醇、三萜類、生物堿、木質素、黃酮類、甾體類和脂肪酸等?？嗄舅仡惢衔锸区f膽子主要的藥理活性物質，是一類由苦楝烷型三萜降解得到的高度氧合三萜及其苷類化合物，多為四環三萜及五環三萜內酯，是鴉膽子的代表性藥理活性成分，其中鴉膽子苦醇（brusatol）和鴉膽丁（bruceantin）是該藥材中代表性抗腫瘤成分；苦木素還具有抗炎和抗醫學原蟲等藥理活性。鴉膽子油是鴉膽子中的脂肪油，抗腫瘤療效確切且毒副作用小，被廣泛應用于各種腫瘤的輔助治療。鴉膽子中蛋白質含量為17.47%，鴉膽子多肽具有抗菌活性，本課題組前期研究發現鴉膽子球蛋白酶解物具有顯著的抑制腫瘤細胞增殖的活性。

目前關于鴉膽子的研究主要集中于新化學成分的發現及作用，對于其次生代謝物合成途徑的研究較少，缺乏鴉膽子基因和蛋白質序列信息，極大地限制了對鴉膽子生長發育研究以及功能分析和利用。轉錄組測序（RNA-seq）是利用高通量測序技術將細胞或組織中mRNA、small RNA或no-codingRNA 進行測序分析的技術，可用于發現新基因、提供基因表達和調控信息。該技術已廣泛用于藥用植物，如半夏、黃芪以及薏苡仁等研究中。因此本文對鴉膽子進行轉錄組測序，并分析苦木素類化合物萜類骨架生物合成以及其衍生化修飾途徑，為今后藥用成分開發和藥用植物的遺傳改良提供理論支撐。

1 材料

1.1 藥材

鴉膽子采集于廣西南寧市武鳴區（108.27719°N，23.15643°E），采集不同發育時期的根、莖、葉、花序及果實等組織部位，液氮速凍后，保存于－80℃備用。

1.2 儀器與試藥

Aglient RNA 6000 Nano Reagents 試劑盒、Agilent 2100 生物分析儀（美國Aglient 公司）；BGISEQ-500 測序儀（深圳華大智造科技有限公司）；NanoDrop 超微量分光光度計（美國Nanodrop 公司）。

2 方法

2.1 RNA 提取與轉錄組文庫構建

利用Aglient RNA 6000 試劑盒提取鴉膽子各組織總RNA。采用NanoDrop 超微量分光光度計進行樣本RNA 純度檢測，通過Agilent 2100生物分析儀進行樣本濃度及完整性RIN（RNA integrity number）進行檢測。質檢合格后，等量混合樣品，用帶有Oligo（dT）的磁珠富集mRNA，采用fragmentation buffer 將mRNA 片段化，隨后使用隨機N6 引物進行反轉錄合成雙鏈cDNA，然后將雙鏈DNA 黏性末端補平、5'端磷酸化、添加并連接測序接頭，經PCR 擴增后構建cDNA 文庫，質檢合格的文庫使用BGISEQ-500進行測序。

2.2 轉錄組測序分析及功能注釋

測序原始序列（raw reads）使用SOAPunke軟件進行過濾，去除接頭污染、未知堿基N 含量大于5%、低質量的reads（質量低于10 的堿基占總堿基數20%以上），獲得clean reads。由于鴉膽子基因組測序尚未完成，缺少參考序列，因此使用Trinity 軟件對clean reads 進行

de novo

組裝，將組裝獲得的轉錄本去除冗余，得到非重復序列基因（Unigene）用于后續分析。使用單拷貝直系同源數據庫BUSCO（benchmarking universal single-copy orthologs）對組裝的轉錄本進行質量評估。

組裝得到轉錄組數據后，使用Blast 軟件將Unigene 序列與NR（non-redundant protein sequence database）、NT（nucleotide sequence database）、Swissprot、KEGG（kyoto encyclopedia of genes and genomes）和KOG（clusters of orthologous groups for eukaryotic complete genomes）數據庫進行比對，進行功能注釋?；贜R 注釋，使用Blast2GO 軟件對Unigene 進行GO（gene ontology）注釋，使用Hmmscan 軟件進行Pfam 注釋。

2.3 生物信息學分析

Getorf 軟件檢測Unigene 的開放閱讀框（ORF）后使用Hmmsearch 將ORF 比對到轉錄因子蛋白結構域（數據來源于TF），然后根據PlantTFDB 描述的轉錄因子家族特征對Unigene進行能力鑒定。使用TransDecoder 軟件識別Unigene 中的候選編碼區域，通過Blast 比對SwissProt 數據庫和Hmmscan 搜索Pfam 蛋白同源序列，從而進行編碼序列（CDS）預測。經MISA軟件對Unigene 進行分析后使用Primer3 對檢測到的簡單重復序列（SSR）進行引物設計。

2.4 基因表達量計算

基因表達量用FPKM（fragments per kilobase of transcript per million mapped reads）表示，即每百萬reads 中來自比對到某一基因每千堿基長度的reads 數目。本實驗使用Bowtie2 軟件將clean reads 比對到基因序列上，然后使用RSEM 計算各個樣品的基因表達水平。計算公式如下：

FPKM ＝（cDNA fragments）/[mapped fragments（millions）/transcript length（kd）]

式中，cDNA fragments 表示比對到某一轉錄本上的片段數目，即雙端reads 數目；mapped fragments（Millions）表示比對到轉錄本上的片段總數，以10為單位；transcript length（kd）表示轉錄本的長度，以10個堿基為單位。

3 結果與分析

3.1 轉錄組測序與組裝

采用NanoDrop 超微量分光光度計和Agilent 2100 生物分析儀對各樣本RNA 進行質量檢測，結果RIN ＞5.8，基線平整，表明完整性較好，樣品質量合格。使用BGISEQ-500 平臺共獲得鴉膽子114.55 Mb Raw reads，過濾去除低質量、接頭污染以及未知堿基N 含量過高的reads后，獲得110.67 Mb（96.61%）clean reads，其中質量值大于20（Q20）的堿基數目占總堿基數目的98.73%，且無明顯AT 或GC 分離現象（見圖1A），表明測序質量較高，能滿足后續分析。

對clean reads 進行組裝獲得轉錄本，使用BUSCO 軟件進行質量評估，其完全匹配的BUSCOs占比98.02%，表明組裝效果良好。使用Tgicl 軟件對轉錄本進行聚類并去除冗余得到85 828 條Unigenes，平均長度為1775 bp，GC 含量為40.07%，N50 值為2608 bp，N70 值為1887 bp，N90 值1005 bp。對Unigenes 的長度進行統計（見圖1B），其中長度≥3000 bp 的序列數目最多，為14 663 條。

圖1 Clean reads 的堿基含量（A）與Unigenes 長度（B）分布圖Fig 1 Distribution of clean reads base content（A）and Unigenes length（B）

3.2 鴉膽子Unigenes 的功能注釋

為獲得更全面的基因功能信息，將Unigenes與NR、NT、Swissprot、KEGG、KOG、Pfam和GO 七大功能數據庫進行比對，共有68 806 條Unigenes 得到注釋，占總數的80.17%。如表1 所示，有67 087（78.16%）條Unigenes 在NR 數據庫匹配成功，數量最多，而在GO 數據庫中，獲得注釋的Unigenes 數量最少，僅有43 042 個，占50.15%；有25 734 條Unigens 在七大功能注釋庫同時得到注釋，占Unigenes 總數的29.98%。

表1 Unigenes 功能注釋表
Tab 1 Function annotation of Unigenes

數據庫數量占比/%NR67 08778.16 NT59 85069.73 Swissprot49 80158.02 KEGG52 41661.07 KOG53 75562.63 Pfam 52 21460.84 GO43 04250.15 Intersection25 73429.98

3.2.1 NR 數據庫對比分析 將Unigenes 與NR 數據庫進行相似序列匹配，結果如圖2A 所示，與鴉膽子Unigenes 匹配數最多的物種是甜橙（

Citrus sinensis

），約占47.07%。此外，克萊門柚（

Citrus clementina

）和可可樹（

Theobroma cacao

）同源性分別占27.7%和2.08%。

3.2.2 GO 數據庫功能分類注釋 將鴉膽子Unigenes 與GO 數據庫進行比對，共有43 042 條Unigenes 得到注釋，獲得了99 940 條對應關系，結果如圖2B 所示，分屬于生物學過程、細胞組分和分子功能三大類。在生物學過程分類中，與細胞過程（cellular process）相關的Unigenes 為11 377 條，所占比例最高；在細胞組分分類中，與膜部分（membrane part）相關的Unigenes 最多，為12 940 條；在分子功能分類中，注釋上的Unigenes 涉及抗氧化、結合、催化、分子功能調節、轉運、轉錄調節等多方面功能活性，其中，涉及到結合（binding）功能的Unigenes 有21 235條，所占比例最高；對應蛋白標簽（protein tag）功能的Unigenes 最少，為7 條。

3.2.3 KEGG 代謝通路分析 鴉膽子52 416 條Unigenes 被注釋到KEGG 數據庫的134 條代謝通路中（見圖2C），其中涉及包含細胞過程（cellular processes，2282 條）、環境信息加工（environmental information processing，3191 條）、遺傳信息加工（genetic information processing，10 076條）、新陳代謝（metabolism，29 789 條）及生物體系統（organismal systems，1775 條）五大類。其中，與新陳代謝相關的Unigenes 數量最多，占總數的63.23%，主要參與影響氨基酸代謝、糖代謝、脂代謝等多個方面。

對鴉膽子藥用活性成分合成相關信號通路進行了分析，涉及不飽和脂肪酸生物合成（Ko01040）、萜類骨架生物合成（Ko00900）、類固醇生物合成（Ko00100）相關Unigenes 數量分別為206、276、8 條。根據鴉膽子KEGG 注釋結果，篩選276 條基因參與萜類骨架生物合成，其中涉及甲羥戊酸（MVA）以及非甲羥戊酸（MEP）途徑的分別有68、83 條Unigenes（見表2）。萜類化合物的生物合成涉及甲羥戊酸途徑和脫氧木酮糖-5-磷酸途徑兩條生物合成途徑，羥甲基戊二酰 CoA 還原酶（HMGR）是萜類代謝過程中的關鍵調控位點，有8 條Unigenes 調控此位點，其FPKM 值為85.37。1-脫氧木糖-5-磷酸合酶（DXS）與1-脫氧木糖-5-磷酸還原異構酶（DXR）是MEP 途徑關鍵限速酶，分別有53 和4條基因編碼，其FPKM 值分別為121.74、64.03。

表2 鴉膽子萜類骨架生物合成中編碼關鍵酶的Unigenes
Tab 2 Unigenes encoding the key enzymes related with terpenoid backbone biosynthesis in (L.) Merr

合成途徑酶的名稱酶編號Unigenes 數目FPKM 值MVAacetyl-CoA acetyltransferase（ACAT）2.3.1.96124.71 MVAhydroxy methylglutaryl-CoA synthase（HMGS）2.3.3.10540.20 MVAhydroxymethylglutaryl-CoA reductase（HMGR）1.1.1.34885.37 MVAmevalonate kinase（MVK）2.7.1.36217.02 MVAphosphomevalonate kinase（PMK）2.7.4.2311.56 MVAmevalonate diphosphate decarboxylase（MVD）4.1.1.33333.39 MVAisopentenyl-diphosphate delta-isomerase（IPPs）5.3.3.211339.72 MVAfarnesyl diphosphate synthase（FDPS）2.5.1.130239.09 MEP1-deoxy-D-xylulose-5-phosphate synthase（DXS）2.2.1.753121.74 MEP1-deoxy-D-xylulose-5-phosphate reductoisomerase（DXR）1.1.1.267464.03 MEP2-C-methyl-D-erythritol 4-phosphate cytidylyltransferase（ispD）2.7.7.60820.36 MEP4-diphosphocytidyl-2-C-methyl-D-erythritol kinase（ispE）2.7.1.148548.69 MEP2-C-methyl-D-erythritol 2，4-cyclodiphosphate synthase（ispF）4.6.121120.13 MEP（E）-4-hydroxy-3-methylbut-2-enyl-diphosphate synthase（ispG）1.17.7.1188.64 MEP4-hydroxy-3-methylbut-2-en-1-yl diphosphate reductase（ispH）1.17.7.411203.29

萜類骨架形成后需要衍生修飾，細胞色素P450（cytochrome P450，CYP450）和糖基轉移酶（UDP-glycosyltransferase，UGT）主要起氧化、置換和糖基化修飾的重要作用。根據鴉膽子轉錄組SwissProt 數據庫注釋結果，有284 條 Unigenes被注釋為 CYP450，隸屬于20 個CYP450 家族，其中 CYP71 家族的Unigenes 為71 條，占比最多，其次是CYP82 和CYP704，分別為34 條和22 條。共找到屬于18 個UGT 亞家族的129 個UGTs，包括22 個 UGT74、17 個 UGT83 等。對各個樣品中基因表達量進行聚類表達分析，其中CYP450家族 FPKM 值較高的有 CYP82A1（FPKM ＝284.88）、CYP71A1（FPKM ＝112.44）和CYP71A2（FPKM ＝105.64），UGT 家族FPKM 值較高的有UGT73C3（FPKM ＝65.93）、UGT71A15（FPKM ＝53.73）和UGT71A16（FPKM ＝36.44）。

3.2.4 KOG 分類 利用Blast 軟件，共有53 755條Unigenes 被注釋到KOG 數據庫中，得到基因同源物的分類信息（見圖2D）。KOG 根據其功能可分為25 類，其中，一般功能預測（general function prediction only）類別中Unigenes 數量最多，占總量的21.94%；細胞運動（cell motility）類Unigenes 數量最少， 為1 個； 未知功能（function unknown），翻譯后修飾、蛋白質轉換、分子伴侶（posttranslational modification，protein turnover，chaperones）、信號轉導機制（signal transduction mechanisms）和轉錄（transcription）等類別的Unigenes 表達豐度也相對較高，分別占總量的9.32%、7.92%、10.78%和6.50%。

圖2 Unigene 功能注釋圖Fig 2 Functional annotation of Unigenes

3.3 CDS 預測

CDS 預測可為鴉膽子基因組圖譜繪制及基因功能研究提供關鍵依據，結果共預測得到56 624個CDS 序列，總長度為64 922 145 bp，N50 值為1446 bp，GC 含量為43.48%。課題組前期工作表明，鴉膽子球蛋白酶解物具有顯著的體外抗腫瘤活性，通過注釋，共獲得10 個與球蛋白相關的CDS序列，結果見表3，可為今后高球蛋白鴉膽子藥材的篩選及抗腫瘤生物肽的開發提供文本信息。

表3 鴉膽子球蛋白相關CDS 序列
Tab 3 Globulin related CDS of (L.) Merr

編號基因IDFPKM 值功能注釋數據庫來源1 Unigene 1078627.63Basic 7S globulinSwissprot 2 Unigene 123901.69Basic 7S globulinSwissprot 3 Unigene 1523759.38Basic 7S globulinSwissprot 4 Unigene 1564050.42Basic 7S globulinSwissprot 5 Unigene 20499.46Basic 7S globulinSwissprot 6 Unigene 2744610.94Basic 7S globulin 2Swissprot 7 Unigene 718641.50Basic 7S globulinSwissprot 8 Unigene1305220.96Citrus sinensis 12S seed storage globulin 1-likeNT 9 Unigene 108290.632S globulin [Corchorus olitorius]NR 10Unigene 308863.03basic 7S globulin-like [Manihot esculenta]NR

3.4 鴉膽子轉錄組SSR 分析

鴉膽子轉錄組共鑒定出6 種SSR 重復類型，共檢測出36 045 個SSR 分布于25 747 個Unigenes中，單核苷酸至六核苷酸重復類型均存在，其中有14 631 個單核苷酸，11 594 個二核苷酸，7554個三核苷酸，649 個四核苷酸，574 個五核苷酸，1043 個六核苷酸，單核苷酸重復類型數量最多。

4 討論

RNA-seq 技術能夠對真核生物復雜的轉錄本結構和表達水平進行分析，還可用于發現非編碼RNA 和全新的基因轉錄本，有助于揭示生物生長發育、逆境應答機制和次生代謝產物的富集調控機制。此外，該技術還在低豐度基因檢測、基因家族鑒定等方面發揮重要作用。有學者利用Illumina 高通量測序發現三江源地區灌木亞菊中與藥用活性相關的基因主要富集在甾體、黃酮類、苯丙素類及萜類化合物；通過與烏頭轉錄組測序分析比較，證實鐵棒錘是烏頭堿生物合成途徑的理想材料；梅瑜等利用高通量測序PacBio Sequel 平臺深入了解甘葛藤轉錄組的整體水平，挖掘其中黃酮類生物合成相關候選基因；林江波等利用Illumina HiSeq 4000 高通量測序技術對鐵皮石斛莖和葉進行轉錄組測序，并分析了植物甾醇生物合成關鍵酶基因的表達。

本研究首次采用BGISEQ-500 測序平臺，進行了鴉膽子轉錄組序列分析，共獲得85 828 條Unigenes，數據量大且組裝效果良好，填補了鴉膽子轉錄組信息的空白，確定了參與萜類骨架生物合成的基因，文本信息將有助于進一步研究苦木素類化合物及其衍生化在鴉膽子中生物合成的分子機制，并促進鴉膽子功能基因組學的研究。MEP 和MVA 通路在植物萜類化合物的生物合成中發揮著重要作用。在次級代謝物的衍生化修飾過程中，需要CYP450 和糖基轉移酶的參與。Zheng 等研究表明CYP450 參與了真菌甲基的羥基化反應；UGT73AD1 可使羥基積雪草酸的C-28 羧基部分糖基化。本研究發現鴉膽子轉錄組中有較多CYP450、UGT 序列，推測其可能在苦木素類衍生修飾中起重要作用。

本研究共獲得56 624 個CDS 序列，構建鴉膽子基因編碼數據庫，豐富了該物種蛋白質編碼基因的注釋，為鴉膽子及其近緣物種蛋白序列鑒定提供了更可靠的質譜數據庫信息。后期可通過蛋白質組學等手段來進行數據校正，并確定新的編碼基因；與蛋白質組學聯用還可以有效地量化生物體中重要的蛋白質，確定基因在植物生長發育過程中的功能。鴉膽子球蛋白及其酶解物具有明確的抗腫瘤活性，在獲得的CDS 序列中有10 個與球蛋白相關，今后可借助分子標記輔助育種技術，進行高球蛋白親本的篩選，增強鴉膽子球蛋白抗腫瘤活性。通過植物轉錄組測序，還可反向預測和開發多肽類藥物。Rodríguez-Decuadro 等根據

Peltophorum dubium

幼苗轉錄組測序結果，采用同源序列搜索和半胱氨酸殘基模式匹配的方法，對其中的抗菌肽進行了預測。

鴉膽子SSR 位點重復類型以單核苷酸和雙核苷酸為主，占72.7%，與蘋果、川穹等植物研究結果一致。短重復序列的大量存在表明該物種的進化水平相對較高，表明鴉膽子處于相對較高的生物進化與分類地位，其基因組經歷了較長的進化時間或具有較高的突變頻率，還可為近緣物種SSR 標記的開發及其遺傳分析提供便利。

5 結論

本研究基于高通量測序平臺進行了鴉膽子全長轉錄組分析，獲得了高質量轉錄本，共有68 806 條Unigenes 注釋到七大功能數據。萜類骨架是苦木素生物合成的重要前體，KEGG 分析共有276 條Unigenes 參與，其中涉及MEP 和MVA途徑的15 種酶；還發現涉及鴉萜類骨架修飾的CYP450 和UGT 的Unigenes 分別有284 條和129條，為提高苦木素類化合物產量的分子生藥學研究提供了重要文本信息。通過CDS 預測共獲得10 個CDS 序列與球蛋白相關，進一步豐富了鴉膽子球蛋白源抗腫瘤肽的數據庫。本研究獲得了鴉膽子轉錄組的數據，可為今后鴉膽子藥用價值開發與利用提供理論基礎。