尹明華 譚 鑫 鄭亞嬌 徐迎昕 邱夢琴曾清華 蔡 紅 陳榮華

(1 上饒師范學(xué)院生命科學(xué)學(xué)院，江西 上饒 334001；2 上饒市紅日農(nóng)業(yè)開發(fā)有限公司，江西 上饒 334700)

馬鈴薯(SolanumtuberosumL.)是世界各地廣泛種植的糧食作物和經(jīng)濟(jì)作物，擴(kuò)大馬鈴薯生產(chǎn)對保障全球糧食安全至關(guān)重要[1]。繼稻米、小麥、玉米之后，馬鈴薯已成為我國第四大主糧作物[2]。江西懷玉山平均海拔1 000 m，屬于高海拔地區(qū)，溫差較大，馬鈴薯的生長周期延長，因此該地區(qū)通常選擇一些早熟的馬鈴薯品種，如懷玉山高山馬鈴薯和懷玉山本地農(nóng)家薯[3]。懷玉山高山馬鈴薯，又稱麻籽洋芋，是2013年農(nóng)業(yè)部核準(zhǔn)的國家地理標(biāo)志農(nóng)產(chǎn)品，營養(yǎng)豐富，且不含還原糖，富含膳食纖維，有和中養(yǎng)胃、健脾利濕、降糖降脂、美容養(yǎng)顏、防治胃癌等功效。此外，研究發(fā)現(xiàn)懷玉山高山馬鈴薯雖對土壤條件要求較低，但在平原地區(qū)種植易出現(xiàn)變種現(xiàn)象[4]，而懷玉山本土農(nóng)家薯在高海拔地區(qū)和平原地區(qū)均能種植，且不會(huì)出現(xiàn)變種現(xiàn)象。目前關(guān)于懷玉山高山馬鈴薯與懷玉山本土農(nóng)家薯適應(yīng)高海拔生境的內(nèi)在機(jī)理的研究尚未見報(bào)道。

轉(zhuǎn)錄組(transcriptome)是植物在某種狀態(tài)或環(huán)境下所有基因轉(zhuǎn)錄產(chǎn)物的總和，受到內(nèi)外源因子的影響和調(diào)控，具有顯著的時(shí)間性和空間性[5]。轉(zhuǎn)錄組測序具有成本低、通量高、靈敏度高、操作簡單、應(yīng)用領(lǐng)域廣等優(yōu)點(diǎn)，且對于未知基因組序列的物種或已知其基因組序列信息的物種均可得到低轉(zhuǎn)錄水平的表達(dá)基因[6]。通過轉(zhuǎn)錄組數(shù)據(jù)的分析，能夠發(fā)掘新的功能基因[7]，篩選差異基因，對懷玉山高山馬鈴薯的保種育種具有重要意義[8]。

本研究利用Illumina Hiseq 2500對懷玉山高山馬鈴薯和懷玉山本土農(nóng)家薯進(jìn)行轉(zhuǎn)錄組測序分析，通過差異基因表達(dá)及KEGG通路分析，對懷玉山高山馬鈴薯高海拔生境下種植的特異性進(jìn)行研究，以期解析懷玉山高山馬鈴薯與懷玉山本土農(nóng)家薯適應(yīng)高海拔生境的內(nèi)在機(jī)理，為其適應(yīng)高海拔生境相關(guān)基因的挖掘與應(yīng)用提供參考。

1 材料與方法

1.1 試驗(yàn)材料

以懷玉山高山馬鈴薯(編號(hào)HAP)和懷玉山本土農(nóng)家薯(編號(hào)LAP)采后50 d且萌發(fā)出芽的塊莖為試驗(yàn)材料，由上饒市紅日農(nóng)業(yè)開發(fā)有限公司提供。

1.2 試驗(yàn)方法

1.2.1 RNA提取 準(zhǔn)備65℃預(yù)熱的含15% β-巰基乙醇的細(xì)胞裂解液溶液(CLB)。取馬鈴薯塊莖萌發(fā)芽于液氮中研磨，取1 g放入裝有3 mL CLB溶液(含15%β-巰基乙醇)的5 mL離心管中，立即振蕩60 s，13 000 r·min-1離心10 min。取上清液，加入等體積氯仿/異戊醇溶液，振蕩60 s，4 500 r·min-14℃離心 10 min。再取上清液加入等體積酚/氯仿/異戊醇溶液，振蕩60 s，13 000 r·min-1離心2 min，重復(fù)以上步驟一次，收集上清液，加入2.5倍體積的預(yù)冷無水乙醇和0.1倍體積5 mol·L-1乙酸鉀4℃過夜沉淀，然后 12 000 r·min-1離心20 min，棄去上清液，吹除乙醇，并用60 μL RNase-Free H2O復(fù)溶，即得到RNA。

