許哲 張晨光 張思遠(yuǎn) 王先裕 高建昌

摘要：植物染色體加倍后的表型變化與基因的表達(dá)密切相關(guān)。本研究對1份二倍體和3份四倍體西瓜材料進(jìn)行轉(zhuǎn)錄組測序并分析差異表達(dá)基因。結(jié)果表明，二倍體和四倍體西瓜的差異表達(dá)基因顯著富集的GO term為細(xì)胞氮化合物生物合成過程、光系統(tǒng)、半胱氨酸型內(nèi)肽酶抑制劑活性、氧化還原酶活性以及類囊體等，KEGG富集分析中發(fā)現(xiàn)差異表達(dá)基因在光合作用、光合作用天線蛋白、葉綠素和卟啉代謝、乙醛酸和二羧酸代謝等通路較為活躍。二倍體和四倍體西瓜中高水平差異基因大都是上調(diào)表達(dá)，氧化應(yīng)激蛋白、脫落酸受體、乙烯響應(yīng)轉(zhuǎn)錄因子等差異表達(dá)基因可能是西瓜多倍體形成的關(guān)鍵因子，對相關(guān)基因的深入研究有助于解析植物倍性變化的分子機(jī)制。

關(guān)鍵詞：多倍體;西瓜;轉(zhuǎn)錄組;差異表達(dá)基因;GO富集分析;KEGG富集分析

中圖分類號： S651.01 ?文獻(xiàn)標(biāo)志碼： A ?文章編號：1002-1302（2020）24-0053-07

西瓜（Citrullus lanatus）原產(chǎn)于非洲，屬于葫蘆科西瓜屬，為一年生蔓生藤本植物，是世界上重要的園藝作物。我國是世界上西瓜生產(chǎn)與消費大國，2016年西瓜播種面積為189.08萬hm2，總產(chǎn)量為7 940萬t[1]。

四倍體西瓜是普通二倍體西瓜染色體加倍的結(jié)果，在田間其最明顯的特征是巨大性，葉片大且肥厚，莖粗且節(jié)間短，葉色濃綠，花和果實也較為巨大，綜合抗病能力較強(qiáng)。染色體加倍獲得四倍體是選育三倍體西瓜（無籽西瓜）的基礎(chǔ)。我國從20世紀(jì)50年代開始對多倍體西瓜誘導(dǎo)方法的研究，如今已建立秋水仙素誘導(dǎo)、體細(xì)胞雜交、胚乳培養(yǎng)等多種四倍體誘導(dǎo)方法;而且建立了染色體計數(shù)、葉肉細(xì)胞分析、流式細(xì)胞儀計數(shù)等染色體倍性檢測方法。利用上述方法選育出很多無籽西瓜品種，為我國西瓜產(chǎn)業(yè)的發(fā)展做出了巨大貢獻(xiàn)。但是，對四倍體西瓜表型變化與染色體倍性的關(guān)系還缺乏分子層面的認(rèn)知。

轉(zhuǎn)錄組是研究細(xì)胞表型和功能的一個重要手段，是指某一特定的生理條件下，細(xì)胞、組織或生物體內(nèi)所有轉(zhuǎn)錄產(chǎn)物的集合，即轉(zhuǎn)錄后所有mRNA的總稱[2]。對二倍體和四倍體西瓜進(jìn)行轉(zhuǎn)錄組水平的分析，為解析四倍體表型的成因提供了新的手段。自2012年西瓜全基因組測序完成后，西瓜轉(zhuǎn)錄組分析亦成為研究熱點。Wechter等利用832個表達(dá)序列標(biāo)簽研究了果實發(fā)育過程中的基因表達(dá)情況[3];Guo等在西瓜果實發(fā)育和成熟過程中鑒定出3 023個差異表達(dá)基因[4];Zhu等通過比較2個紅色和黃色果肉的栽培品種，確定了797個新基因[5];龍婭麗等通過研究二倍體西瓜及其同源四倍體葉片sRNA表達(dá)譜，提供了多倍體抗逆性強(qiáng)的理論依據(jù)[6]。

