不同硬度西瓜果皮的轉錄組測序及相關基因表達分析

張敬敬，李 冰，史宇凡，高秀瑞，潘秀清，宋 雪，武彥榮

(河北省農林科學院 經濟作物研究所，河北 石家莊 050051)

西瓜在我國農業(yè)生產中占有重要地位，我國是西瓜的主要產區(qū)之一，每年西瓜產量約7 000萬t[1]，占全球西瓜產量的67%[2]。河北省西瓜呈現區(qū)域化種植，優(yōu)勢產區(qū)明顯，自“十二五”以來，衡水阜城和保定清苑西瓜播種面積在6 666.67 hm2以上[3]，西瓜生產區(qū)域化種植、上市時間集中，對品種的耐裂性、貯藏性、貨架期品質提出了更高的要求。

西瓜果皮硬度在果實的耐裂性、貯藏性、貨架期品質等方面起著重要作用[4]，生產上許多西瓜品種存在優(yōu)質不耐裂、耐裂不優(yōu)質現象，是造成西瓜生產過程中和采后產量損失的主要原因。通過對果皮硬度研究，在保證產量和品質的基礎上改良果皮硬度，創(chuàng)制優(yōu)質、果皮高硬度、耐裂的種質資源，在培育耐裂西瓜新品種上意義重大。

西瓜裂果是造成田間和產后減產的主要原因，果實開裂是一個復雜的性狀，遺傳和環(huán)境因素對西瓜裂果影響較大[5]，1994年趙尊練等[4]首次發(fā)現西瓜果皮硬度是果實耐貯性的主要內在因素，之后學者們分別從生理生化[6]、細胞水平[7-8]、分子水平[9]對西瓜果皮硬度進行研究，不同硬度西瓜的果皮含水量、纖維素、半纖維素、果膠含量、細胞排列方式、石細胞和外果皮層數相差較大，Liao 等[9]采用BSA-seq 技術和SNP 分析對西瓜果皮硬度進行研究，發(fā)現CIERF4基因與果皮硬度密切相關，該基因可能參與木質素的生物合成和細胞壁修飾，從而調控果皮硬度。而隨著西瓜全基因組測序完成[10]，促進了轉錄組測序技術在西瓜果肉硬度[11-12]、果皮色澤[13]等方面的應用，但是，轉錄組測序技術在西瓜果皮硬度方面的研究鮮有報道。

本試驗以河北省農林科學院經濟作物研究所瓜類室創(chuàng)制的西瓜果皮高硬度的耐裂資源(901)4-1-1-M和低硬度易裂資源BSH為研究對象，通過對果皮硬度差異最大時期果皮的RNA-seq轉錄組測序，分析調節(jié)果皮硬度高的主要生物學過程和代謝通路，探討果皮硬度差異大的基因表達調控機制，鑒定與果皮硬度相關的候選基因，從轉錄組水平解釋西瓜果皮硬度差異的形成機制，旨在為耐貯運西瓜品種培育提供理論依據。

1 材料和方法

1.1 試驗材料

供試西瓜為高硬度耐裂材料(901)4-1-1-M和低硬度易裂材料BSH，由河北省農林科學院經濟作物研究所選育而成。西瓜開花后掛牌標記授粉日期，分別取授粉后第10，14，18，22，26，30天的西瓜，取樣時選取果實大小一致的果實，使用質構儀測定西瓜陽面果皮硬度，3 次生物學重復；選擇檢測部位兩側切割大小為3 cm×3 cm的果皮，厚度1 cm，每個生物學重復包含一個西瓜的2塊西瓜果皮，3 次生物學重復，取樣后用錫箔紙包裹后迅速放入液氮，保存于-80 ℃冰箱備用。選取果實成熟果皮硬度差異最大時期即授粉后第30天的西瓜果皮，用于后續(xù)的轉錄組測序和qRT-PCR驗證。

1.2 果皮轉錄組分析

