PEG 模擬干旱條件下黃瓜種子發芽前后轉錄組分析

謝洋 周國彥 蘇航 邢雨蒙 閆立英

（河北科技師范學院園藝科技學院 河北省高校特色園藝植物生物育種應用技術研發中心，秦皇島 066004）

黃瓜（Cucumis sativusL.）是世界性蔬菜作物，屬于葫蘆科黃瓜屬一年生蔓生草本植物［1］。受全球氣候變暖的影響，我國干旱、半干旱化的耕地面積持續增加，全國大部分地區干旱水平呈上升趨勢［2］。黃瓜根系淺（吸水能力弱）、葉片大（蒸騰作用強），喜濕且不耐旱［3］。干旱脅迫對黃瓜的生長發育、產量和品質形成會產生諸多不良影響，如種子發育不良［4］、植株萎蔫［5］、開花結果異常［6］，甚至死亡［7］等現象。因此，解析黃瓜干旱脅迫響應分子作用機制對于改良其品質和產量具有重要意義。

一定干旱條件下，植物體自身具有維持細胞水分平衡從而保持生存的策略。一方面，植物受到干旱脅迫后，體內產生活性氧，當活性氧超過植物體可處理能力后就會啟用抗氧化防御系統，主要是超氧化物歧化酶（SOD）、過氧化物酶（POD）和過氧化氫酶（CAT）等抗氧化酶的調節作用，通過幾種抗氧化酶的動態變化清除植物體內過量的活性氧，以維持植物體自身平衡［8］；另一方面，干旱脅迫還會對植物產生直接的水分脅迫，而植物細胞會通過積累滲透調節物質如K+、Na+、Ca2+等無機離子，脯氨酸、甜菜堿、可溶性糖和多元醇等有機溶質及誘導產生的滲透調節蛋白等以降低細胞的滲透勢和防止細胞脫水［9］。此外，ABA 途徑是調節植物干旱響應并合理利用水分的重要策略，干旱脅迫促使植物器官中的ABA 產生和積累進而激活下游信號傳導［10］。

植物對干旱脅迫的應答是一個涉及多基因、多信號途徑和代謝過程的復雜反應機制［11］。轉錄組測序可以從mRNA 水平全面地研究植物基因功能及特定的生物學過程，為園藝植物重要性狀基因挖掘和分子標記開發提供強大技術支持［12］。

本研究以抗旱型旱黃瓜‘KS30’與敏感型旱黃瓜‘KS30’為試材，分別通過種子萌發期的PEG 模擬干旱誘導進行發芽前期和發芽后期取樣，構建8個cDNA 文庫，利用轉錄組測序技術和生物信息學分析技術挖掘旱黃瓜抗旱相關基因，并基于PEG 模擬干旱誘導下抗旱型與敏感型黃瓜轉錄組數據開發全轉錄組水平的SNP 分子標記，進一步通過抗旱相關差異基因與SNP 變異位點信息，篩選旱黃瓜抗旱相關SNP 標記。

1 材料與方法

1.1 材料

本研究選用的抗旱型旱黃瓜‘KS30’（DT）和敏感型旱黃瓜‘KS30’（DS）為雜種一代品種。

1.2 方法

1.2.1 樣品準備與轉錄組測序 將DT 和DS 種子分別進行去離子水（CK）和15% PEG6000 溶液（T）處理36 h 后進行第1 批次取樣，樣品命名為DT1?T、DT1?CK、DS1?T、DS1?CK；繼續觀察12 h 后進行第2 批次取樣，命名為DT2?T、DT2?CK、DS2?T、DS2?CK。樣品取樣后立即存入超低溫冰箱，備用。每個處理15 粒種子，3 次生物學重復。

