余靜雅，夏銘澤，張發(fā)起

（1.中國科學(xué)院 西北高原生物研究所 高原生物適應(yīng)與進(jìn)化重點(diǎn)實(shí)驗(yàn)室，西寧 810001； 2.中國科學(xué)院大學(xué) 生命科學(xué)學(xué)院，北京 100049）

甘青蒿（Artemisia tanguticaPamp.）屬于菊科（Asteraceae）蒿屬（ArtemisiaL.）蒿亞屬（SubgenusArtiemisia）腺毛蒿組（SectionViscidipubes），廣泛分布于青藏高原東南邊緣中高海拔地區(qū)，為青藏高原特有種之一［1］，民間入藥具有多種功效［2］。蒿屬中多個(gè)物種被發(fā)現(xiàn)具有驅(qū)寒除濕、活血通經(jīng)的功效，可用于治療瘧疾和感冒發(fā)燒等疾?。?］。

蒿屬植物的揮發(fā)油具有抗菌、抗病毒和殺蟲的作用［4-6］，其次生代謝產(chǎn)物與各種藥理特性密切相關(guān)［7］。作為揮發(fā)油、樹脂和蠟的重要成分，許多萜類次生代謝產(chǎn)物同時(shí)也在植物生理過程中起著重要作用［8］。青蒿素作為一種內(nèi)過氧化物倍半萜內(nèi)酯（Endoperoxide sesquiterpene lactone），在蒿屬多個(gè)物種中均有存在［9］。甘青蒿整株密被腺毛或粘質(zhì)柔毛，葉片中含有大量的揮發(fā)油。然而，目前對于蒿屬植物的研究主要集中于生態(tài)功能、化學(xué)和藥理學(xué)研究［10-11］，其生物分子研究遠(yuǎn)遠(yuǎn)落后［4］，僅有部分蒿屬植物的葉綠體和轉(zhuǎn)錄組被測序［12-13］，以及黃花蒿（Artemisia annuaL.）的基因組被測序［14］。作為青藏高原廣布的蒿屬植物之一，由于地理分布和缺乏基因組序列，甘青蒿相關(guān)研究受到限制。對于缺乏生物分子信息的非模式物種，高通量RNA測序和生物信息學(xué)是挖掘特定生物功能基因的有力工具［15］。

二代測序技術(shù)具有通量高和準(zhǔn)確性較高的優(yōu)勢，被廣泛應(yīng)用，但是由于其讀長較短（＜ 300 bp），基于第二代高通量測序平臺(tái)往往不能準(zhǔn)確得到和組裝出完整轉(zhuǎn)錄本。相對于一代和二代測序技術(shù)，第三代測序技術(shù)具有超長讀長（平均讀長可達(dá)15 kb）和直接讀取目標(biāo)序列的優(yōu)勢［16］，可以彌補(bǔ)二代測序技術(shù)的不足，提高基因注釋的準(zhǔn)確性［17］。目前，全長轉(zhuǎn)錄組測序技術(shù)已廣泛應(yīng)用于動(dòng)植物研究中［18-20］，如利用PacBio測序技術(shù)揭示了艾蒿（Artemisia argyiL.）不同部位中萜類化合物生物合成相關(guān)基因的活性差異，并篩選出12種編碼重要酶的候選基因［12］；利用牛津納米孔技術(shù)（Oxford nanopore technology，ONT）獲得茶樹［Camellia sinensis（L.） O. Ktze.］不同部位的全長轉(zhuǎn)錄組，并檢測到茶葉風(fēng)味合成相關(guān)的途徑［21］；對花生（Arachis glabrataBenth.）的全長轉(zhuǎn)錄組分析揭示了其對干旱、寒冷和鹽等非生物脅迫的基因表達(dá)譜［19］；地黃［Rehmannia glutinosa（Gaert.） Libosch. ex Fisch. et Mey.］的全長轉(zhuǎn)錄組檢測到了毛蕊花糖苷生物合成途徑的關(guān)鍵基因［22］。目前尚無關(guān)于甘青蒿的類似研究，因此，本研究以甘青蒿葉片為材料，利用三代測序平臺(tái)ONT對全長轉(zhuǎn)錄組進(jìn)行測序及組裝，鑒定和注釋新轉(zhuǎn)錄本，尋找可能與蒿屬萜類化合物生物合成相關(guān)的候選基因，并對甘青蒿全長轉(zhuǎn)錄組的轉(zhuǎn)錄因子（Transcription factor，TF）進(jìn)行鑒定和分析，為甘青蒿生物資源利用提供一定的科學(xué)依據(jù)。

1 材料與方法

