王建秋， 王曉麗， 曹子林

(西南林業(yè)大學(xué)， 云南 昆明 650224)

滇白前(SileneviscidulaFranch)又名黏萼蠅子草，為石竹科(Cayophyllaceae)蠅子草屬的一種多年生草本植物，主要生長在海拔1 200～3 200 m的四川、云南、貴州、西藏(東南)等地[1]。據(jù)研究，滇白前對鉛、鋅、鎘具有很強的耐受性，是一種新的Pb/Zn/Cd共超富集植物，為土壤重金屬復(fù)合污染地修復(fù)提供了新的種質(zhì)資源[2]。目前，有關(guān)滇白前的研究主要集中于藥用成分及抗癌作用上[3-4]。種子是植物生命周期下一代的載體[5]，對于這一具有重要藥用價值的重金屬超富集植物，其種子基因組結(jié)構(gòu)及其功能是否有其特殊性未見報道。轉(zhuǎn)錄組是特定的細(xì)胞或組織在某一發(fā)育階段或功能狀態(tài)下轉(zhuǎn)錄出的所有RNA總和[6]，功能基因研究的主要方法是基于高通量測序的轉(zhuǎn)錄組分析[7]。轉(zhuǎn)錄組測序分析不僅能迅速獲得基因轉(zhuǎn)錄組信息、功能注釋和樣品間的差異表達(dá)，而且能揭示特定生物學(xué)過程下，器官、細(xì)胞的分子機制[8]。本研究擬進(jìn)行滇白前種子的轉(zhuǎn)錄組測序，分析測序質(zhì)量并進(jìn)行基因功能注釋及轉(zhuǎn)錄因子預(yù)測，為深入挖掘滇白前的功能基因及分子機制等方面的研究奠定一定的基礎(chǔ)。

1 材料與方法

1.1 樣品采集

樣品于2020年1月1日采自云南省怒江州蘭坪縣金頂鎮(zhèn)鳳凰山金頂鉛鋅礦區(qū)成熟的滇白前種子，地處北緯26°21′～27°02′，東經(jīng)98°23′～99°28′，此礦區(qū)具有儲量大、品位高、儲存集中、埋藏淺等特點，是亞洲最大的鉛鋅礦床[2]。

1.2 RNA提取、文庫的構(gòu)建及測序

從種子中提取總RNA后，利用帶有Oligo(dT)的磁珠與ployA進(jìn)行A-T堿基配對，從總RNA中分離出mRNA；富集得到完整的RNA序列，通過加入fragmentation buffer,將mRNA隨機斷裂、磁珠篩選分離出300 bp左右的小片段；在逆轉(zhuǎn)錄酶作用下，加入六堿基隨機引物，以mRNA為模板反轉(zhuǎn)錄合成一鏈cDNA，隨后進(jìn)行二鏈合成，形成穩(wěn)定的雙鏈結(jié)構(gòu)；連接adaptor，用2%的瓊脂糖膠回收目的條帶大小合適的片段并通過PCR技術(shù)進(jìn)行擴增；最后，用Illumina Novaseq 6000對短序列片段進(jìn)行測序。

1.3 測序原始數(shù)據(jù)過濾及組裝

測序獲得原始數(shù)據(jù)后，使用SeqPrep軟件將原始數(shù)據(jù)中含有大量的接頭序列、低質(zhì)量序列末端以及adapter和N的比例大于10%的堿基去除，從而得到高質(zhì)量的質(zhì)控數(shù)據(jù)(clean data)。使用Trinity軟件[9]對所有clean data進(jìn)行從頭組裝，使用TransRate[10]對從頭拼接結(jié)果常見的錯誤(包括嵌合體、結(jié)構(gòu)錯誤、組裝不完整、堿基錯誤等)進(jìn)行評估，給出每條contig的質(zhì)量得分，并可將這些contigs評分整合獲得整個組裝結(jié)果的綜合分?jǐn)?shù)，對轉(zhuǎn)錄組獲得的序列進(jìn)行過濾和優(yōu)化。使用CD-HIT軟件[11]通過序列比對聚類(Cluster)的方法去除冗余、相似的序列，最后獲得非冗余(non-redundant)的序列。BUSCO[12]利用單拷貝直系同源基因，評估基因組或轉(zhuǎn)錄組的組裝完整性。基因組組裝是與BUSCO的一致序列進(jìn)行tBLASTn比對，然后再用Augustus對基因結(jié)構(gòu)進(jìn)行預(yù)測，而轉(zhuǎn)錄組組裝則先尋找轉(zhuǎn)錄本的ORF編碼框，然后再進(jìn)行HMMER 3比對。