多倍體化是植物進(jìn)化中的一種普遍現(xiàn)象，也是新物種形成的重要途徑[7]。目前研究者多集中在對二倍體西瓜的研究，關(guān)于二倍體西瓜與四倍體西瓜轉(zhuǎn)錄組差異的研究相對較少。本研究利用Illumina二代高通量測序平臺對二倍體西瓜以及3份四倍體西瓜葉片進(jìn)行轉(zhuǎn)錄組測序分析，比較二倍體西瓜和四倍體西瓜中基因表達(dá)的差異，以期通過轉(zhuǎn)錄組測序（RNA-Seq）技術(shù)獲得與倍性相關(guān)的基因序列，為揭示西瓜的倍性機(jī)制提供參考。

1 材料與方法

1.1 試驗材料

試驗材料為1份二倍體，記為WMA，3份四倍體，記為WMB、WMC和WMD。WMB為WMA染色體加倍獲得，為了增加樣品的遺傳多樣性，加入另外2份不同來源的四倍體材料WMC、WMD，材料均來自中國農(nóng)業(yè)科學(xué)院蔬菜花卉研究所西瓜育種課題組。本試驗未設(shè)生物學(xué)重復(fù)。試驗地點為中國農(nóng)業(yè)科學(xué)院蔬菜花卉研究所日光溫室，2019年6月1日將西瓜種子播種于穴盤中，20 d后長出真葉時取2 cm2大小新鮮葉片，立即放入液氮中保存?zhèn)溆谩?/p>

1.2 總RNA的提取及檢測

使用天根生化科技（北京）有限公司RNA試劑盒對試驗材料進(jìn)行總RNA提取，后通過Nanodrop檢測RNA純度（D260 nm/280 nm、D260 nm/230 nm），通過Agilent 2100對RNA片段長度進(jìn)行檢測。

1.3 cDNA文庫構(gòu)建和測序

樣品檢測合格后，使用帶有Oligo（dT）的磁珠富集mRNA，之后加入破碎緩沖液將mRNA打斷成短片段。再用六堿基隨機(jī)引物（random hexamers）以mRNA為模板進(jìn)行反轉(zhuǎn)錄合成一鏈cDNA，加入緩沖液、dNTPs和DNA聚合酶Ⅰ合成二鏈cDNA。接著，利用AMPure XP beads純化雙鏈cDNA。對純化后的雙鏈cDNA進(jìn)行末端修復(fù)、加A（特異堿基）、加接頭。通過AMPure XP beads核酸純化試劑盒對雙鏈cDNA進(jìn)行片段大小選擇，最后進(jìn)行PCR擴(kuò)增以構(gòu)建cDNA文庫。

文庫構(gòu)建完成后用Agilent 2100對文庫的插入片段大小進(jìn)行檢測，當(dāng)片段大小符合預(yù)期時，再使用定量PCR（Q-PCR）方法對文庫的有效濃度進(jìn)行精確定量（文庫有效濃度>4 nmol/L），以達(dá)到高質(zhì)量文庫標(biāo)準(zhǔn)，最后通過Illumina二代高通量測序平臺，采用PE150雙末端測序策略，完成RNA-Seq。

1.4 數(shù)據(jù)質(zhì)控與基因功能注釋

下機(jī)數(shù)據(jù)（raw data）通常會含有少量的接頭污染及低質(zhì)量的Reads，如果不對其進(jìn)行過濾處理會對后續(xù)分析造成影響，為此我們過濾掉了帶有測序接頭adapter的Reads、N（不確定堿基）含量比例大于10%的Reads以及低質(zhì)量堿基（Q≤20）含量大于50%的Reads。使用TopHat2軟件將數(shù)據(jù)過濾后的Clean Reads與西瓜參考基因組[8]進(jìn)行序列比對，之后用DEGseq[9]進(jìn)行差異表達(dá)分析，再分別與基因本體（GO）、京都基因與基因組百科全書（KEGG）數(shù)據(jù)庫比對，獲得相關(guān)注釋信息。