1.1 材料

用于測序的甘青蒿鮮嫩葉片采集于青海省海西蒙古族藏族自治州烏蘭縣（地理坐標(biāo)：N 36°56′37"，E 98°27′20"；海拔2 939 m），植株葉片采集后迅速放入液氮中保存，用于后續(xù)RNA提取。憑證標(biāo)本（chen2019104）存放于中國科學(xué)院西北高原生物研究所青藏高原生物標(biāo)本館（HNWP）。

1.2 方法

1.2.1 RNA提取與建庫測序

采用Trizol法［23］提取甘青蒿葉片總RNA，瓊脂糖凝膠電泳用于檢查提取的RNA是否存在污染，并通過Nanodrop檢測RNA濃度，使用NanoPhotometer spectrophotometer和安捷倫2 100生物分析儀（Agilent Technologies，Palo Alto，CA）評估RNA質(zhì)量和完整性。檢測合格的RNA樣品使用Oligo DT作為引物對目標(biāo)mRNA進(jìn)行反轉(zhuǎn)，通過低循環(huán)PCR擴(kuò)增全長cDNA。NEBNext End Repair/dA-Tailing Module被用于末端修讀及加A，ONT SQKLSK109試劑盒及NEBNext Quick Ligationg Module用于測序接頭的連接。測序平臺(tái)為PromethION（Oxford Nanopore Technologies公司，英國）。

1.2.2 全長轉(zhuǎn)錄組序列分析

對ONT原始測序數(shù)據(jù)進(jìn)行質(zhì)量控制，使用NanoFilt v2.8.0［24］過濾質(zhì)量小于7和長度小于500 bp的序列（參數(shù)：-q 7 -l 50），過濾得到的有效數(shù)據(jù)使用SeqKit v0.12.0［25］默認(rèn)參數(shù)進(jìn)行統(tǒng)計(jì)。已發(fā)表的黃花蒿全基因組作為參考，使用minimap2 v2.17-r941［26］將過濾后的全長序列與參考基因比對，參數(shù)設(shè)置為：-ax splice -uf -k 14，samtools v1.11［27］用于結(jié)果統(tǒng)計(jì)，參數(shù)為flagstat。使用Pinfish（https：//tracker.debian.org/pkg/pinfish，版本：0.1.0；參數(shù)：default）軟件對全長序列快速構(gòu)建非冗余轉(zhuǎn)錄本集，得到一致序列。將得到的一致序列與參考基因組進(jìn)行比對，使用StringTie v2.1.4［28］對比對結(jié)果去冗余，參數(shù)設(shè)置為：-conservative -L -R，合并僅5’端外顯子有差異的比對，得到非冗余轉(zhuǎn)錄本序列。Gffcompare v0.12.1［29］用于新轉(zhuǎn)錄本及新基因的預(yù)測，補(bǔ)充現(xiàn)有注釋。TransDecoder v5.5.0（https：//github.com/TransDecoder/TransDecoder）對新鑒定的轉(zhuǎn)錄本進(jìn)行蛋白質(zhì)編碼（Coding Sequence，CDS）序列預(yù)測，參數(shù)為：-m 50，-single_best_only。

為獲得轉(zhuǎn)錄本全面的功能信息，基于序列相似性和motif相似性，將新轉(zhuǎn)錄本注釋到公共數(shù)據(jù)庫，包括非冗余蛋白數(shù)據(jù)庫（Non-Redundant Protein atabase，NR）［30］、蛋白質(zhì)真核同源數(shù)據(jù)庫（Eukaryotic Orthologous Groups，KOG）［28］、非冗余蛋白質(zhì)序列數(shù)據(jù)庫（Universal Protein，uniprot）［31］、蛋白質(zhì)原核同源數(shù)據(jù)庫（Cluster of Orthologous Groups of Proteins，COG）［32］、基因本體論數(shù)據(jù)庫（Gene Ontology，GO）［33］、蛋白質(zhì)家族域數(shù)據(jù)庫（Protein Families Database，Pfam）［34］和東京基因與基金組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）［35］。

1.2.3 全長轉(zhuǎn)錄組轉(zhuǎn)錄因子分析

利用iTAK［36］結(jié)合數(shù)據(jù)庫PlantTFDB v5.0［37］實(shí)現(xiàn)轉(zhuǎn)錄因子的預(yù)測。

2 結(jié)果與分析

2.1 全長轉(zhuǎn)錄組測序與組裝