1.4 基因功能注釋

用NCBI將Unigene序列比對到蛋白數(shù)據(jù)庫Nr、SwissProt和COG(evalue<0.000 01)，用KOBAST軟件將Unigene序列進(jìn)行KEGG注釋，用BLAST 2 GO軟件將Unigene序列進(jìn)行GO注釋[13]。得到與給定Unigene具有最高序列相似性的蛋白，從而獲得該Unigene的蛋白功能注釋信息，并將注釋結(jié)果進(jìn)行分類統(tǒng)計。將預(yù)測的蛋白序列同相應(yīng)的TF數(shù)據(jù)庫(plant TFdb/animal TFdb)進(jìn)行hmmscan比對，得到相應(yīng)的轉(zhuǎn)錄因子家族。

2 結(jié)果與分析

2.1 滇白前種子轉(zhuǎn)錄組組裝質(zhì)量分析

對滇白前種子進(jìn)行轉(zhuǎn)錄組測序后，通過嚴(yán)格質(zhì)量控制和數(shù)據(jù)過濾，共獲得50 423 900個高質(zhì)量干凈序列。利用Trinity軟件，對這些高質(zhì)量數(shù)據(jù)進(jìn)行denovo組裝，共獲得有效序列片段32 916 084條，N 50為1 192 bp，平均長度為753.87 bp，GC核苷酸含量(%)為43.45%的滇白前Unigene 43 663條(不含N的組裝片段)，長度分布于201～15 324 bp之間(表1)。通過用長度分布柱狀圖對組裝后的基因做統(tǒng)計(圖1)，長度為200～500 bp的有23 651條Unigene(54%)，長度為501～1 000 bp的有9 625條Unigene(22%)，長度為1 001～2 000 bp的有4 606條Unigene(11%)，長度為2 001～3 000 bp的有1 407條Unigene(3%)，長度超過3 000 bp的有797條Unigene(1%)。由此可得，隨著基因長度的增加，基因數(shù)量呈下降趨勢。

表1 Unigene基本信息Table 1 The basic information Unigenes

2.2 滇白前基因的總體注釋情況

將比對到植物中的基因與公共數(shù)據(jù)庫進(jìn)行功能注釋，分別有26 688、20 818、25 657、19 942、14 192條Unigenes比對到Nr、Swissprot、COG、GO、KEGG數(shù)據(jù)庫，至少有29 177條Unigenes比對到一個數(shù)據(jù)庫，并注釋到功能，有14 486條Unigenes沒有比對到數(shù)據(jù)庫，說明有許多功能不清的基因。其中，Nr數(shù)據(jù)庫注釋到的基因最多，共有26 688條，占比為61.12%；KEGG數(shù)據(jù)庫注釋到的基因最少，共有11 527條，占比為26.4%(表2)。

表2 Unigene注釋統(tǒng)計Table 2 Annotation statistics of Unigenes

2.3 滇白前基因的NR功能注釋

共有26 688條Unigenes比對到NR數(shù)據(jù)庫，與其他物種比對后的結(jié)果如圖2。注釋序列物種中注釋到基因數(shù)量最多的是甜菜(Betavulgaris)，共有6 828條，占比為25.58%；其余依次為藜麥(Chenopodiumquinoa)、栓皮櫟(Quercussuber)、菠菜(Spinaciaoleracea)，分別有5 997、4 738、4 272條Unigenes，占比分別為22.47%、17.75%、16.01%；其他289個物種共有4 853條Unigenes，占18.18%。

2.4 COG功能注釋

