百山祖冷杉葉片轉(zhuǎn)錄組分析

王光炯，柳新紅，許大明，吳義松，李因剛，于明堅(jiān)，張麗芳，*

（1.浙江省林業(yè)科學(xué)研究院，浙江杭州 310023；2.浙江農(nóng)林大學(xué)林業(yè)與生物技術(shù)學(xué)院，浙江杭州 311300；3.浙江鳳陽山百山祖國家級自然保護(hù)區(qū)管理局百山祖管理處，浙江慶元 323800；4.浙江大學(xué)生命與科學(xué)學(xué)院，浙江杭州 310058）

【研究意義】百山祖冷杉（Abies beshanzuensisM.H.Wu）為松科（Pinaceae）冷杉屬（Abies）常綠喬木[1]。我國冷杉屬植物種類豐富，約有22個種3個變種，其中如百山祖冷杉、元寶山冷杉（Abies yuanbaoshanensisY.J.Lu et L.K.Fu）、梵凈山冷杉（Abies fanjingshanensisW.L.Huang，Y.L.Tu&S.Z.Fang）和資源冷杉（Abies beshanzuensisvar.ziyuanensis（L.K.Fu et S.L.Mo）L.K.Fu et Nan Li）等分布狹窄且數(shù)量極少，為國家一級重點(diǎn)保護(hù)植物[2]。百山祖冷杉于1963年由吳鳴翔先生首先發(fā)現(xiàn)，1976年正式命名發(fā)表，打破了許多學(xué)者認(rèn)為我國東南大陸地區(qū)不可能有冷杉分布的論斷[3-4]。它自然分布于浙江省麗水市慶元縣百山祖國家級自然保護(hù)區(qū)海拔1 740～1 750 m處，目前原生地僅存3株野生成熟植株。現(xiàn)階段由于百山祖冷杉現(xiàn)存數(shù)量稀少，難開花結(jié)果，自然繁殖十分困難等原因，急需加強(qiáng)科學(xué)監(jiān)測、研究和有效保護(hù)[8-10]。在百山祖冷杉保護(hù)和研究中對瀕危機(jī)制有初步的了解，除了開“花”失敗、自身花粉無效傳播和無法產(chǎn)生有活力種子等原因以外，幼苗和幼樹由于環(huán)境的變化（氣候變暖和環(huán)境污染）等原因不斷地死亡。然而，到目前為止，對于百山祖冷杉的生物學(xué)-生態(tài)學(xué)特性認(rèn)識嚴(yán)重不足，適應(yīng)脅迫環(huán)境的分子調(diào)控機(jī)制也尚不清楚?！厩叭搜芯窟M(jìn)展】伴隨后基因組時(shí)代來臨，轉(zhuǎn)錄組學(xué)連接著基因組與蛋白質(zhì)組，是率先發(fā)展并應(yīng)用最廣、最活躍的一門技術(shù)。近年來，國內(nèi)外的學(xué)者對許多瀕危植物展開了轉(zhuǎn)錄組測序研究，例如，在建立多葉斑葉蘭（Goodyera foliosa（Lindl.）Benth.）的繁殖體系的研究中同時(shí)結(jié)合高通量轉(zhuǎn)錄組測序技術(shù)，成功建立植株再生體系并得到全面完整的轉(zhuǎn)錄組信息，為其擴(kuò)繁、基因功能、遺傳發(fā)育以及調(diào)控機(jī)制研究奠定基礎(chǔ)[5-7]?！颈狙芯壳腥朦c(diǎn)】目前Illumina、Pacbio、SLAF-seq等技術(shù)被廣泛應(yīng)用于功能基因研究、分子標(biāo)記的開發(fā)、代謝途徑研究以及應(yīng)對各種環(huán)境脅迫反應(yīng)機(jī)制研究等方面[11-14]。在冷杉屬植物朝鮮冷杉（Abies koreana）的研究中，通過二代轉(zhuǎn)錄組分析技術(shù)分析熱脅迫下表達(dá)的基因，提供了對遭受熱脅迫的朝鮮冷杉的第一個綜合表征，為識別耐高溫脅迫的品系提供了重要的參考指標(biāo)[15]?！緮M解決的關(guān)鍵問題】本文基于高通量二代轉(zhuǎn)錄組測序平臺對百山祖冷杉的轉(zhuǎn)錄組進(jìn)行測序，獲得大量的基因數(shù)據(jù)，并對組裝的轉(zhuǎn)錄本進(jìn)行分類統(tǒng)計(jì)、功能注釋，通過分析轉(zhuǎn)錄組數(shù)據(jù)研究百山祖冷杉代謝途徑關(guān)鍵基因，為后續(xù)的環(huán)境脅迫調(diào)控機(jī)理等方面提供研究基礎(chǔ)。另外，通過轉(zhuǎn)錄組測序和組裝，開發(fā)微衛(wèi)星分子標(biāo)記，研究百山祖冷杉母樹和子代個體的遺傳多樣性和親緣關(guān)系。