甘青蒿納米孔測序產(chǎn)出的raw data達(dá)到15.66 Gb。NanoFilt對原始數(shù)據(jù)過濾后獲得10.68 Gb的clean data，共得到14 730 754條序列，N50為999 bp，平均長度為779 bp。將過濾后的全長序列與參考基因組比對，共比對上9 557 431條序列，比對率達(dá)到81.6%。進(jìn)一步過濾冗余全長轉(zhuǎn)錄組序列，得到18 634條非冗余全長轉(zhuǎn)錄本，N50為1 185 bp，最大長度為5 211 bp，平均長度為842 bp。

將非冗余轉(zhuǎn)錄本與基因組一致轉(zhuǎn)錄本進(jìn)行比較，發(fā)現(xiàn)4 358條新轉(zhuǎn)錄本，主要為未知的新轉(zhuǎn)錄本（圖1A）。對新鑒定的轉(zhuǎn)錄本進(jìn)行CDS預(yù)測，獲得1 482條CDS序列，N50為816 bp，最大長度為3 492 bp，平均長度為695 bp，長度主要分布在300 ～ 1 000 bp（圖1B）。

圖1 新鑒定轉(zhuǎn)錄本Fig.1 Newly identified transcript

2.2 新鑒定轉(zhuǎn)錄本功能注釋

將新鑒定的轉(zhuǎn)錄本分別注釋到7個(gè)數(shù)據(jù)庫中，4 358條轉(zhuǎn)錄本中共有2 184條（50.11%）被成功注釋（表1）。其中共有2 172條轉(zhuǎn)錄本比對到NR數(shù)據(jù)庫的122個(gè)物種，比對率較高的前3個(gè)物種都屬于菊科，分別是黃花蒿（979條，48.51%）、除蟲菊［Tanacetum cinerariifolium（Treviranus） Schultz Bipontinus］（649條，32.17%）和向日葵（Helianthus annuusL.）（116條，5.71%）。

表1 新鑒定轉(zhuǎn)錄本注釋結(jié)果Tab. 1 The result of newly identified transcript annotation

選取NR、GO、KEGG、KOG和uniprot這5個(gè)常用的數(shù)據(jù)庫注釋情況繪制韋恩圖（圖2），915條轉(zhuǎn)錄本在5個(gè)數(shù)據(jù)庫中均有注釋。

圖2 新鑒定轉(zhuǎn)錄本功能注釋韋恩圖Fig.2 Venn diagram of function annotated newly identified transcript

2.2.1 GO注釋

GO數(shù)據(jù)庫中生物學(xué)過程、細(xì)胞組分和分子功能3個(gè)分支注釋到的新鑒定轉(zhuǎn)錄本的數(shù)量分別為1 064、1 968和1 322條，共注釋到955個(gè)功能組，選取每個(gè)分支下前20個(gè)注釋最多的功能組繪圖（圖3）。在生物學(xué)過程中注釋數(shù)量最多的是翻譯（82條）、碳水化合物代謝過程（31條）和甲基化（21條）。在細(xì)胞組分中注釋數(shù)量最多的是膜的組成（423條）、細(xì)胞核（134條）和細(xì)胞質(zhì)（100條）。在分子功能中注釋數(shù)量最多的是ATP結(jié)合（121條）、金屬離子結(jié)合（114條）和核糖體結(jié)構(gòu)組成部分（92條）。

圖3 新鑒定轉(zhuǎn)錄本GO注釋結(jié)果Fig. 3 GO annotation classification of newly identified transcript

2.2.2 KEGG注釋

KEGG數(shù)據(jù)庫注釋結(jié)果顯示，新鑒定轉(zhuǎn)錄本共注釋到5條主通路（細(xì)胞過程、環(huán)境信息過程、遺傳信息過程、代謝和組織系統(tǒng)）和18條子通路（表2），其中注釋數(shù)量最多的子通路為全局和概覽圖（673條）、碳水化合物代謝（161條）和翻譯（140條）。共有70條轉(zhuǎn)錄本被注釋到可能與萜類化合物代謝相關(guān)子通路中，其中29條轉(zhuǎn)錄本被注釋到萜類和多酮類化合物代謝子通路，41條轉(zhuǎn)錄本被注釋到輔助因子和維生素代謝子通路。

表2 KEGG注釋結(jié)果Tab. 2 KEGG annotation classification of A.tangutica