2 結(jié)果與分析

2.1 轉(zhuǎn)錄組測序數(shù)據(jù)質(zhì)量評估

將二倍體西瓜與四倍體西瓜葉片RNA樣本進(jìn)行測序，結(jié)果見表1。二倍體WMA和四倍體WMB、WMC以及WMD有效數(shù)據(jù)均占原始序列數(shù)據(jù)的9558%以上，Q20（堿基識別錯誤率為1%）堿基百分比均大于97.87%，Q30（堿基識別錯誤率為01%）堿基百分比均大于94.52%，GC含量均大于43.40%，能比對到參考序列的reads百分?jǐn)?shù)均大于95.27%，說明測序結(jié)果有較高的準(zhǔn)確率。

2.2 RNA-Seq相關(guān)性檢查

為了檢測轉(zhuǎn)錄組測序的相關(guān)性，對每個樣品進(jìn)行相關(guān)性分析。相關(guān)性分析是基于樣品間基因整體的表達(dá)水平做的皮爾遜相關(guān)系數(shù)分析，相關(guān)系數(shù)越接近1，表明樣品之間表達(dá)模式的相似度越高。Encode計劃建議皮爾遜相關(guān)系數(shù)的平方（r2）大于0.92（理想的取樣和試驗條件下）。而實際項目操作中，要求r2至少要大于0.8，否則須要對樣品做出合理的解釋，或重新進(jìn)行試驗。本次試驗中WMA、WMB、WMC和WMD 4個樣品之間相關(guān)系數(shù)r2均大于0.899，最高達(dá)到了0.958（圖1），表明樣品間相關(guān)性較好。

2.3 不同倍性西瓜差異表達(dá)基因分析

2.3.1 二倍體WMA與四倍體WMB、WMC和WMD差異表達(dá)基因篩選 對于無生物學(xué)重復(fù)的試驗，為避免引入試驗誤差，應(yīng)該對結(jié)果進(jìn)行嚴(yán)格控制，對差異基因進(jìn)行篩選的閾值一般為：|log2（差異倍數(shù)）|>1 且q值<0.005[10]。結(jié)果（圖2）表明，在WMA和WMB中共檢測到1 458條差異表達(dá)基因，其中表達(dá)上調(diào)的基因有909條，表達(dá)下調(diào)的基因有549條;在WMA和WMC中共檢測到317條差異表達(dá)基因，其中表達(dá)上調(diào)的基因有172條，表達(dá)下調(diào)的基因有145條;在WMA和WMD中共檢測到 1 007 條差異表達(dá)基因，其中表達(dá)上調(diào)的基因有645條，表達(dá)下調(diào)的基因有362條。不同組間差異基因維恩圖顯示，有141個不同組間的共同差異表達(dá)基因。

2.3.2 差異基因的聚類分析 從圖3可以看出，WMB和WMD高表達(dá)量和低表達(dá)量基因模式大體相似。WMA、WMB、WMC和WMD的差異表達(dá)聚合為5類。A組中，WMA和WMD基因表達(dá)量以上調(diào)為主，WMB和WMC基因表達(dá)量以下調(diào)為主;B組中，WMA基因表達(dá)量呈下調(diào)趨勢，WMB、WMC和WMD基因表達(dá)量以下調(diào)為主;C組中WMA和WMC基因表達(dá)量以上調(diào)為主，WMB和WMD基因表達(dá)量以下調(diào)為主;D組中，WMA和WMC基因表達(dá)量以下調(diào)為主，WMB和WMD基因表達(dá)量以上調(diào)為主;E組中，WMA、WMB、WMC和WMD基因表達(dá)量各不一致。二倍體WMA和四倍體WMB、WMD之間的差異表達(dá)基因主要分布在C和D組中，這部分可能調(diào)控西瓜倍性的差異基因。

