李東興 錢家連 許慧慧 張國偉 任俊杰 王茜 許洋 王利兵 于海燕 李迎超

(中國林業(yè)科學研究院林業(yè)研究所，北京，100091) (河北省洪崖山國有林場)(河北邢臺市信都區(qū)農業(yè)農村局) (中國林業(yè)科學研究院林業(yè)研究所)

栓皮櫟(Quercusvariabilis)為殼斗科(Fagaceae)櫟屬(Quercus)落葉喬木，耐貧瘠和干旱，是營造防護林、水源涵養(yǎng)林和防風林的優(yōu)良樹種，具有較高的經(jīng)濟和生態(tài)價值[1-2]。栓皮櫟樹皮也稱軟木，作為一種可再生的生物質材料，被廣泛用于葡萄酒瓶塞和保溫隔熱等諸多領域[3]；栓皮櫟種子中營養(yǎng)豐富，淀粉含量高，可用于釀酒和提取漿紗，也是生產燃料乙醇的“非糧”原料[2]。同時，栓皮櫟殼斗也是制備活性炭和糖醛等的優(yōu)質原料。

然而，絕大多數(shù)櫟屬植物無性繁殖困難[4-5]，其苗木繁育主要依靠播種[6]，種子質量決定了成活率及苗木質量。櫟屬植物種子成熟后期都不經(jīng)歷成熟干燥后獲得脫水耐性的過程，這導致它們對脫水敏感[7]。栓皮櫟種子是典型的脫水敏感頑拗性種子，其成熟脫落時含水量高，干燥后在低溫(-18℃左右)環(huán)境中貯藏的常規(guī)方法會使種子迅速失去活力，進而死亡[8-9]。頑拗性種子脫水敏感性是植物重要的功能特征和由多基因共同響應的次級進化[10-11]。種子脫水過程中，細胞感受到水分流失信號，信號經(jīng)傳導后，在轉錄和翻譯等不同水平上做出響應，以調控相關基因表達和蛋白合成，應對脫水造成的傷害。頑拗性種子脫水過程中，與轉錄因子(MYB1R1、ERF106、ABI3、RL40和TGA10)、激素信號傳導(PP2C、JAZ1、JAZ2和JAZ3)和抗氧化酶(SUPEROXIDE DISMUTASE和CATALASE)相關的基因差異表達，它們可能在種子脫水響應中發(fā)揮重要調控作用，參與種子脫水敏感性的分子調控[13-15]。不同頑拗性種子脫水敏感性不盡相同[16-17]，其調控機制也不相同。雖然已有研究對栓皮櫟頑拗性種子脫水敏感性進行探究，指出種子最佳貯藏溫度為(0～2)℃[8]，但是種子脫水敏感性分子調控機制在很大程度上仍然是未知的。

因此，本研究以從同一栓皮櫟優(yōu)良單株上收集的種子為試驗材料，進行脫水和萌發(fā)試驗，然后對脫水敏感關鍵時期的種子進行轉錄組測序，對獲得差異表達基因的功能和其參與的代謝通路進行分析，篩選種子脫水敏感性的候選基因，分析其可能的調控作用。以期為栓皮櫟種子脫水敏感性分子調控機制的研究提供參考，也為其苗木繁育、種子貯藏與種質資源長期保存提供一定理論依據(jù)。

1 材料與方法

供試種子來源于陜西省周至縣樓觀臺林場中的同一栓皮櫟優(yōu)良單株，2019年9月中旬獲得種子后，首先通過目視法手工去除雜質和肉眼可見的劣質種子，然后用55 ℃溫水浸泡水選，5 min后所有仍然漂浮在水面上的種子被丟棄[18]，剩下的種子經(jīng)瀝水且用濾紙擦干表面水分后，平鋪在陰涼處陰干(12～15 h)。在此期間，每隔2～3 h翻動1次種子。陰干完成后，挑選健康、均勻且沒有萌發(fā)的種子，在國家林業(yè)和草原局北方林木種子檢驗中心內進行試驗。