1.2.2 Hiseq 2500測序 利用Bioanalyzer 2100(Agilent，Germany)控制總RNA的質(zhì)量。首先，5 μg總RNA通過Dynabeads Oligo(dT)(Life Technologies，USA)對mRNA進(jìn)行捕獲，并進(jìn)行片段化處理，然后用SuperScriptIII cDNA Synthesis Kit(Life Technologies，USA)進(jìn)行cDNA的反轉(zhuǎn)錄，再進(jìn)行cDNA末端修復(fù)、3′平末端加A、Adapters連接反應(yīng)最終連接至illumina的測序接頭并進(jìn)行PCR擴(kuò)增反應(yīng)。根據(jù)HiSeq SBS Kit and Cluster Kit v4(Illumina，Sandiego)說明書制備Total RNA測序文庫。建立的文庫進(jìn)一步用Qubit 2.0(Life Technologies，USA)和Bioanalyzer 2100(Agilent，Germany)進(jìn)行質(zhì)量控制。最后，按照 Hiseq 2500(Illumina，Sandiego)的操作說明對形成的cDNA文庫進(jìn)行2*125 bp的高通量測序(由蘇州帕諾米克生物科技有限公司完成測序)，設(shè)3次生物學(xué)重復(fù)。利用CASAVA 1.8軟件進(jìn)行圖像分析，堿基轉(zhuǎn)換，并根據(jù)不同的INDEX進(jìn)行多通路篩選。

1.3 數(shù)據(jù)分析

鑒于illlumina Hiseq 2500錯(cuò)誤率對結(jié)果的影響，對原始數(shù)據(jù)進(jìn)行質(zhì)量控制(quality control，QC)，包括去除低質(zhì)量序列，去除接頭，以獲得高質(zhì)量的序列進(jìn)行后續(xù)工作。質(zhì)控方法為ngsqctoolkit V2.33，參數(shù)默認(rèn)[cutoff read length for HQ (high-quality)=70%, cutoff quality score=20]。參考基因組為馬鈴薯數(shù)據(jù)：PGSC Version 4.03 Pseudomolecule Sequence，來源為PGSC(http://solanaceae.plantbiology.msu.edu)。利用軟件Tophat2分別對每個(gè)樣本進(jìn)行reads mapping，參數(shù)為-i 10-Ⅰ 15000。使用HTseq計(jì)算轉(zhuǎn)錄本的表達(dá)豐度和mapping到轉(zhuǎn)錄本的reads數(shù)，并用bioconductor edgeR包(www.bioconductor.org)進(jìn)行差異分析，根據(jù)edgeR中的exactTest(exprSet)函數(shù)進(jìn)行差異轉(zhuǎn)錄本的分析，并進(jìn)一步根據(jù)fold change和P值(FDR)方法找到差異轉(zhuǎn)錄本(與fold change閾值和FDR有關(guān)，使用4倍上調(diào)/下調(diào)和FDR<0.05)。對獲得的差異基因進(jìn)行樣本和代謝物的雙向聚類，即用層次聚類方法進(jìn)行雙向聚類。聚類軟件為R語言包(www.r-project.org)Pheatmap。應(yīng)用Bioconductor GOseq包的超幾何檢驗(yàn)(phyper)，找出與基因組的所有表達(dá)基因背景相比，在差異表達(dá)基因中顯著性富集的GO term和KEGG代謝通路，并進(jìn)行假陽性檢驗(yàn)(BH方法)。對P值<0.05的GO term/Pathway定義為在差異表達(dá)基因中顯著富集的GO term/Pathway。其中GO和Pathway的注釋信息來源于uniprot數(shù)據(jù)庫。

2 結(jié)果與分析

2.1 測序數(shù)據(jù)預(yù)處理

原始數(shù)據(jù)進(jìn)行質(zhì)量預(yù)處理后，HAP組和LAP組以及這兩組所有樣本的有效比分別為93.6%、93.9%和93.7%(表1)。上述表明，此次轉(zhuǎn)錄組測序數(shù)據(jù)質(zhì)量較好，可以滿足轉(zhuǎn)錄組分析的基本要求。

2.2 reads mapping統(tǒng)計(jì)

由表2可知，HAP組總mapping reads數(shù)為17 496 592，總mapping reads數(shù)占所有reads數(shù)的58.5%，唯一mapping reads數(shù)為16 354 632，唯一mapping reads數(shù)占總mapping reads的比例為93.5%；LAP組總mapping reads數(shù)為17 014 468，總mapping reads數(shù)占所有reads數(shù)的54.1%，唯一mapping reads數(shù)為16 235 240，唯一mapping reads數(shù)占總mapping reads的95.4%。

2.3 Transcript差異研究

2.3.1 豐度分析 由表3可知，HAP組和LAP組表達(dá)量最高的前10個(gè)基因?yàn)镻GSC0003DMT400045665、PGSC0003DMT400075611、PGSC0003DMT400039842、PGSC0003DMT400053553、PGSC0003DMT400053554、PGSC0003DMT400013321、PGSC0003DMT400026257、PGSC0003DMT400069185、PGSC0003DMT400052673和PGSC0003DMT400083968。因此，可以推測HAP組和LAP組的上述10個(gè)基因表達(dá)的膜聯(lián)蛋白、過氧化氫酶、多聚腺苷酸結(jié)合蛋白、伸長因子1α、泛素-40s核糖體蛋白S27a、絲氨酸蛋白酶抑制劑7、60S核糖體蛋白L27、酮醇酸還原異構(gòu)酶、40S核糖體蛋白S3a等蛋白較多，對HAP組高山適應(yīng)性的影響越大。