1 材料與方法

1.1 試驗(yàn)材料

于2019年6月在浙江省慶元百山祖保護(hù)區(qū)內(nèi)（27°45′N，119°11′E）采集百山祖冷杉測序材料，選取當(dāng)年萌發(fā)枝條上部幼嫩葉片，3個生物學(xué)重復(fù)，采集后放置于干冰內(nèi)短期保存，之后保存于-80 ℃超低溫冰箱，用植物RNA提取試劑盒提取RNA。

1.2 試驗(yàn)方法

1.2.1 RNA提取、文庫構(gòu)建以及Illumina測序 從百山祖冷杉葉片組織中提取總RNA，并且檢測RNA的純度和濃度[16-17]。構(gòu)建基于Illumina Hiseq高通量測序的cDNA文庫。文庫構(gòu)建后，本研究為無參百山祖冷杉轉(zhuǎn)錄組分析，先拼接由測序所得到的序列，以拼接好的轉(zhuǎn)錄本為參考序列，利用GO、KO、KOG、NR、等7大數(shù)據(jù)庫與轉(zhuǎn)錄組測序結(jié)果比對并進(jìn)行功能注釋，再分析其GO富集和KEEG富集[18]。

1.2.2 基因表達(dá)量分析 利用BioKanga‘maploci’將讀長進(jìn)行FPKM（Reads per kilobase per million reads）計(jì)數(shù)標(biāo)準(zhǔn)化[19]。FPKM能消除基因長度和測序量差異對計(jì)算基因表達(dá)的影響。利用Bowtie軟件將Reads與Unigene庫進(jìn)行比對并結(jié)合RSEM軟件對表達(dá)量水平進(jìn)行評估，通過FPKM值來表示對應(yīng)Unigene的表達(dá)豐度。

1.2.3 基因功能注釋 通過BLAST軟件將所得到的Unigene序列與NR、Swiss-Prot、GO、COG、KOG、egg‐NOG4.5、KEGG 7大公共數(shù)據(jù)庫進(jìn)行比對，得到所需的Unigene注釋信息[20]。

1.2.4 轉(zhuǎn)錄因子分析 將測序所得到的基因序列與轉(zhuǎn)錄因子數(shù)據(jù)庫在BLAST軟件中進(jìn)行比對[21]。

1.2.5 SSR位點(diǎn)篩選和引物設(shè)計(jì) 通過MISA（http:∕∕pgrc.ipkgatersleben.de∕misa∕）軟件對篩選得到的1 kb以上的Unigene進(jìn)行SSR分析，根據(jù)堿基的重復(fù)次數(shù)鑒定SSR類型[22]。通過Primer 3引物設(shè)計(jì)軟件對找到的SSR位點(diǎn)序列進(jìn)行引物設(shè)計(jì)[23]。通常引物設(shè)計(jì)需要注意的有以下幾點(diǎn):（1）擴(kuò)增的產(chǎn)物大小在100～300 bp，所設(shè)計(jì)引物的長度盡量在18～25 bp，引物序列中的GC含量控制在40%～65%；（2）退火溫度56～64 ℃，并且設(shè)計(jì)的上下游引物退火溫度的差值小于5 ℃。設(shè)計(jì)引物時(shí)，引物不能出現(xiàn)SSR位點(diǎn)序列。（3）設(shè)計(jì)時(shí)應(yīng)避免出現(xiàn)發(fā)卡結(jié)構(gòu)、二聚體以及錯配等情況。隨機(jī)選取10對引物，委托上海捷瑞生物有限公司合成。

