千針萬線草根轉錄組測序及黃酮類化合物合成相關基因挖掘

孫詩瑤，王曉麗，曹子林 ，張博林，郭紫微，趙思懿

（1.西南林業大學生態與環境學院，云南 昆明 650224；2.西南林業大學林學院，云南 昆明 650224）

0 引言

【研究意義】千針萬線草（Stellaria yunnanensis）為石竹科（Caryophyllaceae）繁縷屬（Stellaria）的一種多年生草本植物，其根具有補氣健脾、養肝活血的藥用價值[1]。據研究，繁縷屬植物普遍含有大量的黃酮類化合物，從中分離出的化合物主要有黃酮類、異黃酮類和黃酮醇類物質3類，具廣譜抗病毒藥物的功效[2]。因此，研究千針萬線草根的轉錄組和黃酮類代謝，對其開發利用具重要意義?！厩叭搜芯窟M展】譚瑞璞等[3-4]優化了千針萬線草中總黃酮的提取工藝，并測定不同產地千針萬線草中牡荊素的含量。在模式植物中，與類黃酮合成相關的功能基因已被克隆并進行廣泛的研究[5]。目前，黃酮類化合物的生物合成路徑已較為明確[6]。而就缺乏基因組信息的物種而言，對于植物次生代謝中類黃酮合成途徑的研究，多停留在紅花、中麻黃、半枝蓮、金線蓮、鐵皮石斛、藤茶等藥用植物的基因挖掘和鑒定上[7-12]?！颈狙芯壳腥朦c】千針萬線草基因組和轉錄組數據缺乏，制約了該物種種質資源鑒定、遺傳多樣性分析和分子機制等研究的開展。因此，對千針萬線草根進行高通量轉錄組測序及基因功能注釋，可為后續研究奠定基礎?！緮M解決的關鍵問題】利用Illumina Novaseq 6000高通量測序平臺對千針萬線草根進行轉錄組測序，應用公共數據庫對其基因組進行功能注釋及分類，借助KEGG挖掘出該物種黃酮類化合物合成相關基因，為開展黃酮類化合物合成的調控機制研究提供依據。

1 材料與方法

1.1 樣品采集

千針萬線草（Stellaria yunnanensis）根采自會澤縣迤車鎮箐口村山坡上，埋于土中帶回學校，種于花盆放在苗圃，每天正常澆水萌出根蘗苗后，經西南林業大學曾覺民教授鑒定為千針萬線草。于2021年7月取其根部以無菌水洗凈，濾紙吸去水分后裝入凍存管，液氮速凍、備用。

1.2 RNA提取、文庫的構建及測序

提取根的total RNA后，檢測RNA的濃度、純度和完整性。參照文獻[13]的方法，將帶有Oligo（dT）的磁珠和ployA配對A-T堿基。mRNA從total RNA中分離后，富集得到完整的RNA序列，添加fragmentation buffer隨機斷裂mRNA，利用磁珠篩選分離獲得約300 bp小片段；加入六堿基隨機引物，在逆轉錄酶存在下，以mRNA作為模板，反轉錄合成一鏈cDNA，再進行二鏈合成，進而形成穩定的雙鏈結構。連接adaptor，對目的條帶大小合適的片段用2%的瓊脂糖膠回收后，再用PCR技術進行擴增；最后，對短序列片段用Illumina Novaseq 6000進行測序。

1.3 測序原始數據過濾及組裝

測序得到原始數據后，用軟件Seq Prep將數據中帶有的大量接頭序列、低質量序列末端、adapter以及N的比例高于10%的堿基去除，獲得高質量的clean reads。使用Trinity軟件[14]組裝denove，用TransRate[15]評估從頭拼接結果中普遍存在的錯誤，獲取每一條contig的質量評分，并將這些contigs評分整合后，得到整個組裝結果的分數，篩選優化轉錄組獲得的序列。利用CD-HIT軟件[16]對序列比對聚類，剔除冗余、相似序列，最終獲得非冗余的序列。用BUSCO軟件[17]評價基因組或者轉錄組的組裝完整性。

1.4 基因功能注釋

利用NCBI將Unigene序列比對到Nr、Pfam、SwissProt和 egg NOG數據庫中（evalue＜0.000 01），用BLAST2 GO軟件對Unigene序列進行GO注釋，利用KOBAST軟件對Unigene序列進行KEGG注釋[18]，并分類統計注釋結果。用TF數據庫對預測的蛋白序列進行hmmscan比對并得到轉錄因子家族。