對注釋到可能與萜類化合物代謝相關(guān)子通路的轉(zhuǎn)錄本進(jìn)行進(jìn)一步篩選，去除不完整的轉(zhuǎn)錄本及注釋到不相關(guān)子通路上的基因（諸如卟啉和葉綠素代謝、煙酸和煙酰胺代謝、核黃素代謝等）。新鑒定轉(zhuǎn)錄本中共有9條完整轉(zhuǎn)錄本涉及萜類生物合成途徑（表3），可能是與甘青蒿葉片萜類化合物生物合成相關(guān)的候選基因，經(jīng)過序列比對之后，其中涉及類單萜生物合成的轉(zhuǎn)錄本可能編碼ctcad3和CTI12_AA295060基因。涉及萜類化合物骨架生物合成的轉(zhuǎn)錄本可能編碼CTI12_AA344030、MCS和hdr基因。涉及泛素酮和其他萜類醌的生物合成的轉(zhuǎn)錄本可能編碼CTI12_AA436250、CTI12_AA496620和CTI12_AA199400基因。

表3 萜類化合物生物合成相關(guān)基因Tab. 3 Genes related to terpenoid biosynthesis

2.3 轉(zhuǎn)錄因子分析

iTAK用于甘青蒿全長轉(zhuǎn)錄組序列的轉(zhuǎn)錄因子預(yù)測，結(jié)果獲得29 154個(gè)轉(zhuǎn)錄因子，分別屬于58個(gè)不同的轉(zhuǎn)錄因子家族，選擇數(shù)量最多的前20個(gè)轉(zhuǎn)錄因子家族進(jìn)行可視化（圖4），前3位分別是bHLH家族（2 936個(gè)，10.07%）、NAC（2 325個(gè)，7.97%）和MYB_related（1 841個(gè)，6.31%）。這些轉(zhuǎn)錄因子信息為甘青蒿次生代謝產(chǎn)物的生物合成和抗逆性研究提供了依據(jù)。

圖4 甘青蒿全長轉(zhuǎn)錄組轉(zhuǎn)錄因子家族分析Fig. 4 Analysis of transcription factor family in full-length transcriptome of A. tangutica

3 討論

本研究利用三代測序平臺(tái)ONT結(jié)合生物信息學(xué)分析，獲得了甘青蒿全長轉(zhuǎn)錄組數(shù)據(jù)。通過聚類、去冗余和校正后最終得到18 634條非冗余全長轉(zhuǎn)錄本，N50為1 185 bp。N50大于1 000 bp且GC含量穩(wěn)定，說明序列組裝完整性較好，能夠滿足后續(xù)生物信息學(xué)分析的要求。

本研究中，共有4 358條新轉(zhuǎn)錄本被鑒定，其中2 184條（50.11%）得到7個(gè)數(shù)據(jù)庫的注釋。相比于艾蒿全長轉(zhuǎn)錄組高達(dá)93.97%的注釋率［12］，甘青蒿注釋率低原因可能為該物種與黃花蒿及艾蒿親緣關(guān)系較遠(yuǎn)。此外，黃花蒿及艾蒿多分布于低海拔區(qū)域，而甘青蒿主要分布于青藏高原，由于生境不同，甘青蒿可能擁有獨(dú)特的生理過程。這一情況也反映出公共數(shù)據(jù)庫中菊科植物遺傳信息的匱乏。因此，本研究獲得的甘青蒿全長轉(zhuǎn)錄組能夠豐富菊科植物遺傳數(shù)據(jù)庫；為探究青藏高原植物獨(dú)特適應(yīng)性提供基礎(chǔ)數(shù)據(jù)；為挖掘相關(guān)功能基因提供參考依據(jù)。