2 結(jié)果與分析

2.1 百山祖冷杉轉(zhuǎn)錄組測序結(jié)果以及數(shù)據(jù)組裝

基于二代測序技術(shù)Illumina Hiseq高通量測序平臺對百山祖冷杉的葉片進(jìn)行轉(zhuǎn)錄組測序，最終獲得大量原始數(shù)據(jù)（Raw data），原始測序序列去除一些低質(zhì)量以及含有測序引物、接頭等人工序列的reads后，共獲得7.02 Gb Clean Data，Q30堿基百分比在95.43%及以上。利用Trinity軟件對Clean Data進(jìn)行組裝，組裝得到73 103條轉(zhuǎn)錄本序列，總長度為9 303 348 bp，N50長度為1 930 bp，平均長度為1 272.64 bp。得到41 228條Unigene，總長為42 458 236 bp，Unigene的N50為1 690 bp，平均長度1 029.84 bp，組裝完整性較高（表1）。

表1 拼接長度分布Tab.1 Splicing length distribution

2.2 百山祖冷杉轉(zhuǎn)錄組基因功能注釋

通過BLAST（E-value＜1e-5）軟件將百山祖冷杉測序所得到的所有Unigene與Nr、eggNOG4.5、Swiss-Prot、KOG、COG、GO、KEGG數(shù)據(jù)庫進(jìn)行比對，得到有注釋信息的Unigene為27 738個，占所有Unigene的67.28%。其中注釋到Nr數(shù)據(jù)庫的數(shù)量最多為27 418，占所有Unigene的66.50%。注釋到COG數(shù)據(jù)庫的Unigene數(shù)量最少，僅有9 411條占所有Unigene的22.82%（表2）。

表2 Unigene注釋統(tǒng)計(jì)Tab.2 Unigene annotation statistics

Nr數(shù)據(jù)庫中E值分布的分析結(jié)果顯示，E值分布于0～1e-150的Unigene最多，高達(dá)6 711條，占注釋到Nr數(shù)據(jù)庫的Unigene的24.48%。同一性分布顯示，匹配的同一性為80%～90%的Unigene最多，達(dá)到4 908條，占17.90%。另外，對Nr注釋到基因同源序列的物種分布進(jìn)行了分析（圖1）。

2.3 Unigene的功能分類

將百山祖冷杉的全部Unigene與GO數(shù)據(jù)庫進(jìn)行比對，并將比對到的功能分類進(jìn)行統(tǒng)計(jì)。根據(jù)結(jié)果表明有16 020條Unigene在GO數(shù)據(jù)庫中注釋，匹配包括了3大類，分別為細(xì)胞組分、分子功能和生物學(xué)過程。注釋到細(xì)胞組分大類共15個亞類的基因有33 892個，注釋到分子功能大類共15個亞類的數(shù)量為18 407個，注釋到生物學(xué)過程大類共20個亞類的基因共有29 809個。在細(xì)胞組分的大類中，細(xì)胞（cell）和細(xì)胞組分（cell part）占比最高，分別為6 891和6 890個。在分子功能的大類中催化活性（catalytic activity）亞類注釋的數(shù)量最多為8 185個，其次綁定功能（binding）注釋數(shù)量為7 349。生物學(xué)過程大類中占比最高的亞類為代謝過程（metabolic process）注釋數(shù)量為8 080，其次細(xì)胞過程（cellular process）注釋數(shù)量為7 107（圖2）。

圖1 Unigene在Nr數(shù)據(jù)庫中的各項(xiàng)分布Fig.1 Distribution of Unigene in Nr Database

圖2 Unigene的GO功能分類Fig.2 GO function classification diagram of Unigene