2 結果與分析

2.1 千針萬線草根轉錄組組裝質量分析

對千針萬線草根進行轉錄組測序后，經過嚴格的質量控制及數據過濾得到55 400 486個高質量且干凈的序列。用Trinity軟件將這些高質量數據進行de novo組裝，共獲得37 331 463條有效序列片段，N50為1 799 bp，平均長度為1 093.58 bp，GC核苷酸的含量（GC%）為41.02%，千針萬線草根Unigene34 137條（不含N的組裝片段）長度分布于201～14 683 bp（表1）。用長度分布柱狀圖分析組裝后的基因如圖1所示，長度為200～500 bp、501～1 000 bp、1 001～1 500 bp、1 501～2 000 bp、2 001～2 500 bp、2 501～3 000 bp、3 001～3 500 bp、3 501～4 000 bp、4 001～4 500 bp、大于4 500 bp的unigene分別占40.2%、20.2%、13.3%、10.3%、6.7%、3.8%、2.2%、1.3%、0.8%、1.2%。

圖1 千針萬線草根轉錄組Unigenes的長度分布Fig.1 Distribution of length of unigenes from S.yunnanensis roots

表1 Unigene基本信息表Table 1 Basic information on unigenes

2.2 千針萬線草根基因的總體注釋情況

將組裝得到的基因在Nr、Swiss-Prot、Pfam、egg NOG、GO、KEGG數據庫中逐一進行比對，所得到注釋的基因條數分別為21 510、17 208、17 702、19 414、17 191、10 285條。至少有22 369條unigenes在一個數據庫中進行了匹配和功能注釋，有11 768條unigenes未比對到數據庫，這表明功能未知的基因還較多。相比之下，注釋到Nr數據庫的基因數最多有21 510條，占比為63.01%；注釋到KEGG數據庫的基因最少，共有10 285條，占比為30.13%（表2）。

表2 Unigene注釋統計Table 2 Annotation of unigenes

2.3 千針萬線草根基因的NR功能注釋

共有21 510條unigenes比對到NR數據庫，與其他物種比對后的結果如圖2所示。注釋序列物種中注釋到基因數量最多的是甜菜（Beta vulgaris），共有6 473條，占比為30.09%；其余依次為藜麥（Chenopodium quinoa）、菠菜（Spinacia oleracea），分別有 6 124、3 954條unigenes，占比分別為28.47%、18.38%，這3種植物占比高達76.94%；其他297種物共有4 959條unigenes，占比23.06%。

圖2 千針萬線草根基因比對物種分類Fig.2 Species classification of unigenes

2.4 egg NOG功能注釋

將千針萬線草根基因比對到egg NOG數據庫中，得到注釋到的Unigenes有19 980條，占總Unigenes的58.53%，被分為23類（圖3）。有11 394個未知的功能注釋信息，其準確的生物學功能未確定，占所有功能注釋信息的57.03%；翻譯后修飾、蛋白開關和分子伴侶注釋到1 264條unigenes，所占比例為6.33%；轉錄注釋到1 144條unigenes，所占比例為5.73%；信號傳導機制注釋到895條unigenes，占4.48%；胞內運輸、分泌和囊泡運輸注釋到702條unigenes，所占比例為3.51%；碳水化合物運輸和代謝注釋到632條unigenes，所占比例為3.16%；翻譯、核糖體結構與生物發生注釋到615條unigenes，所占比例為3.08%；最少的是細胞運動和細胞核結構的unigenes，分別僅有6條（0.03%）和4條（0.02%）。這些結果表明，千針萬線草在翻譯后修飾、蛋白開關和分子伴侶、信號轉導機制、胞內運輸、分泌和囊泡運輸、碳水化合物運輸和代謝、翻譯、核糖體結構與生物發生等基因表達豐度較高。

圖3 千針萬線草unigene的egg NOG分類Fig.3 Egg NOG function classification of S.yunnanensis

2.5 GO功能注釋