萜類化合物能夠抑制多種細(xì)菌和蚜蟲，具有顯著的藥理活性，如免疫調(diào)節(jié)、抗腫瘤、降血脂等［38］，廣泛應(yīng)用于制藥及食品工業(yè)。蒿屬植物富含萜類化合物，對其萜類化合物進(jìn)行制取，可以提高蒿屬植物的經(jīng)濟(jì)價(jià)值［12］。通過將轉(zhuǎn)錄本注釋到KEGG數(shù)據(jù)庫，共篩選出9個(gè)可能參與萜類生物合成的完整候選基因，其中6個(gè)為新鑒定出的候選基因。2-C-甲基-D-赤蘚糖醇4-磷酸（2-C-methyl-D-erythritol 4-phosphate，MEP）途徑為生物合成萜類化合物的主要途徑之一。MCS基因和hdr基因參與MEP途徑，同時(shí)也是MEP途徑中的兩個(gè)關(guān)鍵酶［39］。其中，MCS基因則參與MEP途徑的第五步，該基因編碼2-C-甲基-D-赤蘚糖醇2，4-環(huán)二磷酸酯合成酶（2-C-methyl-D-erythritol 2，4-cyclodiphosphate synthase），負(fù)責(zé)催化2-C-甲基-D-赤蘚糖醇2，4-環(huán)二磷酸酯（2-C-methyl-D-erythritol 2，4-cyclodiphosphate）的合成，是二萜類化合物生物合成的關(guān)鍵基因之一［40］。hdr基因參與MEP途徑的最后一步，該基因編碼4-羥基-3-甲基丁-2-烯基二磷酸還原酶（4-hydroxy-3-methylbut-2-enyl diphosphate reductase），負(fù)責(zé)將4-羥基-3-甲基丁-2-烯基焦磷酸酯（4-Hydroxy-3-methylbut-2-enyl pyrophosphate）催化成焦磷酸異戊酯（Isopentenyl pyrophosphate）和二甲基烯丙基焦磷酸酯（Dimethyl allyl pyrophosphate）。研究表明該基因的過量表達(dá)不僅可以提高青蒿素的產(chǎn)量，也增加了黃花蒿中倍半萜類化合物和單萜類化合物的含量［41］。同時(shí)，該基因在艾蒿全長轉(zhuǎn)錄組分析中也被鑒定為萜類化合物生物合成的關(guān)鍵基因［12］。CTI12_AA199400基因編碼細(xì)胞色素P450s，主要參與萜類化合物生物合成過程中多種修飾酶的合成，在植物萜類生物合成中起著重要的作用［42］。然而，其他候選基因?qū)祁惢衔锷锖铣捎绊懙难芯縿t少見報(bào)道。

在此前的研究中，黃花蒿基因組鑒定出2 717個(gè)轉(zhuǎn)錄因子［14］，艾蒿全長轉(zhuǎn)錄組鑒定出5 604個(gè)轉(zhuǎn)錄因子［12］，本研究共鑒定出29 154個(gè)轉(zhuǎn)錄因子，豐富了蒿屬植物的轉(zhuǎn)錄因子信息。已發(fā)現(xiàn)AP2、bHLH、MYB、NAC、WRKY和bZIP轉(zhuǎn)錄因子家族通過激活或抑制萜類生物合成途徑中關(guān)鍵基因的表達(dá)來調(diào)節(jié)萜類代謝［43-44］。丹參（Salvia miltiorrhizaBunge）MYB轉(zhuǎn)錄因子能夠促進(jìn)相關(guān)基因的表達(dá)，增加丹參酮（Tanshinone）的產(chǎn)量［45］。黃花蒿bHLH轉(zhuǎn)錄因子通過激活關(guān)鍵酶的啟動(dòng)子來調(diào)節(jié)青蒿素的生物合成［46］，并參與植物的萜類化合物的生物合成［47］；WRKY［48］、HD-ZIP［49］和MYB［50］與青蒿素的合成密切相關(guān)。本研究中，同樣鑒定到了這些轉(zhuǎn)錄因子家族的表達(dá)，由此推測甘青蒿體內(nèi)可能存在青蒿素的生物合成，具體過程還需要進(jìn)一步研究。

青藏高原不同尋常的極端環(huán)境通常會(huì)誘導(dǎo)植物的各種抗逆反應(yīng)［51］。本研究涉及萜類生物合成的候選基因中，ctcad3編碼乙醇脫氫酶，其過量表達(dá)能夠提高植物對鹽、干旱、寒冷和病原體感染的抗逆性［52］。CTI12_AA295060基因編碼的葡萄糖/核糖醇脫氫酶作為一種醇-多糖脫氫酶（Alcohol-polyolsugar dehydrogenase），可能具有碳水化合物代謝和抗旱能力［53］，在干旱和鹽脅迫下會(huì)大量表達(dá)，以增強(qiáng)植物的耐鹽性［54］。此外，有研究發(fā)現(xiàn)，在非生物脅迫下，植物中NAC、bZIP、MYB和WRKY家族的轉(zhuǎn)錄因子會(huì)大量表達(dá)［55-56］。本研究相應(yīng)基因及轉(zhuǎn)錄因子的表達(dá)，可能反映了甘青蒿對青藏高原極端環(huán)境的適應(yīng)性。

本研究利用ONT測序技術(shù)對甘青蒿全長轉(zhuǎn)錄組進(jìn)行測序和組裝，鑒定新基因及新轉(zhuǎn)錄本，并對其進(jìn)行功能注釋，篩選出的9個(gè)與萜類化合物生物合成相關(guān)候選基因，有助于為甘青蒿中新基因的發(fā)現(xiàn)和功能基因分析提供可靠數(shù)據(jù)。分析甘青蒿全長轉(zhuǎn)錄組的轉(zhuǎn)錄因子，能夠?yàn)楦是噍锕δ芑蛲诰蚝蜕镔Y源利用提供參考。