表1 質(zhì)量預(yù)處理前后數(shù)據(jù)結(jié)果統(tǒng)計(jì)Table 1 Statistics of data results before and after quality preprocessing

注：統(tǒng)計(jì)的數(shù)據(jù)量均是read1+read2的結(jié)果。原始數(shù)據(jù)為各樣本的測序reads數(shù)和堿基總數(shù)；有效數(shù)據(jù)為各樣本經(jīng)質(zhì)量預(yù)處理后的reads數(shù)和堿基總數(shù)；有效數(shù)據(jù)比例為有效reads與原始reads的比值。

Note: The statistical data amount is the result of read1+read2. The original data is the reads number and the total number of bases for each sample. The effective data is the number of reads and the total number of bases after each sample is preprocessed, and the ratio of effective data is the ratio of effective reads to original reads.

表2 reads mapping統(tǒng)計(jì)Table 2 Statistics of reads mapping

表3 表達(dá)量最高的前10個(gè)基因Table 3 The top 10 genes with the highest expression

2.3.2 差異分析 對獲得的差異基因進(jìn)行統(tǒng)計(jì)，結(jié)果表明，以LAP組為對照，HAP組上調(diào)的差異基因有956個(gè)，下調(diào)的差異基因有2 262個(gè)，總數(shù)為3 218個(gè)，應(yīng)用edgeR程序plotSmear函數(shù)對log-fold-changes對log counts-per-million(CPM)進(jìn)行作圖(類似于芯片的MA-plot)，差異的基因用紅色點(diǎn)顯示(圖1)。HAP組和LAP組豐度(RPKM)的分布范圍均為0～5(圖2)。其中，差異最顯著的10個(gè)基因分別為PGSC0003DMT400002035、PGSC0003DMT400026625、PGSC0003DMT400048684、PGSC0003DMT400050249、PGSC0003DMT400007870、PGSC0003DMT400050232、PGSC0003DMT400026364、PGSC0003DMT400003877、PGSC0003DMT400076023和PGSC0003DMT400068129，且10個(gè)基因均為未知蛋白基因(表4)。

表4 差異表達(dá)的前10個(gè)mRNATable 4 The top 10 differentially expressed mRNA

注：紅點(diǎn)表示差異基因。Note:Reddots indicate differential genes.圖1 差異基因分析Fig.1 Differential gene analysis

圖2 基因表達(dá)值分布盒圖Fig.2 Distribution box map of gene expression value

2.3.3 熱圖分析 進(jìn)一步對3 218個(gè)差異基因做雙向聚類熱圖分析，雙向聚類熱圖分析主要用于依據(jù)樣本的相似代謝譜和代謝方式將其聚集在一起。由圖3可知，HAP組和LAP組在整個(gè)表達(dá)譜有明顯的差異，表明二者的代謝譜和代謝方式存在一定的差異，從而影響塊莖的品質(zhì)。

圖3 差異基因熱圖分析Fig.3 Analysis of differential gene heatmap

2.4 GO和KEGG富集分析