1.1 脫水試驗設計

參照程繼銘等[19]的方法，采用硅膠快速脫水法進行脫水試驗。脫水處理梯度為1、2、3、4、5、6、7、8、9、10、11、12、13、14、15 d(T1～T15)，共15個處理，以不經(jīng)脫水處理的種子為對照(CK)。每個處理均隨機挑選45粒種子，重復3次，共135粒種子，分別置于裝有硅膠的大號自封袋中(m(種子)∶m(硅膠)=1∶3)，混合均勻后，在室溫約25 ℃，相對濕度為40%～45%的實驗室中進行脫水試驗。到達預脫水時間后，每個處理分別隨機取100粒種子進行萌發(fā)試驗，30(每組10粒分3組)粒種子經(jīng)液氮速凍后，迅速放置于冰箱里-80 ℃儲存，用于轉錄組測序分析。

1.2 種子含水量測定

參照GB 2772—1999《林木種子檢驗規(guī)程》，采用稍加改進后的低恒溫烘干法測定種子初始含水量(CIM)。首先將樣品盒和蓋烘干至恒質量后，把種子迅速切成片狀放入樣品盒中，然后在103 ℃的烘箱中烘17 h(達到預設溫度后開始計時)。3次重復，每次重復10粒種子，CIM按照以下公式計算：

CIM=((L-N)/(L-P))×100%。

式中：L為樣品盒和蓋以及樣品的烘前質量(g)；N為樣品盒和蓋以及樣品的烘后質量(g)；P為樣品盒和蓋的質量(g)。

不同脫水處理后種子含水量用種子初始含水量減去脫水量計算得到。

1.3 萌發(fā)試驗

選擇直徑為12 cm的培養(yǎng)皿，雙層濾紙加一層紗布做為發(fā)芽床，然后置于光暗周期為30 ℃光照8 h，20 ℃黑暗16 h的光照培養(yǎng)箱中進行萌發(fā)試驗[20]，5次重復，每個重復20粒種子。以胚根露出至少2 mm作為萌發(fā)標志[21]，每天觀察并記錄1次萌發(fā)情況，萌發(fā)持續(xù)時間為28 d[22]。種子萌發(fā)率(PG)按照以下公式計算：

式中：Dt表示在t日時的萌發(fā)數(shù)量；N表示種子總數(shù)。

1.4 種子總RNA提取、轉錄組測序和分析

用Trizol試劑(Invitrogen,CA,USA)提取栓皮櫟種子總RNA后，分別用Agilent 2100 Bioanalyzer和RNA 6000 Nano LabChip Kit(Agilent,CA,USA)檢測所提樣品RNA的質量和純度。樣品合格后在杭州聯(lián)川生物技術股份有限公司Iilumina測序平臺(Illumina HiSeqTM4000)進行高通量測序，得到轉錄組數(shù)據(jù)。然后使用Cutadapt軟件對測序產生的原始數(shù)據(jù)進行預處理，過濾掉不合格的序列后得到有效數(shù)據(jù)，接著用HISAT2將獲得的有效數(shù)據(jù)與參考基因組[23]比對。使用edgeR對String Tie組裝和定量完的基因進行差異分析，獲得差異表達基因(DEGs)。對差異表達基因進行GO和KEGG富集分析，并以P<0.05作為標準，篩選差異表達基因顯著富集到的GO類別和KEGG代謝途徑。

1.5 實時熒光定量PCR

為驗證RNA-seq結果的準確性，選擇4個差異表達基因進行qRT-PCR驗證。使用PrimeScriptTMReagent Kit with gDNA Eraser試劑盒將RNA反轉錄成第一條鏈的cDNA。然后按照制造商的說明，用KAPA SYBR FAST qPCR Master Mix(KapaBiosystems,USA)進行qRT-PCR試驗。使用NCBI上的Primer-BLAST設計引物，引物序列如表1所示，以PP2A[24]為內參基因，采用2-ΔΔCt計算相對表達量[25]。