將百山祖冷杉的所有Unigene與COG數(shù)據(jù)庫進(jìn)行比對，發(fā)現(xiàn)有9 411條Unigene得到注釋。分配至26類功能中，基因的數(shù)量在不同的功能分類存在明顯的差異。其中，僅一般功能預(yù)測類基因數(shù)量最多，達(dá)到1 261條，其次糖運(yùn)輸和新陳代謝功能類別，注釋的基因多達(dá)1 028條。此外，翻譯、核糖體結(jié)構(gòu)、生物合成功能類別（970條）與翻譯后修飾、蛋白質(zhì)折疊和分子伴侶功能類別（813條）注釋得到的基因也較多（圖3）。

圖3 COG功能分類統(tǒng)計(jì)Fig.3 COG function classification statistics

2.4 基因轉(zhuǎn)錄因子分析

在植物的生長發(fā)育過程中常會遇到干旱、低溫、洪澇等多種非生物因子的脅迫作用。在植物的正常生長發(fā)育中常常會遭受逆境的影響，影響嚴(yán)重時(shí)甚至對植物的傷害是不可逆的，最終導(dǎo)致植株的死亡[24]。植物中含有大量的轉(zhuǎn)錄因子，當(dāng)植物參與逆境反應(yīng)時(shí)相關(guān)轉(zhuǎn)錄因子調(diào)控脅迫應(yīng)答基因表達(dá)，使植物產(chǎn)生抗逆反應(yīng)，通過生理生化的變化來減少自身在逆境中受到的傷害。百山祖冷杉轉(zhuǎn)錄組測序共獲得1 872個轉(zhuǎn)錄因子，分布在208個轉(zhuǎn)錄因子家族中（圖4）。其中C2H2、AP2∕ERF-ERF、MYB-related、bHLH、MYB、C3H等轉(zhuǎn)錄因子家族所包含的轉(zhuǎn)錄因子較多，分別占4.70%、3.63%、3.25%、2.99%、2.93%、2.88%，這些轉(zhuǎn)錄因子家族主要對植物生長發(fā)育的調(diào)控起到重要作用。

圖4 轉(zhuǎn)錄因子家族分布Fig.4 Transcription factor family distribution

表3 部分轉(zhuǎn)錄家族及相應(yīng)Unigene數(shù)量Tab.3 Partial transcription factor family and corresponding Unigene numbers

圖5 Unigene的代謝通路分析Fig.5 Unigene’s metabolic pathway analysis

2.5 Unigene的代謝通路分析

將所得到的Unigene與KEGG數(shù)據(jù)庫進(jìn)行比對后進(jìn)行代謝通路的分析，結(jié)果顯示共有9 976條Unigene成功獲得注釋。根據(jù)功能將這些所注釋的Unigene分為6個大類，其中有21個亞類涉及126條代謝途徑。6大類分別為新陳代謝、遺傳信息處理、細(xì)胞過程、有機(jī)系統(tǒng)、環(huán)境信息處理、人類疾病。新陳代謝大類涉及95條代謝通路3 102個基因，遺傳信息處理大類涉及21條代謝途徑1 949個基因，細(xì)胞過程大類涉及3條代謝途徑471個基因，環(huán)境信息處理大類涉及3條代謝途徑272個基因，有機(jī)系統(tǒng)和人類疾病大類均涉及2條代謝途徑所包含的基因分別為162個和15個。其他代謝通路的詳細(xì)信息見圖5。

2.6 百山祖冷杉轉(zhuǎn)錄組SSR位點(diǎn)分析

將所得到的數(shù)據(jù)進(jìn)行篩選，選出長度大于1 kb的Unigene，通過MISA軟件進(jìn)行SSR分析，共得到1 794個SSR位點(diǎn)，相比其他物種的SSR位點(diǎn)數(shù)量，百山祖冷杉的位點(diǎn)數(shù)量顯著減少。利用Primer 3軟件對SSR位點(diǎn)兩側(cè)的序列進(jìn)行引物設(shè)計(jì)，共獲得100多對SSR引物，部分引物信息見表4。