將以上得到的差異基因作為前景，全部轉(zhuǎn)錄本作為背景，進(jìn)行GO和KEGG的富集分析，利用超幾何分布算法(phyper)計(jì)算前景轉(zhuǎn)錄本同GO/Pathway分類中某個(gè)特定分支的P值，并用FDR進(jìn)行校正。HAP組和LAP組的差異結(jié)果作GO、KEGG富集分析，結(jié)果表明，在GO分析中，P≤0.05的數(shù)量為28，P≤0.01的數(shù)量為20；在KEGG分析中，P≤0.05和P≤0.01的數(shù)量均為4。

2.5 差異基因的GO富集分析

GO數(shù)據(jù)庫包含了基因參與的生物過程，所處細(xì)胞位置及具有的分子功能三個(gè)方面的功能信息。通過GO中的注釋信息，可以對基因的功能進(jìn)行預(yù)測。Gene Ontology顯著性富集分析以GO term為單位，通過GO功能顯著性富集分析能確定差異表達(dá)基因行駛的主要生物學(xué)功能。通過富集發(fā)現(xiàn)，以LAP組為對照，HAP組的葉綠素結(jié)合、光系統(tǒng)Ⅰ、光合作用和光捕獲、光系統(tǒng)Ⅱ、血紅素結(jié)合、氧化還原酶活性，作用于配對的供體，分子氧的混入或還原等光合系統(tǒng)相關(guān)的GO term富集顯著(表5)

表5 富集的GO termTable 5 Enriched GO term

表5(續(xù))

2.6 差異基因的Pathway富集分析

在生物體內(nèi)，不同基因相互協(xié)調(diào)行使其生物學(xué)功能，通過Pathway顯著性富集能確定差異表達(dá)基因參與的最主要生化代謝途徑和信號(hào)轉(zhuǎn)導(dǎo)途徑。KEGG(Kyoto Encyclopedia of Genes and Genomes)是有關(guān)Pathway的主要公共數(shù)據(jù)庫。Pathway顯著性富集分析以KEGG Pathway為單位。以LAP組為對照，HAP組部分極其顯著的Pathway有光合作用-天線蛋白、丙烷合成、光合作用、苯丙氨酸代謝(表6)。

由圖4可知，HAP組和LAP組與光合作用-天線蛋白、丙烷合成、光合作用、苯丙氨酸代謝等Pathway可以轉(zhuǎn)換光能，提供ATP和H+，用于阿魏酸、芥子酸、芥子醇、松柏醇等次生代謝物的合成。

3 討論

近年來，隨著高通量測序技術(shù)的迅猛發(fā)展，轉(zhuǎn)錄組測序技術(shù)已被廣泛應(yīng)用于生命科學(xué)的研究。目前通過該技術(shù)已經(jīng)完成了諸多植物的轉(zhuǎn)錄功能基因組測序，且一些重要功能基因被挖掘，獲得大量在某一特定組織或細(xì)胞中優(yōu)勢表達(dá)的基因及在某個(gè)生物學(xué)過程中顯著富集的基因[9-10]。此外，轉(zhuǎn)錄組分析還可解析品種之間的差異基因表達(dá)，如黃娟等[11]對甜蕎根和金蕎根進(jìn)行轉(zhuǎn)錄組分析，獲得了4 709個(gè)差異基因，其中 2 460 個(gè)上調(diào)表達(dá)，2 249個(gè)下調(diào)表達(dá)。Zhang等[12]和Ai等[13]對早實(shí)枳和普通枳春梢進(jìn)行轉(zhuǎn)錄組分析，發(fā)現(xiàn)獲得的差異基因與柑橘成花相關(guān)。Zhang等[14]對柱型和普通型蘋果進(jìn)行轉(zhuǎn)錄組分析，發(fā)現(xiàn)了一些新的選擇性剪切位點(diǎn)和轉(zhuǎn)錄本。Ikegami等[15]對無花果栽培種和野生種的果實(shí)進(jìn)行轉(zhuǎn)錄組分析，發(fā)現(xiàn)與性別形成單性結(jié)實(shí)有關(guān)的MADS-box基因、查爾酮合酶基因、赤霉素調(diào)控蛋白基因差異表達(dá)。本研究對懷玉山高山馬鈴薯和懷玉山本土農(nóng)家薯轉(zhuǎn)錄組分析結(jié)果表明，HAP組LAP組上調(diào)的差異基因956個(gè)，下調(diào)的差異基因2 262個(gè)，總數(shù)為3 218個(gè)；HAP組和LAP組差異最顯著的10個(gè)基因均為未知蛋白基因。

表6 富集的KEGG PathwayTable 6 Enriched KEGG Pathway