根據總RNA 提取試劑盒（TIANGEN）操作說明，分別提取8 份樣品總RNA。利用瓊脂糖凝膠電泳分析RNA 降解程度和是否有污染，Nanodrop檢測RNA 純度（OD260/280）。總RNA 樣品檢測合格后，用磁珠寡核苷酸Oligo（dT）（Invitrogen）提取mRNA 并用裂解緩沖液將其片段化；以mRNA 為模板，用隨機引物合成第一條cDNA 鏈，隨后加入buffer、dNTPs 和DNA 聚合酶I 合成第二條cDNA。雙鏈cDNA 純化和末端修復以及3'末端單核苷酸A（腺嘌呤）添加；將此片段連接到測序接頭上，通過PCR 擴增富集獲得最終的cDNA 文庫；最后，利用Illumina HiSeqTM2500 平臺進行文庫測序［13］。RNA提取、cDNA 文庫構建與上機測序委托北京諾禾致源科技股份有限公司。

1.2.2 參考基因組比對與新基因預測 將測序獲得的原始數據raw reads 經過濾、測序錯誤率檢查、GC含量分布檢查，獲得后續分析使用的clean reads。使用HISAT2 軟件將clean reads 與參考基因組進行快速精確的比對，獲取reads 在參考基因組上的定位信息［14］。

1.2.3 定量分析與差異基因鑒定 根據基因比對在參考基因組上的位置信息，采用subread 軟件中的featureCounts 工具統計每個基因從起始到終止范圍內覆蓋的readscount。采用FPKM（expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced）方法進行基因表達定量［15］。基因表達定量完成后，對其表達數據進行統計學分析，獲得假設檢驗概率（P?value）和多重假設檢驗校正FDR 值（Padj），篩選樣本在不同狀態下表達水平顯著差異的基因。利用edgeR 軟件，篩選標準設定為|log2（FoldChange）|≥2 且Padj ≤0.05［16］。