表4 百山祖冷杉SSR部分引物信息Tab.4 Information of part pairs of development from Abies beshanzuensis

2.7 百山祖冷杉基因表達(dá)量的分析

百山祖冷杉的轉(zhuǎn)錄組測序中共得到23 540 099條reads，其中有20 271 952條reads比對到基因圖上，占總量的86.12%，有13 041 470條reads比對到基因圖唯一位置，占總量的64.33%，比對到多個位置的reads為7 230 482占35.67%。測序得到73 103條轉(zhuǎn)錄本序列，41 228條Unigene。41 228條Unigene的FPKM值的平均值為17.28，范圍在0～28 221.41，其中，6 786條Unigene的FPKM值在平均值以上。

3 結(jié)論與討論

生物多樣性喪失和物種滅絕是當(dāng)今世界面臨的重大生態(tài)危機(jī)之一，許多珍稀瀕危植物是自然生態(tài)系統(tǒng)中的關(guān)鍵種，在維持自然生態(tài)系統(tǒng)運(yùn)行中發(fā)揮著重要作用，一個物種的滅絕，有可能發(fā)生連鎖反應(yīng)，對自然生態(tài)系統(tǒng)的穩(wěn)定性造成一定的影響，導(dǎo)致災(zāi)難性后果[25-27]。百山祖冷杉從1976年發(fā)現(xiàn)以來，基礎(chǔ)研究取得了一些成果，但是在外界環(huán)境中，百山祖冷杉對環(huán)境因子的脅迫反應(yīng)機(jī)制一直不是很清楚。在植物的逆境脅迫中，經(jīng)常受到非生物逆境（如干旱、鹽、極端溫度及重金屬脅迫等）的侵害，在應(yīng)答干旱脅迫方面，DREB、WRKY、NAC、bHLH及bZIP等轉(zhuǎn)錄因子受干旱信號誘導(dǎo)，調(diào)節(jié)下游抗旱基因的表達(dá)，從而提高蔬菜作物抗旱能力。同時(shí)，水分運(yùn)輸相關(guān)功能基因（PIP、TIP）、E3連接酶SIZ1及脫水蛋白DHN也被報(bào)道受干旱誘導(dǎo)，并通過調(diào)節(jié)水勢、滲透勢及ROS積累抵御干旱脅迫。在抵御鹽脅迫方面，SOS途徑至關(guān)重要。Sl SOS2能夠通過調(diào)節(jié)Sl SOS1和Na+∕H+逆向轉(zhuǎn)運(yùn)蛋白Le NHX2∕4的表達(dá)維持離子平衡和調(diào)節(jié)植物器官中Na+的分配。蔬菜抗鹽研究中NAC、ERF、MYB等轉(zhuǎn)錄因子響應(yīng)鹽脅迫并激活抗逆相關(guān)基因表達(dá)，從而提高蔬菜作物抗鹽能力[28]。

本文利用RNA-seq技術(shù)[29]，對百山祖冷杉的轉(zhuǎn)錄組進(jìn)行測序、組裝，并與Nr、GO、KEGG等7大公共數(shù)據(jù)庫進(jìn)行比對。GO分析顯示，獲得注釋的基因共參與了43個通路，最富集的通路主要是在生物學(xué)過程中，包括免疫系統(tǒng)形成、發(fā)育過程和代謝等，而代謝過程獲得的注釋是最多的。而在COG數(shù)據(jù)庫比對分析結(jié)果中，糖運(yùn)輸和新陳代謝功能注釋的基因最多，植物在干旱脅迫條件下，體內(nèi)水分代謝與碳代謝會發(fā)生失衡現(xiàn)象，光合速率降低、蒸騰速率降低，導(dǎo)致生長降低，使植物陷入饑餓狀態(tài)[30]。糖作為重要的代謝物質(zhì)，生物合成不僅為植物的生長發(fā)育提供碳源和能源，而且分解產(chǎn)生的小分子單糖可以調(diào)控植物細(xì)胞的滲透壓以抵抗脅迫[31]。GO分析和COG數(shù)據(jù)庫分析的結(jié)果一致，進(jìn)一步說明百山祖冷杉在生長過程中代謝活動旺盛，而且一些次生代謝物的生物合成途徑，包括生物堿、萜類化合物和黃酮類化合物等在抵抗病原物的入侵方面起著重要的作用[32]。