注釋到GO數據庫中的unigenes一共有17 191條，獲得69 356個功能注釋，共分為生物過程、細胞組分和分子功能三類，依次得到21 481、27 214和20 661個功能注釋，分別占比為30.97%、39.24%、29.79%。3個本體細分為53個功能亞類，分子功能類包括16個功能亞類，結合（8 993個注釋）所占比例最多，其次是催化活性（8 393個注釋）；生物學過程類占其中23個功能亞類，細胞進程（6 958個注釋）和代謝過程（5 972個注釋）占比較大，其次是生物調控（2 581個注釋）；細胞組分類包括14個亞類，“細胞”所占比例最多（8 164個注釋），其次是“膜”（5 897個注釋）和“細胞器”（4 792個注釋）。只有少數轉錄本被注釋參與細胞殺傷、毒素活性、碳利用、介導分子活動、細胞增殖及發信號。這些結果表明，以結合、催化活性、細胞、細胞過程、代謝過程相關的基因較多（圖4）。

2.6 KEGG功能注釋

在KEGG數據庫共注釋到6 505個unigenes，可歸為6類，與代謝有關的通路共100條、注釋到的unigenes有3 190條，所占比例最高（49.04%）；與遺傳信息處理有關的通路有21條，注釋到的unigenes有2 115條（占32.51%）；環境信息處理相關的通路有4條，注釋到的unigenes有409條（占6.29%）；細胞過程相關的通路有6條，注釋到的unigenes有501條（占7.70%）；生物體系統相關的通路有2條，注釋到的unigenes有259條（占3.98%）；而人類疾病相關的通路有2條，注釋到的unigenes有31條（占0.48%）。代謝相關的通路細致分成10個亞類，其中排名前三的是碳水化合物代謝、氨基酸代謝、脂類物質代謝，分別占比為11.79%、7.29%、6.41%，遺傳信息加工分為4個亞類，環境信息處理和人類疾病各分為2個亞類，細胞過程和生物體系統的通路各分為1個亞類。其中“翻譯”“折疊、分類和降解”“運輸和代謝”“能量代謝”以及“信號轉導”分別占14.05%、9.99%、7.70%、5.70%、5.21%，其余12個亞類所占比例均小于5%（圖5）。

2.7 千針萬線草根中黃酮類物質生物合成相關基因挖掘

基于KEGG信號通路富集分析結果如表3所示，在千針萬線草根的轉錄組數據中共篩選獲得80條與黃酮類生物代謝相關的基因。其中 18條unigene 映射到了苯丙烷代謝通路（ko00940），45條unigene映射到類黃酮生物合成代謝通路上（ko00941），10條unigene映射到黃酮醇生物合成代謝通路（ko00944），7條unigene映射到異黃酮生物合成代謝通路上（ko00943）。

表3 千針萬線草根轉錄組中黃酮類化合物合成相關基因Table 3 Flavonoid biosynthesis-related gene in transcriptome of S.yunnanensis roots

在苯丙烷代謝通路中，18條unigene分別編碼苯丙烷途徑中的3個關鍵酶：PAL、C4H、4CL。其中5條unigene編碼該途徑中第1個催化酶PAL，5條unigene編碼第2個催化反應酶C4H，該酶屬于細胞色素P450（CY450）家族。8條unigene編碼該途徑的第3個酶4CL，其控制苯丙烷類代謝走向不同的代謝途徑，是一個重要分支點。在類黃酮生物合成代謝通路中，45條unigene共編碼10種該通路下的關鍵酶，包括查爾酮合成酶（CHS）、查耳酮異構酶（CHI）、黃烷酮-3-羥化酶（F3H）、肉桂酸-4-單加氧酶（CYP73A）、類黃酮-3′單加氧酶（CYP75B1）、根皮苷合酶（PGT1）、花青素還原酶（ANR）、莽草酸羥基肉桂轉移酶（HCT）、香豆酰脂3′羥化酶（C3′H）、咖啡酰輔酶A-O甲基轉移酶（CCOAOMT）。其中，CHS是催化黃酮類生物合成的第一個限速步驟，是植物聚酮合成酶超家族的成員之一，是植物類黃酮合成途徑中的第一個關鍵結構。CHI制約合成黃酮醇的含量。F3H是黃烷酮分支點的一個核心酶。在黃酮醇生物合成代謝通路中，分別編碼該通路下的2種關鍵酶：類黃酮-3′單加氧酶（CYP75B1）和黃酮醇-3-0-葡萄糖苷葡萄糖基轉移酶（FG3）。在異黃酮生物合成代謝通路中，7 條 unigene均編碼異黃酮2′-羥化酶（CYP81E）。

2.8 轉錄組Unigene的轉錄因子預測