2.3.3 差異基因GO富集分析 GO是一個用于描述生物體中基因和蛋白質(zhì)的功能分類體系。GO分為分子功能（molecular function）、生物過程（biological process）、和細(xì)胞組成（cellular component）等3個方面。GO的基本單位為term，每個term對應(yīng)一個功能或?qū)傩?。在分析中采用的軟件為GOseq[11]，前30個富集最顯著的GO term見圖4，如果不足30個，則全部展示。結(jié)果表明，WMA和WMB之間富集的差異表達(dá)基因與生物過程、細(xì)胞組成和分子功能相關(guān)的分別有21、3、6個GO term。說明大部分差異基因與生物過程相關(guān)，其中與生物過程相關(guān)富集最顯著、基因數(shù)量最多的term為小分子代謝過程（small molecule metabolic process）和細(xì)胞氮化合物生物合成過程（cellular nitrogen compound biosynthetic process）;與細(xì)胞組成相關(guān)，富集最顯著的term為光系統(tǒng)Ⅰ（photosystem Ⅰ），富集基因數(shù)量最多的類別為光系統(tǒng)（photosystem）;與分子功能相關(guān)富集最顯著的term為半胱氨酸型內(nèi)肽酶抑制劑活性（cysteine-type endopeptidase inhibitor activity），富集基因數(shù)量最多的term為氧化還原酶活性（oxidoreductase activity）。

在WMA和WMC之間，富集的差異表達(dá)基因與生物過程、細(xì)胞組成和分子功能相關(guān)的分別有15、0、15個GO term。其中與生物過程相關(guān)富集最顯著的term為發(fā)病機(jī)制（pathogenesis），富集基因數(shù)量最多的term為胺生物合成過程（amine biosynthetic process）;與分子功能相關(guān)富集最顯著、基因數(shù)量最多的term為內(nèi)肽酶抑制劑活性（endopeptidase inhibitor activity）、內(nèi)肽酶調(diào)節(jié)劑活性（endopeptidase regulator activity）、肽酶抑制劑活性（peptidase inhibitor activity）、肽酶調(diào)節(jié)劑活性和酶抑制劑活性（enzyme inhibitor activity）。

在WMA和WMD之間，富集的差異表達(dá)基因與生物過程、細(xì)胞組成和分子功能相關(guān)的分別有19、5、6個GO term。其中與生物過程相關(guān)富集最顯著、基因數(shù)量最多的term為電子轉(zhuǎn)運（electron transport），其次為細(xì)胞氮化合物生物合成過程（cellular nitrogen compound biosynthetic process）;與細(xì)胞組成相關(guān)，富集最顯著、基因數(shù)量最多的term為類囊體（thylakoid）和類囊體組分（thylakoid part）;與分子功能相關(guān)富集最顯著、基因數(shù)量最多的term為氧化還原酶活性（oxidoreductase activity）。

2.3.4 差異基因的KEGG注釋 KEGG是一個整合了基因組、化學(xué)和系統(tǒng)功能信息的數(shù)據(jù)庫，是系統(tǒng)分析基因產(chǎn)物在細(xì)胞中的代謝途徑以及基因產(chǎn)物功能的數(shù)據(jù)庫[12]。本研究以KEGG代謝途徑數(shù)據(jù)庫為依據(jù)[13]，差異表達(dá)基因顯著富集的前20條通路見圖5。在WMA和WMB之間，差異表達(dá)基因富集最顯著的途徑有光合作用、光合作用天線蛋白、乙醛酸和二羧酸代謝、卟啉和葉綠素代謝等;在WMA和WMC之間，差異表達(dá)基因富集最顯著的途徑有異喹啉生物堿生物合成、?；撬岷蛠喤；撬岽x、托烷和哌啶、吡啶生物堿的生物合成等;在WMA和WMD之間，差異表達(dá)基因富集最顯著的途徑有光合作用天線蛋白、光合作用、卟啉和葉綠素代謝、乙醛酸和二羧酸代謝等。