將滇白前基因比對到COG數(shù)據(jù)庫發(fā)現(xiàn)，有26 500條注釋到了COG數(shù)據(jù)庫中，占總Unigenes的61.12%，這些Unigenes被分為23類(圖3)。有13 917個功能注釋信息未知，未確定其準(zhǔn)確的生物學(xué)功能，占所有功能注釋信息的52.51%；翻譯后修飾與轉(zhuǎn)運注釋到1 807條Unigenes，所占比例為6.82%；翻譯、核糖體結(jié)構(gòu)與生物合成注釋到1 320條Unigenes，所占比例為4.98%；轉(zhuǎn)錄注釋到1 238條Unigenes，所占比例為4.67%；胞內(nèi)運輸、分泌和囊泡運輸注釋到1 224條Unigenes，所占比例為4.62%；信號傳導(dǎo)機制注釋到1 073條Unigenes，占4.05%；最少的是細(xì)胞能動性和核酸結(jié)構(gòu)的Unigenes，分別僅有7條(0.03%)和3條(0.01%)。這些結(jié)果表明，滇白前在蛋白質(zhì)翻譯后修飾與轉(zhuǎn)運、翻譯、核糖體結(jié)構(gòu)與生物合成、轉(zhuǎn)錄、胞內(nèi)運輸、分泌和囊泡運輸、信號轉(zhuǎn)導(dǎo)機制等基因表達(dá)豐度較高。

2.5 GO功能注釋

一共有19 942條Unigenes被注釋到了GO數(shù)據(jù)庫中，得到79 481個功能注釋，全面性地描述了不同生物基因的生物學(xué)特征，共分為三大功能，分別是生物過程、細(xì)胞組分和分子功能，依次得到24 770、30 335、24 376個功能注釋，占比分別為31.16%、38.17%、30.67%。表3中展示的是注釋上的GO二級分類條目情況。3個本體細(xì)分為52個功能亞類，生物學(xué)過程類占其中23個功能亞類，細(xì)胞進(jìn)程(8 260個注釋)和代謝過程(7 051個注釋)占比較大，其次是生物調(diào)控(2 466個注釋)；分子功能類包括16個功能亞類，結(jié)合(10 327個注釋)所占比例最多，其次是催化活性(10 240個注釋)；細(xì)胞組分類包括13個亞類，“細(xì)胞”所占比例最多(9 017個注釋)，其次是“膜”(6 205個注釋)和“細(xì)胞器”(5 066個注釋)。只有少數(shù)轉(zhuǎn)錄本被注釋參與細(xì)胞殺死、行為、蛋白質(zhì)標(biāo)簽、翻譯調(diào)節(jié)器活動、營養(yǎng)物儲存活動及細(xì)胞外基質(zhì)。結(jié)果表明，以細(xì)胞過程、代謝過程、結(jié)合和催化活性相關(guān)的基因較多。對GO功能注釋進(jìn)一步分析，梳理得到許多抗逆性相關(guān)注釋，分別是環(huán)境應(yīng)激響應(yīng)、蛋白激酶、激素代謝、滲透調(diào)節(jié)、氧化還原反應(yīng)、受體和離子轉(zhuǎn)運相關(guān)的注釋(表4)。

表4 滇白前轉(zhuǎn)錄組中與逆境脅迫相關(guān)的基因Table 4 Adversity stress-related genes in the transcriptome of S. viscidula

表3 滇白前GO注釋結(jié)果Table 3 GO Annotation of S. viscidula

2.6 KEGG功能注釋

有11 527個Unigenes比對到KEGG數(shù)據(jù)庫的參考代謝通路中，可將滇白前Unigenes歸為6個類別，代謝相關(guān)的通路Unigenes有5 155條，占比最多為44.72%；其次是遺傳信息處理相關(guān)的通路Unigenes有3 582條，占比為31.07%；環(huán)境信息處理相關(guān)的通路Unigenes有491條，占比為4.26%；細(xì)胞過程相關(guān)的通路Unigenes有968條，占比為8.40%；生物體系統(tǒng)相關(guān)的通路Unigenes有794條，占比為6.89%；而人類疾病相關(guān)的通路Unigenes有537條，占比為4.66%。代謝相關(guān)的通路細(xì)致分成11個亞類，其中排名前三的是碳水化合物代謝、氨基酸代謝、能量代謝，占比分別為11.46%、7.37%、6.03%，遺傳信息加工和生物系統(tǒng)各分為4個亞類，環(huán)境信息處理和細(xì)胞過程各分為兩類，人類疾病相關(guān)的通路分為3個亞類，其中“翻譯”、“折疊、分類和降解”、“運輸和代謝”以及“環(huán)境適應(yīng)性”分別占15.47%、9.33%、6.11%、5.66%，其余11個亞類所占比例均小于5%(表5)。