千針萬線草轉錄組共預測到724個轉錄因子，可分為33個轉錄因子大家族。其中，MYB_superfamily（116個，16.02%）屬于最大家族，其次是ERF（59個，8.15%），接著為bHLH（57個，7.87%）、C2C2（53個，7.32%）、WRKY（52個，7.18%）、bZIP（47個，6.49%）、NAC（39個，5.39%）、C2H2（35個，4.83%）、B3-superfamily（31個，4.28%）、GRAS（31個，4.28%）、LBD（30個，4.14%）、C3H（29個，4.01%），其他21個家族共有基因145個，占20.03%（圖6）。

圖6 千針萬線草根Unigene的轉錄因子預測Fig.6 Predicted transcription factors of unigenes of S.yunnanensis roots

3 討論

開展藥用植物轉錄組研究是綜合詮釋植物基因組信息、發掘植物次生代謝生物合成途徑中關鍵酶以及植物次生代謝成分生物合成途徑與調控的新途徑[19]。本研究首次對藥用植物千針萬線草根進行轉錄組測序，測序產生clean reads的Q20值為98.42%，Q30值為95%以上，說明測序質量合格。經denove拼接組裝，去除冗余后，共得到34 137條Unigene，N50為1 799 bp，平均長度為1 093.58 bp。組裝結果與杭白芷（Angelica dahurica）根（N50為1 703 bp，平均長度1 164 bp）較為接近[20]。N50是評價組裝序列完整性的重要指標，本研究中N50＞800 bp，表明組裝片段完整，所得Unigenes 數量基本涵蓋全部轉錄信息，序列信息量龐大，可以滿足轉錄組分析的要求。

甜菜（Beta vulgaris）、藜麥（Chenopodium quinoa）、菠菜（Spinacia oleracea）同屬黎科，在Nr數據庫注釋中，注釋到基因占比高達76.94%，可能黎科與石竹科在親緣關系上比較接近。而石竹科蠅子草屬的滇白前注釋到黎科這三種植物基因占比達64.06%[21]，千針萬線草相比于滇白前，在與黎科親緣關系上可能較近。在GO功能注釋方面，可以分為3個大類和53個亞類，其中結合、催化活性、細胞過程及代謝過程的Unigene數目較多。通過egg NOG功能分類可知，參與翻譯后修飾、蛋白開關和分子伴侶、轉錄、信號傳導機制的Unigene較多，但有11 394條未知功能基因，占57.03%。這是由于egg NOG更新比較快，收錄來源很廣泛，數據庫收錄的未知序列占比就高，所以注釋的結果未知序列占比較高。轉錄因子是能從DNA上讀取并解釋遺傳“藍圖”的蛋白質組之一，它和DNA相結合，有助于啟動調節基因轉錄的程序，對植物應對非生物逆境脅迫具有重要影響[22]。研究表明，MYB、AP2/ERF、bHLH、bZIP、NAC、WRKY和C2H2轉錄因子對植物應對干旱、高鹽和寒冷脅迫均有顯著影響[21]。對這七類轉錄因子進行統計，千針萬線草根中有405個（55.94%），可能與其抗逆適應相關。在野外觀察發現，千針萬線草在向陽干旱貧瘠的陡坡也能正常生長，這是其抗逆性強的一個方面，是否能適應其他逆境，還有待進一步研究。

通過KEGG數據庫和通路分析發現，涉及代謝途徑的Unigene（3 190個）最多，通路數量最多的也是涉及代謝途徑的，共100條。從本研究獲得的基因注釋結果來看，有5個生物合成代謝途徑80個unigene涉及到藥用成分的合成，其中包括黃酮、黃酮醇、異黃酮等物質。通過進一步篩選，這些基因分別編碼16種關鍵酶，包括苯丙烷代謝途徑所需的3種酶（PAL、C4H、4CL），類黃酮代謝所需的10種關鍵酶（CHS、CHI、F3H、CYP73A、CYP75B1、PGT1、ANR、HCT、C3′H、CCOAOMT）和黃酮醇代謝途徑所需的2種關鍵酶（CYP75B1、FG3）和異黃酮代謝途徑關鍵酶CYP81E。千針萬線草根的藥用成分主要是黃酮類物質，通過對千針萬線草根轉錄組測序、注釋、分析以及黃酮類化合物合成相關基因的挖掘，豐富了千針萬線草根的遺傳信息，可以為進一步研究該物種藥用成分合成的關鍵基因及調控機制、克隆關鍵基因及功能分析奠定基礎，也為千針萬線草品質的形成提供理論依據。