3 結(jié)論與討論

多倍體化后，隨著染色體的加倍，染色體上每個位點等位基因的數(shù)量也會發(fā)生增倍，并可能由此導(dǎo)致基因表達(dá)質(zhì)和量的變化[7]。本試驗中，通過對二倍體西瓜和四倍體西瓜的差異表達(dá)基因進(jìn)行GO功能注釋，主要GO term為細(xì)胞氮化合物生物合成過程、光系統(tǒng)、半胱氨酸型內(nèi)肽酶抑制劑活性、氧化還原酶活性以及類囊體。Compton等認(rèn)為葉長與葉寬是西瓜倍性水平很好的標(biāo)志，四倍體子房直徑是二倍體子房直徑的1.4倍[14]。蔡力研究發(fā)現(xiàn)，同一發(fā)育時期四倍體紫錐菊葉綠素a、葉綠素b和總?cè)~綠素含量均高于二倍體紫錐菊[15]。本次試驗發(fā)現(xiàn)的差異表達(dá)基因大都與植株的光合作用有關(guān)，一定程度上驗證、解釋了前人的研究，說明四倍體植株形成需要調(diào)控光合作用的相關(guān)基因有較高的表達(dá)量。

通過KEGG顯著性富集分析發(fā)現(xiàn)差異表達(dá)基因主要在光合作用、光合作用天線蛋白、葉綠素和卟啉代謝、乙醛酸和二羧酸代謝等通路中較為活躍，說明它們在四倍體植株形成和生長過程中發(fā)揮重要作用。

多倍體化不僅導(dǎo)致植物基因組大小以及結(jié)構(gòu)發(fā)生改變，還影響了基因的表達(dá)[16]。為了挖掘在二倍體和四倍體西瓜中與倍性相關(guān)的基因，筆者查找了幾個顯著富集的代謝通路下的差異表達(dá)基因，結(jié)合基因定量FPKM，發(fā)現(xiàn)一些在四倍體中表達(dá)量明顯高于二倍體的基因：Cla97C10G205730、Cla97C10G187010、Cla97C01G019450、Cla97C02G026280、Cla97C01G004920和Cla97C07G140200，并對它們進(jìn)行了GO功能注釋。

Cla97C10G187010、Cla97C02G026280、Cla97C01G004920、Cla97C07G140200是調(diào)控細(xì)胞核成分的基因，控制多倍體西瓜細(xì)胞遺傳與代謝。研究發(fā)現(xiàn)西瓜四倍體葉片比二倍體的葉片大而且厚，顏色也較深[17]，Cla97C10G205730和Cla97C01G019450是調(diào)控西瓜細(xì)胞壁、液泡膜、內(nèi)質(zhì)網(wǎng)、高爾基體等膜結(jié)構(gòu)的基因，它們在四倍體西瓜中的表達(dá)量上調(diào)較為明顯，促進(jìn)了細(xì)胞中膜結(jié)構(gòu)的發(fā)育，可能是多倍體西瓜葉片大而厚的主要原因。

另外，對3組對比中差異倍數(shù)達(dá)到2倍以上的差異表達(dá)基因進(jìn)行了注釋分析。結(jié)果表明二倍體和四倍體西瓜中高水平差異基因大都是上調(diào)表達(dá)，主要包括氧化應(yīng)激蛋白、脫落酸受體、乙烯響應(yīng)轉(zhuǎn)錄因子等。這些差異表達(dá)基因可能是西瓜多倍體形成的關(guān)鍵因子。