轉(zhuǎn)錄因子在植物的生長和調(diào)控中起著重要的作用[33]。在植物對抗環(huán)境脅迫過程中，調(diào)控相關(guān)生理反應(yīng)基因表達(dá)的蛋白在防衛(wèi)反應(yīng)和逆境信號轉(zhuǎn)導(dǎo)中發(fā)揮作用，使植物適應(yīng)外界不良環(huán)境[34-36]。百山祖冷杉轉(zhuǎn)錄因子分析表明，C2H2、AP2∕ERF-ERF、MYB-related、bHLH、MYB、C3H等轉(zhuǎn)錄因子家族所包含的轉(zhuǎn)錄因子較多，這些轉(zhuǎn)錄因子家族在植物發(fā)育中起到關(guān)鍵調(diào)控作用。已有的研究表明，C2H2鋅指轉(zhuǎn)錄因子家族主要參與植物葉的發(fā)生，花器官的調(diào)控，側(cè)枝的形成及逆境脅迫等生命過程[37]，意味著C2H2鋅指轉(zhuǎn)錄因子家族在百山祖冷杉生長各階段、各組織中都參與表達(dá)。另外，有研究表明AP2∕ERF-ERF家族成員能夠參與調(diào)控植物花的形成和發(fā)育，種子的大小，種子的重量，種子中脂肪酸的合成和抗逆等。該家族基因在種子不同部位存在特異性表達(dá)，在種子發(fā)育的各個過程中起著不同的作用[38]。各轉(zhuǎn)錄因子分析研究為百山祖冷杉種子敗育的機(jī)理研究提供參考，同時(shí)為今后的幼苗繁育及生境改良工作提供分析參考。

通過轉(zhuǎn)錄組測序得到的數(shù)據(jù)開發(fā)SSR引物，應(yīng)用于百山祖冷杉遺傳多樣性等方面研究。前人研究中已有對百山祖冷杉遺傳多樣性的研究[39-40]，但大多為ISSR、RAPD等分子標(biāo)記[38]，遺傳多樣性研究中存在引物數(shù)量有限等原因，無法將百山祖冷杉的雜交實(shí)生苗或無性繁殖的個體與其親本進(jìn)行鑒別[41]。通過轉(zhuǎn)錄組的高通量測序?qū)?shù)據(jù)進(jìn)行篩選和SSR分析得到1 794個SSR位點(diǎn)，所得到的位點(diǎn)數(shù)量并不多，重復(fù)類型以單核苷酸重復(fù)、雙核苷酸重復(fù)和三核苷酸重復(fù)為主，其中單核苷酸重復(fù)類型最多。根據(jù)這些位點(diǎn)信息設(shè)計(jì)了大量的SSR引物，為今后百山祖冷杉母本與子代親緣關(guān)系鑒別以及進(jìn)一步的遺傳多樣性分析提供基礎(chǔ)信息。

本研究基于百山祖冷杉轉(zhuǎn)錄組測序所獲得的大量數(shù)據(jù)，通過后續(xù)分析系統(tǒng)了解了百山祖冷杉的基因功能注釋、代謝通路等，為今后對其逆境脅迫、遺傳背景、基因功能以及遺傳標(biāo)記等研究提供基礎(chǔ)數(shù)據(jù)，為百山祖冷杉的保護(hù)以及種群的擴(kuò)大繁殖提供理論基礎(chǔ)。但現(xiàn)階段對于百山祖冷杉的研究剛在起步階段，由于技術(shù)手段的限制，未對研究中所注釋的Unigene基因進(jìn)行定量PCR驗(yàn)證，進(jìn)一步確定轉(zhuǎn)錄組數(shù)據(jù)的準(zhǔn)確和可實(shí)用性，而這些存在的問題也是接下來百山祖冷杉研究工作的重點(diǎn)。