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

（河北科技師范學(xué)院園藝科技學(xué)院 河北省高校特色園藝植物生物育種應(yīng)用技術(shù)研發(fā)中心，秦皇島 066004）

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

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

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

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

1 材料與方法

1.1 材料

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

1.2 方法

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

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

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

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

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

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

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

2 結(jié)果

2.1 PEG模擬干旱誘導(dǎo)下黃瓜種子萌發(fā)特征分析

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

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

2.2 測序數(shù)據(jù)質(zhì)量分析

從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 樣本測序數(shù)據(jù)質(zhì)量匯總Table 1 Quality summary of samples sequencing data

2.3 差異表達(dá)基因鑒定

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

圖2 比較組合差異基因數(shù)目統(tǒng)計(jì)柱狀圖Fig. 2 Statistical bar chart for comparison of the number of combinatorial differential genes

2.4 差異基因韋恩圖分析

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

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

2.5 差異基因富集分析

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

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

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

圖5 KEGG 富集分析散點(diǎn)圖Fig. 5 Point diagram of KEGG enrichment analysis

將發(fā)芽前期的抗旱型與敏感型2 個材料，PEG模擬干旱誘導(dǎo)與對照的差異表達(dá)基因進(jìn)行關(guān)聯(lián)分析，篩選到亞油酸代謝通路中3 個，分別是CsaV3_4G 023820、CsaV3_2G006420 和CsaV3_4G023810；戊糖、葡萄糖醛酸轉(zhuǎn)換（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；以及植物激素信號轉(zhuǎn)導(dǎo)（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 基于轉(zhuǎn)錄組SNP分子標(biāo)記開發(fā)與篩選

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

圖6 SNP 變異位點(diǎn)統(tǒng)計(jì)分析Fig. 6 Statistical analysis of SNP variation sites

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

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

2.7 RT?qPCR分析

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

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

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

3 討論

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

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

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

4 結(jié)論

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