總之，通過轉(zhuǎn)錄組初步分析，我們較為全面地了解了西瓜二倍體與四倍體的轉(zhuǎn)錄水平變化，發(fā)現(xiàn)了一些與倍性相關(guān)的基因，進(jìn)一步深入研究這些基因的生理功能，有助于揭示多倍體形成的分子機(jī)制，豐富人類對植物多倍性的認(rèn)知。

參考文獻(xiàn)：

[1]中華人民共和國農(nóng)業(yè)部. 中國農(nóng)業(yè)統(tǒng)計資料2016[M]. 北京：中國農(nóng)業(yè)出版社，2017.

[2]姚 娜，劉秀明，董園園，等. 轉(zhuǎn)錄組的測序方法及應(yīng)用研究概述[J]. 北方園藝，2017（12）：192-198.

[3]Wechter W P，Levi A，Harris K R，et al. Gene expression in developing watermelon fruit[J]. BMC Genomics，2008，9：275.

[4]Guo S G，Liu J A，Zheng Y，et al. Characterization of transcriptome dynamics during watermelon fruit development-sequencing，assembly，annotation and gene expression profiles[J]. BMC Genomics，2011，12（1）：454.

[5]Zhu Q L，Gao P，Liu S，et al. Comparative transcriptome analysis of two contrasting watermelon genotypes during fruit development and ripening[J]. BMC Genomics，2017，18：3.

[6]龍婭麗，江雪飛，周 鵬，等. 二倍體西瓜及其同源四倍體葉片sRNA表達(dá)譜分析[J]. 熱帶作物學(xué)報，2018，39（4）：661-668.

[7]王 濤，陳孟龍，劉 玲，等. 植物多倍體化中基因組和基因表達(dá)的變化[J]. 植物學(xué)報，2015，50（4）：504-515.

[8]Guo S G，Zhang J G，Sun H H，et al . The draft genome of watermelon （Citrullus lanatus） and resequencing of 20 diverse accessions[J]. Nat Genet，2013，45：51-58.

[9]Anders S，Huber W. Differential expression analysis for sequence count data[J]. Genome Biol，2010，11：1-12.

[10]Liu Y，Wei H B，Ma M D，et al. Arabidopsis FHY3 and FAR1 regulate the balance between growth and defense responses under shade conditions[J]. The Plant Cell，2019，31（9）：2089-2106.

[11]Young Matthew D，Wakefield M J，Smyth G K，et al. Gene ontology analysis for RNA-seq：accounting for selection bias[J]. Genome Biology，2010，11：R14.

[12]張少平，洪建基，邱珊蓮，等. 紫背天葵高通量轉(zhuǎn)錄組測序分析[J]. 園藝學(xué)報，2016，43（5）：935-946.

[13]Kanehisa M，Arak M，Goto S，et al. KEGG for linking genomes to life and the environment[J]. Nucleic Acids Research，2008，36：480-484.

[14]Compton M E，Barnett N，Gray D J. Use of fluorescein diacetate （FDA） to determimne ploidy of in vitro watermelon shoots[J]. Plant Cell，Tissue and Organ Culture，1999，58：199-203.

[15]蔡 力. 二倍體和四倍體紫錐菊中葉形態(tài)結(jié)構(gòu)及其光合效率的比較研究[D]. 廣東：華南農(nóng)業(yè)大學(xué)，2016：35.

[16]王家利，王 芳，郭小麗，等. 同源多倍體化效應(yīng)研究進(jìn)展[J]. 中國農(nóng)學(xué)通報，2013，29（12）：22-29.

[17]施先鋒，彭金光，李煜華，等. 西瓜多倍體鑒定方法的研究[J]. 浙江農(nóng)業(yè)科學(xué)，2010（2）：273-274.陶亞軍，尹建國，樊繼偉，等. 蘇北地區(qū)稻麥周年生產(chǎn)力品種組合篩選[J]. 江蘇農(nóng)業(yè)科學(xué)，2020，48（24）：60-66.