表1 qRT-PCR引物序列

2 結果與分析

2.1 種子萌發(fā)率變化

隨著脫水時間的延長，栓皮櫟種子含水量和萌發(fā)率呈現(xiàn)出一直降低的趨勢(圖1)。未經(jīng)脫水處理種子的萌發(fā)率及含水量分別為84.00%和34.27%，脫水12 d(T12)后，種子萌發(fā)率降為0，此時種子含水量為17.17%。脫水初期(CK～T2)種子萌發(fā)率下降最快；經(jīng)過一個較慢的下降階段后(T2～T4)，萌發(fā)率又迅速下降(T4～T6)，最后緩慢降至0。此外，Pearson相關性分析結果表明栓皮櫟種子含水量與萌發(fā)率間存在極顯著的正相關(P<0.01,R=0.976)，這說明其對水分流失高度敏感。此外，臨界含水量和致死含水量分別指50%左右和全部種子無法存活時的含水量，能夠作為頑拗性種子脫水敏感性的標準[26-27]。本研究中，當栓皮櫟種子含水量低于28.20%(T2)時，近一半種子失去生活力，即栓皮櫟種子臨界含水量可能在28.20%左右；當種子含水量低于17.79%(T11)時，種子基本全部失去生活力，即種子致死含水量可能在17.79%左右。因此，為探究栓皮櫟種子脫水過程中的基因表達變化，我們對未脫水、脫水2 d和11 d的栓皮櫟種子樣品進行轉錄組測序分析。

豎線表示均值加減標準誤

2.2 轉錄組測序質量分析

T11、T2和CK3個時期栓皮櫟種子轉錄組測序結果表明共得到了475 939 682條原始序列，過濾掉帶接頭(Adaptor)、含有N(N表示無法確定堿基信息)的比例大于5%和低質量(質量值Q≤10的堿基數(shù)占整個序列的20%以上)的讀取片段后，共得到450 383 524條有效序列。它們的GC含量和Q30的值分別是43.50%～44.00%和98.05%～98.28%(表2)。與參考基因組(GCF_002906115.1)比對后，86.73%～88.33%的測序片段能夠比對到參考基因組上，而且51.03%～52.61%的片段能唯一比對到參考基因組一個位置上。這些結果表明了本次測序數(shù)據(jù)質量較好，能夠進行后續(xù)差異表達分析。

表2 轉錄組測序數(shù)據(jù)統(tǒng)計

2.3 差異表達基因的鑒定

對栓皮櫟種子轉錄組測序結果進行差異表達基因分析。最終，以基因差異表達倍數(shù)≥1并且P<0.05為篩選標準，在T2與CK比較組中獲得4 405個差異表達基因，包括2 724個上調和1 681個下調基因；在T11與T2比較組里鑒定出1 765個上調和2 301個下調，共4 066個差異表達基因；有3 208個上調和2 699個下調，共5 907個差異表達基因，在比較組T11與CK中被獲得(表3)。這些差異表達基因可能在栓皮櫟頑拗性種子脫水過程中發(fā)揮重要作用，與其脫水敏感性密切相關。

表3 不同比較組差異表達基因數(shù)

2.4 差異表達基因的GO和KEGG分析

基因本體論(GO)分析能夠全面描述生物體中基因屬性，其總共有3個主要類別，分別描述基因的分子功能、細胞組分和生物學過程。如圖2所示，在3個比較組中，上調和下調的差異表達基因均主要富集在“防御反應”、“轉錄調控，DNA模板”和“氧化還原過程”等生物學過程類別；“細胞核”、“細胞質”和“質膜”等細胞組分類別以及“蛋白結合”、“ATP結合”和“DNA結合轉錄因子活性”等分子功能類別中。

圖2 差異表達基因的GO富集分析