委托北京諾禾生物科技有限公司完成RNA提取、文庫構建和轉錄組測序。以(901)4-1-1-M和BSH陽面果皮為轉錄組分析材料，利用Illumina HiSeqTM測序平臺進行轉錄組測序分析，以西瓜97103基因組為參考基因組進行序列比對，對獲得的差異基因進行GO富集分析、KEGG富集分析和差異基因功能分析。

1.2.1 RNA提取和文庫構建 以(901)4-1-1-M和BSH陽面果皮為轉錄組分析材料，采用天根DP441多糖多酚植物總RNA提取試劑盒提取西瓜果皮總RNA，使用Agilent 2100 bioanalyzer精確檢測RNA完整性和總量，RNA質控合格后通過Oligo(dT)磁珠富集帶有polyA尾的mRNA，再用二價陽離子將得到的mRNA隨機打斷，以片段化的mRNA為模板、隨機寡核苷酸為引物，在M-MuLV逆轉錄酶體系中合成cDNA第1條鏈，隨后用RNaseH降解RNA鏈，并在DNA polymerase Ⅰ體系下，以dNTPs為原料合成cDNA第2條鏈。純化后的雙鏈cDNA經過末端修復、加A尾并連接測序接頭，用AMPure XP beads篩選250～300 bp的cDNA進行PCR擴增，并再次使用AMPure XP beads純化PCR產物，最終獲得文庫。

1.2.2 轉錄組測序 文庫構建完成，庫檢合格后，利用Illumina HiSeqTM測序平臺進行轉錄組測序分析，在測序的flow cell中加入4種熒光標記的dNTP、DNA聚合酶以及接頭引物進行擴增，在每一個測序簇延伸互補鏈時，每加入一個被熒光標記的dNTP就能釋放出相對應的熒光，測序儀通過捕獲熒光信號，并通過計算機軟件將光信號轉化為測序峰，從而獲得待測片段的序列信息。

1.2.3 測序結果分析 對測序獲得的原始數據去除接頭、無法確定的堿基信息、低質量的信息，獲得高質量的測序數據(Clean data)，以西瓜97103基因組為參考基因組進行有參轉錄組比對分析，采用Subread軟件對基因表達水平進行定量，利用DESeq 2軟件進行差異基因分析，以|log2(FoldChange)|≥1 & padj≤0.05為標準篩選差異基因。采用ClusterProfiler軟件對差異基因集進行GO功能富集分析和KEGG通路富集分析。

1.3 實時熒光定量qRT-PCR驗證

選擇與果皮硬度發(fā)育相關的9個差異基因進行實時熒光定量PCR(qRT-PCR)驗證，采用天根DP441多糖多酚植物總RNA提取試劑盒提取(901)4-1-1-M和BSH西瓜果皮的RNA，再用HiFiScript gDNARemoval cDNA Syntheris Kit反轉錄試劑盒進行反轉錄合成cDNA后作為模板，最后采用ChemoHS qPCR Mix(None ROX)試劑，以西瓜的Actin基因作為內參基因，在CFX connectTM(BIO-RAD Applied Biosystems，Foster City，CA，USA)實時熒光定量PCR儀上進行qRT-PCR，每個基因設置3次生物學重復，按以下程序進行：95 ℃，15 min；40個PCR循環(huán)(95 ℃，10 s；58 ℃，30 s；72 ℃，30 s)。熒光定量數據采用2-ΔΔCt相對定量計算，各引物序列見表1。

表1 qRT-PCR引物序列Tab.1 The primer sequence of qRT-PCR

2 結果與分析

2.1 質構儀測定結果

不同時期的(901)4-1-1-M和BSH西瓜果皮硬度試驗結果表明(圖1)，(901)4-1-1-M果皮硬度在各個時期均極顯著高于BSH(P<0.01)，果實授粉后第18天(901)4-1-1-M果皮硬度達最大值89.02 N，BSH果皮硬度為42.46 N；授粉后30 d(901)4-1-1-M果皮硬度為73.77 N，比BSH果皮硬度高47.21 N，此時果皮硬度差異最大。

A、B為在P<0.01條件下差異顯著。A and B means the significant difference under P<0.01.

2.2 轉錄組測序質量評估