注：A：sot00196；B:sot00195；C:sot00940；D：sot00360。紅色背景表示注釋上表達(dá)上調(diào)差異基因；藍(lán)色背景表示注釋上表達(dá)下調(diào)差異基因；綠色背景表示只注釋上非差異表達(dá)基因。Note:A：sot00196. B:sot00195. C:sot00940. D：sot00360. The red background represents that the up-expression of differential genes. The blue background represents the down-expression of differential genes. The green background represents non differentially expressed genes.圖4 標(biāo)記差異基因與光合作用相關(guān)PathwaysFig.4 Pathways of marker differential genes related with photosynthesis

轉(zhuǎn)錄組分析也可用于化學(xué)成分合成和生理代謝相關(guān)的基因篩選以及高山低溫干旱脅迫的研究[16-20]，如劉莉等[21]分析東亞中高海拔地區(qū)所特有珍稀瀕危蕨類植物巖穴蕨的轉(zhuǎn)錄組，發(fā)現(xiàn)bHLH、C3H、AP2/ERF、bZIP等轉(zhuǎn)錄因子與高山低溫環(huán)境相關(guān)；黃玉蘭等[22]分析薏苡的轉(zhuǎn)錄組發(fā)現(xiàn)，264條 Unigenes涉及植物信號(hào)轉(zhuǎn)導(dǎo)途徑，其中AMPK(AMP-activated protein kinase)、MAPK(mitogen-activated protein kinase)、Ca2+信號(hào)轉(zhuǎn)導(dǎo)途徑與干旱密切相關(guān)。Cao等[23]對野生種三淺裂野牽牛進(jìn)行轉(zhuǎn)錄組測序分析，挖掘了2 848個(gè)包含bZIP、MYB、bHLH、Orphans、ERF、C2H2、NAC、WRKY等的轉(zhuǎn)錄因子家族，證實(shí)轉(zhuǎn)錄因子與干旱、鹽堿等逆境脅迫密切相關(guān)。寧宇等[24]對青藏高原及其周邊地區(qū)的西藏嵩草進(jìn)行轉(zhuǎn)錄組分析，獲得了西藏嵩草分蘗和抗寒相關(guān)的179個(gè)Unigenes。李錦等[25]建立了新疆天山雪蓮轉(zhuǎn)錄組注解知識(shí)庫，有利于低溫功能基因組學(xué)及低溫耐受機(jī)制的研究。張得芳等[26]對不同海拔下花葉海棠葉片轉(zhuǎn)錄組核苷酸變異進(jìn)行了分析和比較，發(fā)現(xiàn)在低海拔地區(qū)花葉海棠葉片中單核苷酸多態(tài)性的位點(diǎn)數(shù)高于高海拔地區(qū)。本研究發(fā)現(xiàn)HAP組和LAP組與葉綠素結(jié)合、光系統(tǒng)Ⅰ、光合作用和光捕獲、光系統(tǒng)Ⅱ、血紅素結(jié)合、氧化還原酶活性等光合系統(tǒng)相關(guān)的GO term富集顯著；HAP組和LAP組與光合作用-天線蛋白、丙烷合成、光合作用、苯丙氨酸代謝等Pathway富集顯著，說明HAP組和LAP組在懷玉山高海拔生境下主要涉及光合作用基因的差異表達(dá)，其原因可能是高海拔的生境涉及到光合作用基因的表達(dá)，影響一些光合產(chǎn)物的積累，最終影響懷玉山高山馬鈴薯(麻籽洋芋)特殊品質(zhì)的形成，光合作用增強(qiáng)也可能是懷玉山高山馬鈴薯適于高山種植和在平原地區(qū)易于變種的根本原因。轉(zhuǎn)錄組是組織或細(xì)胞轉(zhuǎn)錄出的所有RNA集合[27]，Xu等[28]對暗柳橙及其紅肉芽變品種紅暗柳進(jìn)行轉(zhuǎn)錄組分析，也發(fā)現(xiàn)涉及光合作用的基因差異表達(dá)，同樣認(rèn)為光合作用增強(qiáng)是紅暗柳積累番茄紅素的關(guān)鍵原因。

4 結(jié)論

本研究結(jié)果表明，懷玉山高山馬鈴薯(麻籽洋芋)和懷玉山本土農(nóng)家薯在懷玉山高海拔生境下主要涉及到光合作用基因的差異表達(dá)。但本研究僅研究了懷玉山高山馬鈴薯(麻籽洋芋)轉(zhuǎn)錄組的變化，下一步將對其他高山馬鈴薯轉(zhuǎn)錄組變化進(jìn)行系統(tǒng)研究，以獲得對高山馬鈴薯適應(yīng)高海拔生境的全面認(rèn)識(shí)。