1.2.4 差異基因韋恩圖與富集分析 利用在線網站(https://bioinfogp.cnb.csic.es/tools/venny/index.html)制作不同比較組合間的韋恩圖，展示比較組合共有或獨有的差異基因。采用clusterProfile 軟件分別對差異基因集進行GO（gene ontology）功能富集和KEGG（Kyoto encyclopedia of genes and genomes）通路富集分析，均以Padj 小于0.05 作為顯著性富集的閾值，并分別選取最顯著的20 個通路繪制散點圖進行結果的展示［17-18］。

1.2.5 SNP 變異位點檢測與分析 利用PEG 模擬干旱誘導下抗旱型與敏感型黃瓜轉錄組數據開發全轉錄組水平的SNP 分子標記。使用GATK 進行變異位點的檢測，進而根據SnpEff 注釋信息對每個變異位點進行變異位點功能統計（同義突變、錯義突變、無義突變）、變異位點區域統計（EXON、INTRON、INTERGENIC）以及變異位點影響統計（HIGH、MODERATE、LOW、MODIFIER）等方面進行分析和繪圖［19］。聯合轉錄組挖掘到的重要候選基因和SNP 位點檢測信息，篩選抗旱相關SNP 分子標記。

1.2.6 熒光定量PCR 分析 選取轉錄組中差異倍數大、表達豐度高且含有SNP 變異位點的基因進行熒光定量PCR（RT?qPCR）分析。采用Beacon Designer 7.0 軟件設計RT?qPCR 引物，內參基因序列信息參考前人研究報道［20］，均由中美泰和生物技術（北京）有限公司合成。總RNA 提取（型號DP432）、cDNA反轉錄（型號KR116）與RT?qPCR 加樣與擴增程序（型號FP313）均參考天根生化科技（北京）有限公司試劑盒說明書，每個基因3 次重復。利用美國伯樂bio?rad CFX connect 實時熒光定量PCR 儀進行擴增，基因相對表達量計算采用2-??CT方法，導出數據經Excel 和DPS 軟件作圖與差異顯著性分析。

2 結果

2.1 PEG模擬干旱誘導下黃瓜種子萌發特征分析

從圖1 可知，PEG 模擬干旱誘導下，敏感型（DS）旱黃瓜種子萌發力較對照變差，其芽長和芽長+根長分別下降90.90%和42.29%；而抗旱型（DT）旱黃瓜較對照其種子萌發力增強，其芽長和芽長+根長分別增加93.33%和24.50%（圖1）。結果表明，基因型的差異導致抗旱型（DT）和敏感型（DS）旱黃瓜種子響應干旱的機制不同，且抗旱型（DT）更能適應干旱脅迫環境。

圖1 PEG 模擬干旱誘導下旱黃瓜種子萌發特征Fig. 1 Seed germination characteristics of dry cucumber under PEG drought simulation

2.2 測序數據質量分析

從8 個樣本文庫DT1_CK、DT1_T、DS1_CK、DS1_T、DT2_CK、DT2_T、DS2_CK、DS2_T 中分別獲 得43.83、39.02、40.69、42.27、40.79、43.82、41.55 和43.06 M 的clead reads， 與raw reads 比值介于97.94%-99.18%，且其Q20 值分別為97.20%、97.58%、97.10%、97.30%、97.13%、97.23%、97.34%和97.22%，說明測序reads 成功比對到基因組的比例非常高，能夠獲得高比例的reads 在參考基因組上的定位信息（表1）。

表1 樣本測序數據質量匯總Table 1 Quality summary of samples sequencing data

2.3 差異表達基因鑒定

從4 組比較組差異基因數量可看出，抗旱型材料發芽前期PEG 模擬干旱誘導較對照（DT1_TvsDT1_CK）鑒定到727 個上調和345 個下調差異基因，敏感型材料發芽前期PEG 模擬干旱誘導較對照（DS1_TvsDS1_CK）鑒定到1 226 個上調和111個下調差異基因；抗旱型材料發芽后期PEG 模擬干旱誘導較對照（DT2_TvsDT2_CK）鑒定到333 個上調和480 個下調差異基因，敏感型材料發芽后期PEG 模擬干旱誘導較對照（DS2_TvsDS2_CK）鑒定到600 個上調和470 個下調差異基因。同一時期，抗旱型較敏感型的差異基因數量偏少，并且在上調基因數量上尤為明顯。如發芽前期抗旱型誘導與對照相比，上調基因727 個，下調基因345 個；而敏感型上調基因1 226 個，是抗旱型1.69 倍。同樣，發芽后期敏感型上調基因是抗旱型的1.80 倍（圖2）。結果說明，敏感型旱黃瓜對PEG 模擬干旱誘導響應表現更為劇烈。

圖2 比較組合差異基因數目統計柱狀圖Fig. 2 Statistical bar chart for comparison of the number of combinatorial differential genes

2.4 差異基因韋恩圖分析

發芽前期抗旱型和敏感型2 個材料，PEG 模擬干旱誘導較處理組共有差異基因數量為89 個，分別有212 和1 202 個差異基因特異于抗旱型和敏感型比較組中。在發芽后期抗旱型和敏感型2 個材料，PEG 模擬干旱誘導較處理組共有差異基因數量為132 個，分別有681 和938 個差異基因特異于抗旱型和敏感型比較組中（圖3）。從差異基因數量來看，發芽前期響應PEG 模擬干旱誘導的差異基因在抗旱型和敏感型旱黃瓜中差異數目較懸殊，可作為后續抗旱相關基因挖掘的重要數據基礎。

圖3 比較組合差異基因韋恩圖Fig. 3 Venn diagram of comparing the combined differential genes

2.5 差異基因富集分析

發芽前期抗旱型材料，PEG 模擬干旱誘導較處理組差異基因主要顯著富集的GO 條目有氧化還原過程（GO：0055114）、細胞壁（GO：0005618）、氧化還原酶活性（GO：0016491）、陽離子結合（GO：0005506）等；發芽期敏感型材料，PEG 模擬干旱誘導較處理組差異基因主要顯著富集的GO 條目有氧化還原過程（GO：0055114）、響應刺激（GO：0050896）、膜（GO：0016020）、氧化還原酶活性（GO：0016491）、轉運活性（GO：0005215）、DNA 結合轉錄因子活性（GO：0003700）等（圖4）。

圖4 GO 富集分析柱狀圖Fig. 4 Bar diagram of GO enrichment analysis

發芽前期抗旱型材料，PEG 模擬干旱誘導較處理組差異基因主要顯著富集的KEGG 通路有苯丙素的生物合成（csv00940）、戊糖、葡萄糖醛酸轉換（csv00040）和亞油酸代謝（csv00591）（圖5）；發芽期敏感型材料，PEG 模擬干旱誘導較處理組差異基因主要顯著富集的KEGG 通路有苯丙素的生物合成（csv00940）、植物激素信號轉導（csv04075）、戊糖、葡萄糖醛酸轉換（csv00040）、α-亞油酸代謝（csv00592）、玉米素生物合成（csv00908）和亞油酸代謝（csv00591）等（圖5）。

圖5 KEGG 富集分析散點圖Fig. 5 Point diagram of KEGG enrichment analysis

將發芽前期的抗旱型與敏感型2 個材料，PEG模擬干旱誘導與對照的差異表達基因進行關聯分析，篩選到亞油酸代謝通路中3 個，分別是CsaV3_4G 023820、CsaV3_2G006420 和CsaV3_4G023810；戊糖、葡萄糖醛酸轉換（csv00040）通路中4 個，分別 是CsaV3_2G025090、CsaV3_3G011280、CsaV3_3G028700 和CsaV3_2G014800；苯丙素的生物合成（csv00940）通路中6 個，分別是CsaV3_1G042480、CsaV3_6G006890、CsaV3_3G030410、CsaV3_2G03 5150、CsaV3_2G009070 和CsaV3_6G039690；以及植物激素信號轉導（csv04075）通路中7 個，分別是CsaV3_2G013210、CsaV3_5G023170、CsaV3_6G 007970、CsaV3_7G027610、CsaV3_1G030250、Csa V3_2G004130 和CsaV3_3G005590（表2）。

表2 KEGG 顯著富集的差異基因分析Table 2 Differential gene analysis of KEGG enrichment

2.6 基于轉錄組SNP分子標記開發與篩選

8 個轉錄組文庫DT1_CK、DT1_T、DS1_CK、DS1_T、DT2_CK、DT2_T、DS2_CK、DS2_T 分別檢測到外顯子區域的變異數量為5 018、5 101、7 340、7 972、5 252、5 819、8 471 和8 228 個；檢測到同義突變、錯義突變、無義突變3 種突變類型，其中，無義突變數量分別為25、26、29、26、24、27、29和26 個，錯義突變數量分別為2 114、2 166、2 898、3 173、2 132、2 375、3 315 和3 211 個，同義突變數量分別為2 905、2 936、4 441、4 810、3 118、3 443、5 175 和5 033 個；變異位點影響高的數量分別為41、43、49、48、40、47、49 和44 個（圖6）。

圖6 SNP 變異位點統計分析Fig. 6 Statistical analysis of SNP variation sites

4 個 文 庫DS1_T、DS1_CK、DT1_T、DT1_CK 檢測到含SNP 變異位點基因數目分別為6 934、6 427、4 962 和5 141 個。將基于轉錄組中KEGG 顯著富集通路上的差異基因與SNP 變異位點數據進行整合分析，篩選在抗旱關鍵候選基因中含SNP 變異位點的分子標記。結果（表3）顯示，亞油酸代謝通路中CsaV3_4G023820、CsaV3_2G006420，以及戊糖、葡萄糖醛酸轉換通路中CsaV3_2G025090 均檢測到SNP 變異位點，說明該基因和SNP 分子標記可能為旱黃瓜抗旱關鍵基因。

表3 SNP 變異位點與差異基因整合分析Table 3 Integration analysis of SNP variation sites and differential genes

2.7 RT?qPCR分析

通過SNP 變異位點與差異基因整合分析篩選到CsaV3_4G023820、CsaV3_2G006420 和CsaV3_2G02 5090 共計3 個候選基因。利用Beacon Designer 7.0軟件設計出RT?qPCR 引物（表4）。RT?qPCR 分析發現，發芽前期CsaV3_4G023820、CsaV3_2G006420在抗旱型材料PEG 模擬干旱誘導較對照上調表達，而敏感型材料PEG 模擬干旱誘導較對照下調表達；CsaV3_2G025090 在敏感型材料PEG 模擬干旱誘導較對照上調表達，而在PEG 模擬干旱誘導下抗旱型材料中未檢測到相應表達量（圖7）。結果說明，抗旱型材料和敏感型材料干旱響應基因表達模式不同，并且CsaV3_4G023820 和CsaV3_2G006420 可作為響應干旱的候選基因。

表4 SNP 變異位點與差異基因整合分析Table 4 Integration analysis of SNP variation sites and differential genes

圖7 候選基因RT-qPCR 分析Fig. 7 RT-qPCR analysis of candidate genes

3 討論

高通量測序技術為園藝植物重要性狀基因挖掘和分子標記開發提供了強大技術支持。近十年，基因組、轉錄組、蛋白組、代謝組等組學技術迅猛發展，大量相關研究和報道應運而生。如何從海量的組學數據中，整合或篩選有效的生物基因信息，將成為高效且低成本挖掘重要性狀基因、開發分子標記和解析植物復雜生物學機制的一種有效手段。

轉錄組測序可以從mRNA 水平全面地研究植物基因功能及特定的生物學過程，利用該技術能研究植物在干旱下所有基因的表達水平情況，有助于揭示植物對干旱的應答機理［21］。張蕾等［22］以甘藍型油菜‘湘油15’為試材，利用轉錄組分析干旱脅迫下差異基因和關鍵代謝通路。張婷茹等［23］利用轉錄組測序技術分析了角果堿蓬干旱脅迫下差異表達基因，鑒定到426 個上調表達基因和294 個下調表達基因。盧坤等［24］利用轉錄組測序技術鑒定到甘藍型油菜葉片干旱脅迫響應相關基因，從轉錄水平解析了油菜適應干旱環境的分子機制。王玉斌等［25］以高粱干旱敏感材料和耐旱材料為試材，利用轉錄組技術鑒定到耐旱材料4 032 個差異基因，敏感材料8 928 個差異基因，并且推測苯丙素生物合成是高粱耐旱性的主要KEGG 通路。本研究也采用以上轉錄組測序技術方法，以抗旱型和敏感型旱黃瓜為試材，鑒定到PEG 模擬干旱誘導下抗旱型旱黃瓜中1 072 個差異表達基因，敏感型旱黃瓜中1 337 個差異表達基因，并挖掘到亞油酸代謝（csv00591）、戊糖、葡萄糖醛酸轉換（csv00040）、苯丙素的生物合成（csv00940）和植物激素信號轉導（csv04075）等4 個KEGG 顯著富集通路。其中苯丙素的生物合成（csv00940）通路是在很多物種中廣泛被鑒定到為干旱響應通路，這與高粱研究結果相一致［25］。

基于轉錄組數據開發旱黃瓜SNP 分子標記，可為旱黃瓜建立一種高效和低成本的分子標記開發方法。本研究以8 個轉錄組數據庫為基礎開發SNP 分子標記，結合4 個KEGG 富集通路中20 個關鍵候選基因定位信息，篩選出9 個旱黃瓜抗旱性狀相關SNP 分子標記。值得說明的是，亞油酸代謝通路中CsaV3_4G023820、CsaV3_2G006420 在抗旱型和敏感型兩類材料中，檢測到變異位點信息完全一致，均為“C”突變成“G”和“C”突變成“T”，而CsaV3_2G025090 檢測到變異位點信息在抗旱型和敏感型兩類材料中不完全一致，進一步暗示著亞油酸代謝通路及CsaV3_4G023820、CsaV3_2G006420 在旱黃瓜抗旱性方面起重要調控作用，為旱黃瓜分子設計育種提供理論支持。

4 結論

本研究利用轉錄組測序鑒定到富集于亞油酸代謝（csv00591）通路、戊糖、葡萄糖醛酸轉換（csv00040）通路、苯丙素的生物合成（csv00940）通路以及植物激素信號轉導（csv04075）通路上20 個旱黃瓜抗旱相關基因。開發了全轉錄組水平的SNP 分子標記，篩選出抗旱基因CsaV3_4G023820、CsaV3_2G006420及其2 個相關SNP 分子標記。