經轉錄組測序可知，6個樣品共得到275 043 248條原始序列(Raw reads)，為保證信息分析質量，去除Raw reads里面含有帶接頭的、低質量的reads，得到濾后序列(Clean reads)272 077 872條，后續(xù)分析都基于Clean reads；20(Q20)和30(Q30)分別表示堿基識別出錯的概率為1%和0.1%，Q20和Q30值均達到90%以上；GC百分率分布檢查用于檢測有無AT、GC 分離現象，本試驗中西瓜果皮的GC百分率均在43%以上(表2)。

表2 6個樣本轉錄組數據測序質量Tab.2 The sequencing quality of 6 sample transcriptome data

將Clean reads 與西瓜參考基因組比對，由比對結果可知(表3)，完全比對到參考基因組的reads占總reads數的71%以上，唯一比對位置序列比對率在70%以上，多個比對位置序列比對率在1.8%以下，外顯子比對率在75%以上。據此可知，此次測序質量能夠滿足試驗要求。

表3 與參考基因組比對結果Tab.3 Sequence comparison of samples with reference genome and genes

2.3 基因表達量分析

通過西瓜樣本間相關性熱圖(圖2)，對樣本間基因表達量的生物學重復進行評估，樣品間基因表達水平相關性是檢驗試驗可靠性和樣本選擇是否合理的重要指標，相關系數R2越接近1，表明樣品之間表達模式的相似度越高，由圖3可知，6個西瓜樣本的R2最小值為0.889，滿足試驗的生物學重復要求。

圖2 6個西瓜樣本間相關性分析Fig.2 The correlation analysis among 6 watermelon samples

2.4 基因表達數據分析

西瓜(901)4-1-1-M和BSH的共表達基因16 229個，(901)4-1-1-M特有的表達基因431個，BSH所特有的表達基因390個，采用DESeq 2軟件對基因表達進行差異分析，以|log2(FoldChange)|≥1 & padj≤0.05為標準，篩選差異表達基因(DEG)。由圖3可知，(901)4-1-1-M和BSH西瓜果皮中共篩選出1 085個差異表達基因，其中有555個上調表達基因，530個下調表達基因。

圖3 差異基因火山圖Fig.3 Volcanic map of differentially expressed genes

2.5 DEGs的GO富集分析

利用GO(Gene Ontology)數據庫對獲得的1 085個差異基因進行富集分析，1 085個差異基因被注釋在生物學過程(Biological process)、細胞組分(Cellular component)及分子功能(Molecular function)三大功能方面，分別涉及254，41，213條功能分支。對富集程度前30的GO 條目作柱形圖(圖4)，差異基因顯著富集在DNA復制(GO：0006260)(13個)、細胞壁(GO：0005618)(13個)、細胞邊緣(GO：0071944)(16個)、外部封裝結構(GO：0030312)(13個)、細胞外區(qū)(GO：0005576)(7個)、質外體(GO：0048046)(5個)、酶抑制劑活動(GO：0004857)(13個)、血紅素結合(GO：0020037)(31個)、四吡咯結合(GO：0046906)(31個)、氧化還原酶活性，作用于配對供體(GO：0016705)(27個)、鐵離子結合(GO：0005506)(26個)、轉移酶活性，轉移己糖基(GO：0016758)(30個)、果膠酯酶活性(GO：0030599)(8個)、蛋白質二聚化活性(GO：0046983)(27個)、酶調節(jié)劑活性(GO：0030234)(14個)和DNA結合(GO：0003700)(35個)方面。

2.6 DEGs的KEGG富集分析

利用KEGG數據庫，1 085個差異基因被富集到81個代謝通路中，選擇富集程度前20條代謝通路進行散點作圖(圖5)，差異基因主要在苯丙烷代謝(19個)，DNA復制(8個)，二萜生物代謝(5個)，角質、栓木和蠟的生物合成(5個)，戊糖、葡萄糖醛酸轉換(8個)，α亞麻酸代謝(6個)，淀粉和蔗糖代謝(12個)，植物激素信號轉導(17個)等通路富集，其中苯丙烷生物代謝(cmo00940)、DNA復制(cmo03030)、二萜生物合成(cmo00904)通路顯著性富集(padj<0.05)。苯丙烷生物代謝途徑是顯著富集的代謝通路中注釋到差異表達基因數目最多的通路，說明苯丙烷代謝物在西瓜果皮硬度差異方面存在重要作用。