表5 滇白前Unigenes的KEGG注釋結(jié)果及分類Table 5 KEGG annotation and classification of Unigenes in S. viscidula

2.7 轉(zhuǎn)錄組Unigene的轉(zhuǎn)錄因子預(yù)測

滇白前轉(zhuǎn)錄組共預(yù)測到527個轉(zhuǎn)錄因子，可分為32個轉(zhuǎn)錄因子大家族。其中，MYB(65個，12.33%)屬于最大家族，其次是ERF(54個，10.25%)，接著為bHLH(43個，8.16%)、C 2 C 2(38個，7.21%)、bZIP(35個，6.64%)、B 3(35個，6.64%)、C 3 H(29個，5.5%)、GRAS(24個，4.55%)、WRKY(24個，4.55%)和LBD(24個，4.55%)，其他13個家族共有基因132個，占25.05%(圖4)。

2.8 SSR分析結(jié)果統(tǒng)計

采用MISA基于拼接所得轉(zhuǎn)錄本序列信息進(jìn)行SSR分析(表6)，篩選得到43 663個Unigenes，識別的SSRs總數(shù)為4 730個，包含SSR的序列數(shù)目為3 995個，共有6種SSR重復(fù)類型，其中單堿基重復(fù)SSR 2 554條，占54%；二堿基重復(fù)SSR 459條，占9.7%；三堿基重復(fù)SSR 1 538條，占32.53%；四堿基重復(fù)(67條)、五堿基重復(fù)(22條)和六堿基重復(fù)(90條)共僅占3.78%。單堿基重復(fù)SSR中，發(fā)生頻率最高的是A；二堿基重復(fù)SSR中，發(fā)生頻率最高的是AG，最低的是CG；三堿基重復(fù)SSR中，發(fā)生頻率最高的是ATC和AAG，頻率最低的是ACG。

表6 滇白前SSR分析Table 6 Statistics of SSR analysis of S. viscidula

3 結(jié)論與討論

高通量測序技術(shù)深受研究者歡迎，目前該技術(shù)廣泛應(yīng)用到植物特殊功能基因的挖掘與鑒定[14]。本研究利用Illumina Novaseq 6 000對滇白前種子進(jìn)行轉(zhuǎn)錄組測序，產(chǎn)生clean reads的Q 20值為97.79%，Q 30值為93.25%以上。經(jīng)de nove拼接組裝及去冗余處理后，共得到43 663條Unigene，最大長度為15 324 bp，最小長度為201 bp，N 50為1 192 bp，平均長度為753.87 bp。組裝結(jié)果明顯高于白刺花(Sophoraviciifolia)種子(N 50為537 bp，平均長度為462 bp)[15]，而與油菜(Brassicanapus)(N 50為898 bp，平均長度659 bp)較為接近[16]。