生物體內，不同基因相互協(xié)調行使其生物學功能，基于途經(jīng)的分析有助于更進一步了解基因的生物學功能。本研究對獲得的差異表達基因進行了KEGG富集分析，結果發(fā)現(xiàn)在比較組T2與CK中，上調的差異表達基因顯著富集在“蛋白質在內質網(wǎng)中的加工”和“植物激素信號傳導”等途徑；“異黃酮生物合成”、“谷胱甘肽生物合成”和“植物激素信號傳導”等是T11與T2中顯著富集的通路；“谷胱甘肽生物合成”和“植物激素信號傳導”等代謝途徑在T11與CK中顯著富集(圖3)。比較組T2與CK中，下調的差異表達基因在“植物激素信號傳導”和“MAPK信號級聯(lián)”等途徑中顯著富集；“其它多糖降解”和“檸檬烯和蒎烯降解”等途徑是T11與T2中顯著富集的通路；T11與CK中，“植物激素信號傳導”和“淀粉和蔗糖代謝”等是差異表達基因顯著富集的通路(圖3)。這些代謝通路有助于了解栓皮櫟種子應對脫水時的代謝信息，以便更好地探究其脫水敏感性的潛在調控機制。

2.5 轉錄因子對脫水的響應

本研究中，我們在比較組T2與CK、T11與T2和T11與CK分別發(fā)現(xiàn)117、110和183個差異表達的轉錄因子基因(表4)。其中，MYB108(LOC112035993)、MYB1R1(LOC111996955)、WRKY28(LOC111998523)、WRKY72(LOC111987575)和ERF1B(LOC111996147)等14個轉錄因子表達持續(xù)上調；UNE10(LOC112020213)和CRF4(LOC112020213)等5個轉錄因子表達持續(xù)下調；ABI3(LOC112020213)、NGA1(LOC112036666)和RAP2-11(LOC112003881)等10個轉錄因子表達先上調后下調；轉錄因子MYB102(LOC111989495)和SRM1(LOC112023214)表達先下調后上調(表5、圖4)。

表4 與轉錄因子和Ca2+途徑相關的差異表達基因

2.6 Ca2+途徑對脫水的響應

本研究中，我們在比較組T2與CK、T11與T2和T11與CK分別發(fā)現(xiàn)12、4和13個與Ca2+途徑相關的差異表達基因(表4)。其中，有3個基因在3個比較組中均差異表達。鈣調蛋白11(LOC111999077)表達持續(xù)上調；鈣調蛋白3(LOC112017550)表達先上調后下調；鈣調蛋白(LOC111983139)表達先下調后上調(圖5)。

圖3 差異表達基因的KEGG富集分析

表5 與轉錄因子相關的差異表達基因

圖4 與轉錄因子相關的差異表達基因熱圖

圖5 與Ca2+途徑相關的差異表達基因熱圖

2.7 與植物激素信號傳導相關的DEGs

脫水2 d后，我們發(fā)現(xiàn)與茉莉酸(JA)信號傳導相關，編碼茉莉酸—氨基酸合成酶的1個JAR1、茉莉酮酸酯ZIM結構域蛋白的2個JAZ、髓細胞組織增生蛋白(MYCs)的4個基因表達均上調；在比較組T11與T2中，2個JAZ基因表達上調，1個JAZ和5個MYC2基因表達下調(圖6A；圖7A和B)。此外，我們同樣發(fā)現(xiàn)脫水2 d后，與水楊酸(SA)信號傳導相關，編碼TGACG基序結合因子的4個基因表達上調，1個基因表達下調；在比較組T11與T2中，1個TGA基因表達下調，1個編碼病程相關蛋白1的基因表達上調(圖6B；圖7A和B)。比較組T11與CK中，與JA和SA信號傳導相關的DEGs最多(圖7C)。

A顯示了茉莉酸信號傳導的途徑；B顯示了水楊酸信號傳導的途徑。圓圈代表代謝產物，方框代表調節(jié)基因。紅色表示基因表達上調，藍色表示下調，橘色表示該基因家族中既有上調基因又有下調基因。