圖5 1 085個DEGs的前20條KEGG的功能分類Fig.5 The 20 tops of Enriched KEGG analysis for 1 085 DEGs

2.7 苯丙烷生物代謝中DEGs分析

對1 085個差異基因進行分析，(901)4-1-1-M和BSH在苯丙烷生物合成中存在19個差異基因(表4)。通過分析苯丙烷代謝通路發(fā)現(圖6)，下游代謝物質木質素、對羥基苯酚木質素、愈創(chuàng)木酚木質素和5羥基愈創(chuàng)木酚木質素信號通路中富集到7個差異基因上調表達(Cla97C01G025380、Cla97C03G055890、Cla97C04G070210、Cla97C05G097050、Cla97C10G203590、Cla97C02G045070、Cla97C02G049940)，8個差異基因下調表達(Cla97C09G165820、Cla97C02G044950、Cla97C04G075800、Cla97C09G165630、Cla97C10G195660、Cla97C07G139710、Cla97C07G144540、Cla97C02G032910)，其中，下游基因Cla97C01G025380在(901)4-1-1-M果皮中基因表達量是BSH的60.24倍，而該基因在BSH中基因表達量最低，平均含量為0.056(表4)。

圖6 苯丙烷生物合成途徑中的差異基因表達Fig.6 DEGs in phenylpropanoid biosynthesis pathways

表4 苯丙烷生物合成路徑差異基因表達量Tab.4 Expression analysis of genes related to phenylpropanoid biosynthesis

2.8 差異基因功能分析

從苯丙烷代謝路徑中篩選出的19個差異基因，其中，6個過氧化氫酶基因(Cla97C02G049940、Cla97C05G097050、Cla97C03G055890、Cla97C02G045070、la97C04G070210)，4個轉移酶基因(Cla97C02G033110、Cla97C09G165630、Cla97C07G144540、Cla97C02G032910)，3個糖基水解酶基因(Cla97C11G210660、Cla97C06G123590、Cla97C11G224410)，2個輔酶A連接酶基因(Cla97C09G165820、Cla97C10G195660)，脫氨酶、羥化酶、磷脂酶和橋小檗堿酶各1個，依次為Cla97C04G075800、Cla97C02G044950、Cla97C07G139710、Cla97C01G025380。

2.9 實時熒光驗證

為了驗證轉錄組數據的準確性，從富集程度最高的苯丙烷代謝路徑中選擇與果皮硬度相關的差異表達基因進行實時熒光驗證(qRT-PCR)，對轉錄組測序數據(RNA-seq)和qRT-PCR的結果進行相關性分析，相關系數R2為0.791，證明RNA-seq結果試驗數據可靠(圖7)。

圖7 RNA-seq和qRT-PCR相關性分析Fig.7 The correlation analysis results between RNA-seq and qRT-PCR

驗證基因中包括與木質素代謝物質形成過程相關的9個差異基因，苯基丙氨酸在Cla97C04G075800、Cla97C02G044950、Cla97C09G165820、Cla97C10G195660基因作用下形成香豆酰輔酶A，再在Cla97C09G165630、Cla97C07G139710、Cla97C02G032910、Cla97C07G144540、Cla97C02G045070基因作用下最終形成不同種類的木質素；qRT-PCR結果顯示，下游基因Cla97C02G045070在(901)4-1-1-M中表達量較BSH提高7.13倍，這與RNA-seq結果中富集到苯丙烷通路的基因上調表達結果一致。

3 結論與討論