COG、GO和KEGG注釋是基因注釋的重要途徑，GO的功能分類對了解基因功能有重要作用，而KEGG數(shù)據(jù)庫中的參考通路不僅可以推測基因功能，而且可以研究基因在不同代謝通路中的位置及作用[17]。通過COG功能分類可知，參與蛋白質(zhì)翻譯后修飾與轉(zhuǎn)運、翻譯、核糖體結(jié)構(gòu)與生物合成的Unigene最多，分別是1 807條、1 320條，發(fā)現(xiàn)有13 917條未知功能基因，這有可能與序列片段過短、注釋信息缺乏和獲得新基因等有關(guān)系[18]。與COG數(shù)據(jù)庫的比對，為了解基因功能及進(jìn)化研究奠定了一定的基礎(chǔ)[19]。GO功能共分為3個大類和52個亞類，其中參與細(xì)胞過程和代謝過程的Unigene數(shù)目最多，分別是8 260條和7 051條，并深入挖掘到從屬于環(huán)境應(yīng)激響應(yīng)、蛋白激酶、植物內(nèi)源激素、滲透調(diào)節(jié)、氧化還原、受體和離子轉(zhuǎn)運的與植物抗逆相關(guān)的Unigenes。通過KEGG數(shù)據(jù)庫和通路分析發(fā)現(xiàn)，涉及代謝途徑的Unigene(5 155條)最多，通路數(shù)量最多的也是涉及代謝途徑的，共101條。其中，嘌呤代謝(298)和嘧啶代謝(231)是核苷酸代謝的核心，被視為植物界內(nèi)的管家功能，可以通過激活A(yù)BA代謝途徑從而發(fā)揮對非生物逆境脅迫條件的應(yīng)激保護作用[20]。植物激素信號轉(zhuǎn)導(dǎo)代謝途徑在KEGG中注釋到196個Unigenes，植物在非生物逆境脅迫條件下，積累ABA，通過激素激活應(yīng)激相關(guān)的基因參與基因表達(dá)和信號轉(zhuǎn)導(dǎo)，從而獲得承受惡劣環(huán)境的耐受能力[21]。

轉(zhuǎn)錄因子是讀取并解釋DNA中遺傳“藍(lán)圖”的蛋白質(zhì)組之一，其與DNA結(jié)合，可以幫助啟動一個增加或減少基因轉(zhuǎn)錄的程序，因此轉(zhuǎn)錄因子對于許多細(xì)胞過程是至關(guān)重要的。轉(zhuǎn)錄因子在脅迫刺激下不斷合成，將信號傳遞和放大，并通過與啟動子特定區(qū)域相結(jié)合調(diào)控下游響應(yīng)逆境相關(guān)基因的表達(dá)，從而引起植物生理生化的改變，在植物多種非生物逆境脅迫中起著關(guān)鍵作用[22]。根據(jù)轉(zhuǎn)錄因子預(yù)測的結(jié)果，滇白前中MYB_superfamily轉(zhuǎn)錄因子最多，其次是AP 2/ERF，接著是bHLH、C 2 C 2、bZIP、B 3、WRKY、NAC、C 2 H 2等轉(zhuǎn)錄因子。據(jù)研究，植物中的MYB、bZIP、WRKY、NAC、AP 2/ERF、bHLH、C 2 H 2等各類轉(zhuǎn)錄因子在響應(yīng)高鹽、干旱、寒冷等逆境脅迫中均有重要作用[23-28]。對這些轉(zhuǎn)錄因子進(jìn)行統(tǒng)計,大約有261個(49.5%)可能與其抗逆適應(yīng)相關(guān)。這也許是滇白前能在重金屬濃度很高的土壤中正常生長，也是其能分布在3 200 m高海拔地方的原因[1-2]。簡單序列重復(fù)SSR標(biāo)記是目前比較理想的分子標(biāo)記，滇白前轉(zhuǎn)錄組測定共篩選了3 995個SSRs序列，分為6種SSR重復(fù)類型，占比較多的是單堿基重復(fù)和三堿基重復(fù)。這為后續(xù)開展滇白前生長和抗逆分子機制研究提供了有力依據(jù)。

云南省素有“有色金屬王國”之美譽，土壤重金屬污染已成為有色金屬開采行業(yè)面臨的嚴(yán)重問題，而且污染土壤通常表現(xiàn)為多種重金屬復(fù)合的特點[2]。滇白前是鉛、鋅、鎘共超富集植物，是一種修復(fù)土壤重金屬復(fù)合污染的寶貴的物種資源[2]。本研究以滇白前種子為研究材料，利用高通量測序技術(shù)構(gòu)建了轉(zhuǎn)錄組數(shù)據(jù)庫，得到了許多的轉(zhuǎn)錄本序列信息，并進(jìn)行了功能注釋和分類、轉(zhuǎn)錄因子和SSR等分析，揭示了滇白前種子轉(zhuǎn)錄組的整體特征。這為研究滇白前的繁殖機制、環(huán)境適應(yīng)性機制等提供了依據(jù)。