A代表T2與CK;B代表T11與T2;C代表T11與CK。

2.8 qRT-PCR驗證

為驗證RNA-Seq數(shù)據(jù)的準確性，我們對4個DEGs進行了qRT-PCR驗證。盡管表達量并非完全一樣，但是它們的表達趨勢基本一致(表6)。這說明本研究RNA-Seq結果是可靠的，能夠進行后續(xù)的差異表達分析。

表6 對4個DEGs的qRT-PCR驗證

3 結論與討論

頑拗性種子脫水后，測定含水量和萌發(fā)率是研究種子頑拗性的常用方法[19]。與我們的預期一樣，栓皮櫟種子萌發(fā)率隨著含水量的下降而下降，并且含水量與萌發(fā)率之間呈現(xiàn)極顯著的正相關(P<0.01,R=0.976)，相同的結果在其它頑拗性種子的研究中也被得到證實[12-13，20]。本研究中，栓皮櫟種子含水量由34.27%(CK)下降到28.20%(T2)時，近一半的種子已喪失生命力，此時種子已受到嚴重的傷害。此外，當栓皮櫟種子含水量從34.27%(CK)下降到17.17%(T11)時，種子幾乎完全失去生命力，這與Yu et al.[28]對頑拗性七葉樹(Aesculuschinensis)種子的研究類似。

植物中參與調控脅迫誘導型基因表達的各類轉錄因子，尤其是MYB、WRKY、AP2-EREBP和bHLH等轉錄因子家族，是信號傳遞網(wǎng)絡中的關鍵調控因子，在植物對非生物脅迫的響應中發(fā)揮重要作用[29-30]。山茶(Camelliasinensis)種子脫水過程中，差異表達的MYB、WRKY和bHLH等轉錄因子家族基因可能參與種子脫水敏感性調控[12]。本研究中，4個MYB、5個WRKY、3個ERF和2個HSF類轉錄因子在栓皮櫟種子脫水過程中持續(xù)上調。轉錄因子MYB108是ABA誘導細胞死亡的負調控因子，能夠調節(jié)活性氧(ROS)的產生，還可能在轉錄水平上參與JA的信號傳導，調控多種脅迫反應的信號傳導[31-32]。頑拗性種子脫水過程中，ROS的積累是其不耐脫水的原因之一[26]。栓皮櫟種子脫水過程中，持續(xù)上調的MYB108可能增強了ROS的產生，使種子受到脫水傷害。與頑拗性桑寄生(Taxillusichinensis)種子[13]不同，栓皮櫟種子脫水過程中MYB1R1表達持續(xù)上調，這可能是因為不同頑拗性種子對脫水響應不同。轉錄因子ERF1B能夠調控擬南芥(Arabidopsisthaliana)中JA響應基因，并作用于JA信號傳導途徑下游[33]，參與植物對脅迫的響應。栓皮櫟種子脫水過程中，持續(xù)上調的ERF1B轉錄因子可能通過調控激素信號傳導中相關基因的表達，在種子對脫水的響應過程中發(fā)揮作用。此外，我們同樣發(fā)現(xiàn)2個ERF類轉錄因子在栓皮櫟種子脫水過程中持續(xù)下調。CRF轉錄因子能夠調控脅迫信號傳導途徑中的基因表達[34]，持續(xù)下調的CRF4轉錄因子可能干擾了栓皮櫟種子中脅迫信號的傳遞，細胞中防御機制無法啟動，受損程度不斷加深，進而導致種子對脫水敏感。bHLH類轉錄因子能增強細胞清除ROS的能力，降低氧化應激對植物造成的傷害，提高對非生物脅迫耐受性[35]。擬南芥中過表達的AtUNE12基因使其能快速清除ROS，減少脅迫導致的細胞損傷[36]。我們發(fā)現(xiàn)UNE10、HEC1-like和bHLH94-like，3個bHLH轉錄因子基因在栓皮櫟種子脫水過程中持續(xù)下調，推測這可能使種子清除ROS的能力下降，過量積累的ROS使細胞損傷程度不斷加深，進而導致種子對脫水敏感。