果實硬度變化是自然界常見的生理生化現象[14]，近年來，對果實硬度形成機理研究較多如蘋果[15-16]、李子[17]、藍莓[18]、西瓜[19-20]等，但對果皮硬度研究較少，Liao 等[9]采用BSA-seq 技術對西瓜果皮硬度進行研究，將果皮硬度基因精細定位到9 kb 的CIERF4(Cla97C10G187120)基因，該基因可能參與木質素的生物合成和細胞壁修飾從而調控果皮硬度。秦改花等[21]通過對石榴籽粒硬度與木質素、纖維素含量的相關性研究發(fā)現，籽粒木質素、纖維素含量與籽粒硬度顯著相關，其可作為評價石榴籽粒硬度的指標。Qin等[22]通過對軟硬籽石榴不同時期的種皮轉錄組、基因組和代謝組聯合分析，篩選出與硬籽石榴中木質素單體芥子醇大量積累的相關候選基因POD(Pgr011171.1)和木質素單體松柏醇積累的相關候選基因F5H(Pgr013634.1)、CAD(Pgr022328.1，Pgr006310.1，Pgr006318.1和Pgr006319.1)。劉戀等[23]通過采用質構儀、生理生化指標測定和轉錄組分析發(fā)現，高硬度金桔果皮中木質素含量顯著高于低硬度金桔，轉錄組測序結果表明，苯丙烷生物合成代謝引起的木質素含量差異導致了3個品種金桔果皮韌性存在差異。

本試驗采用RNA-seq 技術分析了果皮硬度差異顯著的2份資源(901)4-1-1-M和BSH的果皮轉錄本差異，與果皮低硬度易裂資源BSH相比，(901)4-1-1-M中555個基因上調表達，530個基因下調表達；GO富集分析發(fā)現，1 085個差異表達基因顯著富集在細胞壁、細胞外圍、胞外區(qū)、四吡咯骨架、氧化還原酶活性、果膠酯酶活性等功能類別；KEGG富集分析發(fā)現，19個差異表達基因顯著富集在苯丙烷代謝通路，其中下游基因Cla97C01G025380在(901)4-1-1-M中顯著上調表達，導致木質素單體芥子醇、松柏醇、5羥基松柏醇形成，進而在Cla97C03G055890、Cla97C04G070210、Cla97C05G097050、Cla97C10G203590、Cla97C02G045070、Cla97C02G049940基因作用下形成木質素、5羥基愈創(chuàng)木酚木質素、愈創(chuàng)木酚木質素和對羥基苯酚木質素代謝物。研究結果與Liao 等[9]發(fā)現的西瓜果皮硬度與木質素的生物合成路徑相關的結論一致；與Qin等[22]發(fā)現的硬籽石榴中木質素單體芥子醇大量積累和木質素單體松柏醇積累的研究結果一致；與劉戀等[23]發(fā)現的金桔果皮硬度與木質素正相關，果皮韌性差異富集在苯丙烷代謝通路的結論一致。

因此，推測苯丙烷代謝通路中與木質素代謝物形成相關的差異表達基因可能是造成西瓜果皮硬度差異的主要原因，但與Liao 等[9]發(fā)現的CIERF4(Cla97C10G187120)基因參與木質素代謝結果不一致，主要原因可能是由于西瓜材料不同，但本試驗中發(fā)現，下游基因Cla97C01G025380在BSH中基因表達量低，影響芥子醇、松柏醇、5羥基松柏醇代謝物形成，這與Qin等[22]研究結果一致，本試驗中發(fā)現的西瓜果皮硬度富集到苯丙烷代謝路徑中木質素與劉戀等[23]的研究結果一致。

本研究通過對西瓜果皮硬度差異顯著的2種資源(901)4-1-1-M和BSH成熟期的果皮進行轉錄組分析，發(fā)現了一些與西瓜果皮硬度密切相關的代謝路徑，尤其是苯丙烷代謝途徑最終形成的不同種類的木質素在西瓜果皮硬度形成中起重要作用，其中，19個差異基因被富集到該通路，包括Cla97C04G075800、Cla97C02G044950、Cla97C10G195660、Cla97C09G165820、Cla97C09G165630、Cla97C02G032910、Cla97C07G144540、Cla97C02G045070等基因。該研究從轉錄水平解析了(901)4-1-1-M和BSH西瓜果皮硬度存在差異的分子機制，為挖掘西瓜果皮硬度相關的關鍵基因提供了理論基礎。