植物對各種脅迫因子的響應中，Ca2+作為關鍵的第二信使，傳遞脅迫信號，引發(fā)細胞的防御反應，Ca2+結合蛋白(CMLs)對Ca2+信號傳遞至關重要[37]。CML37能通過影響茉莉酸—異亮氨酸綴合物的合成，進而將Ca2+和JA信號連接起來，CML37基因突變體中，JAR1表達降低，JA信號傳導途徑受到干擾[38]。本研究中，CML11在種子脫水過程中持續(xù)上調，這可能影響了Ca2+和JA信號傳導，導致種子不耐脫水。

植物感受到環(huán)境脅迫信號后，體內會大量合成JA。JA在腺苷酸形成酶的催化下形成高生物活性的JA-Ile，它能與JA受體COIl特異性結合為COIl-JA-Ile復合物，然后通過招募JAZ形成COI1-JA-Ile-JAZ三元復合體，解除對下游MYC2等轉錄因子的抑制作用，激活JA響應基因的表達[39-40]。本研究中，脫水2 d后，JAR1,JAZ和MYC2基因表達上調，這與Luo et al.[41]對棉花(Gossypiumspp.)的研究有些類似。JAR1基因的上調可能促進了COI1-JA-Ile-JAZ復合物的形成，解除JAZ對MYC2的抑制作用，上調的MYC2通過調控下游基因的表達響應種子脫水。擬南芥脫水過程中。MYC2基因的顯著下調可能是JA信號傳導途徑中，下游基因不表達的原因[42]。當種子繼續(xù)脫水時，編碼JAZ的基因表達失衡，MYC2基因表達下調。這可能導致某些下游脅迫響應基因在栓皮櫟種子繼續(xù)脫水時不表達，使種子受到嚴重傷害。JA和SA能相互串擾，共同調控多種信號網(wǎng)絡，進而調控植物防御反應[43]。植物感受到脅迫信號時，SA能和NON-EXPRESSER OF PATHOGENESIS-RELATED GENES 1(NPR1)結合[44]，影響NPR1結構域變化，誘導下游基因轉錄。NPR1與bZIP類轉錄因子TGA家族和熱休克轉錄因子HSFA1間均存在互作關系，進而在植物對SA信號的響應中發(fā)揮重要調控作用[45-46]。栓皮櫟種子脫水過程中，水分流失誘導SA信號通路中TGA相關基因表達失衡，既有上調也有下調，這表明TGA對栓皮櫟種子響應脫水中SA信號的調控作用不同。擬南芥II類TGA，即tga2、tga5和tga6突變體中，SA下游的PR1基因表達量比野生型高，其負調控SA下游基因，但該突變體中SA誘導PR1表達能力喪失，植物不能產生SAR[46]。當水分繼續(xù)流失時，TGA基因表達下調，PR-1基因表達上調，其可能通過調控下游脅迫相關基因，在栓皮櫟種子應對嚴重水分流失過程中發(fā)揮作用。因此，我們推測脫水誘導栓皮櫟頑拗性種子中JA和SA信號傳導途徑的變化是其脫水敏感性的潛在調控機制。

綜上所述，栓皮櫟頑拗性種子脫水過程中，MYB108、MYB1R1、ERF1B、UNE10、CRF4、CML11、JAZ、MYC2和TGA等轉錄因子及Ca2+、JA和SA信號相關的基因差異表達。它們可能參與栓皮櫟種子脫水過程中ROS的產生和清除，以及脅迫信號的傳導過程，進而在種子脫水敏感性的調控中發(fā)揮作用。本研究為栓皮櫟種子脫水敏感性分子調控機制的研究提供理論依據(jù)，有助于種子貯藏和種質資源的長期保存，同時也能為其它頑拗性種子的研究提供一定參考。