李鵬程 畢真真 孫 超 秦天元 梁文君 王一好 許德蓉 劉玉匯 張俊蓮 白江平,*

DNA甲基化參與調(diào)控馬鈴薯響應(yīng)干旱脅迫的關(guān)鍵基因挖掘

李鵬程1,2,**畢真真1,2,**孫 超1,2秦天元1,2梁文君1,2王一好1,2許德蓉1,2劉玉匯1,2張俊蓮1,2白江平1,2,*

1甘肅省干旱生境作物學(xué)重點(diǎn)實(shí)驗(yàn)室/ 甘肅省作物遺傳改良與種質(zhì)創(chuàng)新重點(diǎn)實(shí)驗(yàn)室, 甘肅蘭州 730070;2甘肅農(nóng)業(yè)大學(xué)農(nóng)學(xué)院, 甘肅蘭州 730070

植物受到干旱脅迫時(shí), 會(huì)通過(guò)DNA甲基化做出快速反應(yīng)以幫助其應(yīng)對(duì)脅迫。為探究在干旱脅迫下, DNA甲基化是如何影響基因轉(zhuǎn)錄表達(dá), 本研究對(duì)甘露醇模擬干旱和5-azadC (去甲基化)處理下, 抗旱性不同的2個(gè)馬鈴薯品種(抗旱型, 青薯9號(hào); 干旱敏感型, 大西洋)進(jìn)行轉(zhuǎn)錄組學(xué)分析, 以Fold-change > 2和校正后< 0.01進(jìn)行差異表達(dá)基因(DEG)的篩選。GO富集分析發(fā)現(xiàn), 2種處理都共同顯著富集到氧化應(yīng)激和碳水化合物代謝過(guò)程相關(guān)的GO term。說(shuō)明不同耐旱性馬鈴薯在響應(yīng)干旱脅迫時(shí), 與這些GO term相關(guān)的基因也受DNA去甲基化調(diào)控。對(duì)既響應(yīng)干旱又響應(yīng)DNA去甲基化的1345個(gè)DEG進(jìn)行KEGG功能富集發(fā)現(xiàn), 與植物抗旱相關(guān)的通路有植物MAPK信號(hào)途徑、植物激素信號(hào)轉(zhuǎn)導(dǎo)途徑、植物谷胱甘肽代謝通路、糖酵解與糖異生和磷酸肌醇代謝通路。說(shuō)明這些通路相關(guān)基因在大西洋和青薯9號(hào)2個(gè)抗旱性不同的馬鈴薯品種中, 響應(yīng)干旱的敏感性受DNA甲基化調(diào)控。接著對(duì)DEG上游1500 bp啟動(dòng)子區(qū)域進(jìn)行順式作用原件和甲基化CpG島分析發(fā)現(xiàn), 干旱脅迫下參與植物谷胱甘肽代謝的GST基因通過(guò)DNA去甲基化來(lái)降低啟動(dòng)子區(qū)ABRE和CAAT-box作用元件的甲基化水平, 進(jìn)而激活該基因的表達(dá)以應(yīng)對(duì)干旱脅迫。因此, 利用比較轉(zhuǎn)錄組學(xué)分析干旱和DNA去甲基化處理下的差異基因, 可挖掘到DNA甲基化參與調(diào)控馬鈴薯響應(yīng)干旱脅迫的相關(guān)基因, 為研究馬鈴薯干旱脅迫響應(yīng)的表觀遺傳學(xué)機(jī)理提供新的研究思路。

馬鈴薯; 干旱脅迫; DNA甲基化; 轉(zhuǎn)錄組; 谷胱甘肽S轉(zhuǎn)移酶

干旱脅迫通常發(fā)生在土壤水分不足以維持植物正常的生命活動(dòng)且蒸發(fā)量大的區(qū)域[1]。水分缺乏是限制植物生長(zhǎng)發(fā)育的關(guān)鍵因素, 影響植物的伸長(zhǎng)和膨大生長(zhǎng)[2-3]以及土壤養(yǎng)分供應(yīng)的不平衡等[4]。馬鈴薯作為四大主糧作物之一, 也是水分敏感作物。中國(guó)的馬鈴薯主要種植在西北干旱和半干旱地區(qū), 因此干旱嚴(yán)重影響了馬鈴薯的植株生長(zhǎng)、塊莖大小和商品屬性[5]。

近年來(lái), 除了闡明響應(yīng)干旱脅迫的信號(hào)轉(zhuǎn)導(dǎo)機(jī)制外, 越來(lái)越多的研究表明, 表觀遺傳調(diào)控在植物干旱脅迫響應(yīng)中具有重要的作用[6-7]。DNA甲基化是發(fā)現(xiàn)最早、研究最深入的表觀遺傳修飾中的一種。DNA甲基化的修飾往往與基因的轉(zhuǎn)錄相關(guān), 即DNA甲基化程度越高, 基因的轉(zhuǎn)錄水平越低[8-9]。干旱脅迫會(huì)引起植物體內(nèi)DNA甲基化水平的改變, 激活脅迫相關(guān)基因的表達(dá)、轉(zhuǎn)座子的活性以及調(diào)控啟動(dòng)子區(qū)域的開(kāi)關(guān), 以增強(qiáng)植物適應(yīng)及抵御脅迫的能力[10]。Fei等[11]研究表明, 2個(gè)不同耐旱性的小麥品種中的胞質(zhì)甘油醛-3-磷酸脫氫酶(GAPC)[12]啟動(dòng)子中不同類(lèi)型和水平的DNA甲基化動(dòng)態(tài)變化與干旱脅迫反應(yīng)中基因表達(dá)之間存在顯著關(guān)系。Liang等[13]研究表明, 轉(zhuǎn)錄因子可能通過(guò)DNA甲基化調(diào)控從而在毛果楊干旱脅迫反應(yīng)中發(fā)揮重要作用。對(duì)耐旱(Bachar)和干旱敏感(F177)蠶豆基因型的全基因組DNA甲基化多態(tài)性進(jìn)行分析, 確定了6個(gè)與干旱脅迫相關(guān)的差異甲基化區(qū)域, 且發(fā)現(xiàn)干旱脅迫下Bachar的抗旱相關(guān)基因的表達(dá)水平高于F177。這表明全基因組DNA甲基化變化可能是蠶豆對(duì)響應(yīng)干旱脅迫的重要調(diào)控機(jī)制[14]。本試驗(yàn)前期通過(guò)5-azadC處理不同干旱敏感的馬鈴薯品種發(fā)現(xiàn), 干旱脅迫下DNA去甲基化參與調(diào)控馬鈴薯試管苗表型性狀及生理生化的響應(yīng)[15]。而有關(guān)DNA甲基化如何影響干旱脅迫下馬鈴薯相關(guān)基因的變化尚不清楚。

通過(guò)轉(zhuǎn)錄組測(cè)序比較不同物種中不同基因型的轉(zhuǎn)錄表達(dá)是探索非生物脅迫下抗性基因、闡明代謝途徑以及解析響應(yīng)脅迫的耐受性機(jī)制的最合適技術(shù)之一[16-17]。因此, 本研究通過(guò)比較轉(zhuǎn)錄組學(xué)分析以探討耐旱型品種青薯9號(hào)和干旱敏感型品種大西洋在干旱脅迫和DNA甲基化抑制劑處理下相關(guān)基因的表達(dá)特性, 以期獲得受DNA甲基化調(diào)控響應(yīng)干旱脅迫的相關(guān)基因, 從而為更好地闡明DNA甲基化參與調(diào)控馬鈴薯響應(yīng)干旱脅迫的分子機(jī)制, 擴(kuò)展對(duì)馬鈴薯耐旱機(jī)制的認(rèn)識(shí), 為馬鈴薯抗旱分子育種的新方法提供一定的理論依據(jù)。

1 材料與方法

1.1 試驗(yàn)材料及處理方法

馬鈴薯耐旱品種為青薯9號(hào)(Qingshu 9), 干旱敏感品種為大西洋(Atlantic), 均由甘肅省作物遺傳改良與種質(zhì)創(chuàng)新重點(diǎn)實(shí)驗(yàn)室提供。使用紙橋法將2個(gè)品種的馬鈴薯組培苗莖段接種到MS液體培養(yǎng)基(不含瓊脂), 每個(gè)莖段帶有1個(gè)腋芽, 生長(zhǎng)環(huán)境條件為溫度(22±2)℃, 光強(qiáng)25.0~37.5 μmol m-2s-1, 光照16 h d-1。培養(yǎng)24 d后, 用200 mmol L-1甘露醇和60 μmol L-1DNA甲基化抑制劑(5-azadC)處理馬鈴薯試管苗, 分別取處理2、6、12和24 h以及未進(jìn)行處理(對(duì)照, 0 h)試管苗的莖葉, 迅速使用液氮冷凍, 存儲(chǔ)在-80℃超低溫冰箱備用。每個(gè)處理3個(gè)生物學(xué)重復(fù)。5-azadC (5-aza-2’-deoxycytidine)購(gòu)于Sigma公司。

1.2 莖葉總RNA提取、文庫(kù)構(gòu)建及轉(zhuǎn)錄組測(cè)序

使用植物總RNA提取TRIzol試劑(Invitgen, 美國(guó))分別提取54個(gè)樣品{2個(gè)品種′[對(duì)照+(2種處理′4個(gè)時(shí)間點(diǎn))]′3次生物學(xué)重復(fù)}的總RNA, 用Turbo DNase I (Ambion, 美國(guó))處理30 min后使用RNeasy Plant Mini Kit (QIAGEN, 德國(guó))進(jìn)行純化; 然后使用TruSeq RNA sample Prep V2 kit (Illumina, 美國(guó))構(gòu)建轉(zhuǎn)錄組測(cè)序cDNA文庫(kù), 并利用Agilent 2200 TapeStation (Agilent, 美國(guó))系統(tǒng)檢測(cè)cDNA文庫(kù)質(zhì)量; 最后使用Illumina HiSeqX10測(cè)序儀上進(jìn)行雙端測(cè)序。

1.3 轉(zhuǎn)錄組數(shù)據(jù)分析

首先使用FASTQC對(duì)下機(jī)數(shù)據(jù)進(jìn)行質(zhì)控并計(jì)算Q20、Q30和GC含量。然后從ENSEMBL植物數(shù)據(jù)庫(kù)中下載了馬鈴薯參考基因組序列(SolTub 3.0)和注釋文件(gtf), 利用HISAT2將測(cè)序reads比對(duì)到馬鈴薯參考基因組中, 獲得sam文件; 利用Samtools 工具將sam文件轉(zhuǎn)換為bam文件并排序, 使用Cufflinks對(duì)排序后的bam文件進(jìn)行轉(zhuǎn)錄本的組裝和定量; 使用HTseq-count計(jì)算各樣本轉(zhuǎn)錄本的counts數(shù), 利用Python腳本計(jì)算參考基因組各基因的長(zhǎng)度, 使用R包Dseq2計(jì)算各基因的TPM (Transcripts Per Million, 每百萬(wàn)條reads的轉(zhuǎn)錄本)和樣本間的差異表達(dá)分析。將校正后的<0.01和變化倍數(shù)大于2定義為差異表達(dá)基因(DEG), 通過(guò)比較青薯9號(hào)和大西洋在同一處理同一時(shí)間點(diǎn)的表達(dá)水平來(lái)確定差異表達(dá)基因。采用在線工具g:Profiler[18](https://biit.cs. ut.ee/gprofiler_beta/gost)進(jìn)行GO和KEGG功能富集分析, 以<0.05作為顯著富集GO term和KEGGpathway的篩選條件。

1.4 啟動(dòng)子區(qū)甲基化CpG島預(yù)測(cè)和順式作用元件分析

先使用perl腳本批量提取所分析基因上游1500 bp啟動(dòng)子序列, 然后利用在線軟件MethPrimer[19](http://www.urogene.org/cgi-bin/methprimer/MethPrimer.cgi)對(duì)差異富集基因5￠調(diào)控區(qū)進(jìn)行CpG島預(yù)測(cè), 預(yù)測(cè)參數(shù)設(shè)置為: CpG島長(zhǎng)度>100 bp, CPG檢測(cè)含量/CpG期望含量 (Obs/Exp)>0.60, GC含量>50%。使用Plant Care[20](http://bioinformatics.psb.ug ent. be/webtools/plantcare/html/)預(yù)測(cè)基因啟動(dòng)子的順式作用元件?

1.5 實(shí)時(shí)熒光定量PCR驗(yàn)證

隨機(jī)選擇22個(gè)差異基因, 使用Primer Express 3.0.1設(shè)計(jì)引物(表1), 進(jìn)行實(shí)時(shí)熒光定量PCR分析。采用cDNA合成試劑盒9 (TaKaRa, 日本)將各樣品3次生物學(xué)重復(fù)的RNA反轉(zhuǎn)錄為cDNA, 以為內(nèi)參, 利用QuantStudio5實(shí)時(shí)熒光定量PCR系統(tǒng) (ABI, 美國(guó))進(jìn)行熒光定量檢測(cè)。每個(gè)樣品3次生物學(xué)重復(fù), 按2–ΔΔCT相對(duì)定量法計(jì)算相對(duì)表達(dá)量。

表1 實(shí)時(shí)熒光定量PCR所用的引物

2 結(jié)果與分析

2.1 不同馬鈴薯品種干旱和DNA甲基化處理下差異表達(dá)基因的鑒定與分析

在數(shù)據(jù)過(guò)濾和質(zhì)量評(píng)估后, 54個(gè)樣品產(chǎn)生了大約2547百萬(wàn)的reads數(shù), 平均每個(gè)樣本約產(chǎn)生7.1 G堿基, GC含量均在44.17%~46.83%范圍之間。應(yīng)用RNA-seq分析比對(duì)軟件HISAT2將各樣本所有的reads與馬鈴薯參考基因組序列進(jìn)行比對(duì), 比對(duì)率為75.56%~94.64% (表2)。

為了分析不同耐旱型馬鈴薯品種在甘露醇和5-Azadc處理期間轉(zhuǎn)錄水平的差異, 使用R包Dseq2對(duì)大西洋和青薯9號(hào)在同一時(shí)間樣本的TPM值計(jì)算了所有差異基因的表達(dá)。以Flod-change>2和校正后<0.01作為差異表達(dá)基因 (DEG)的篩選標(biāo)準(zhǔn)。結(jié)果表明, 在對(duì)照 (0 h)下, A vs. Q中篩選到1030個(gè)DEG,其中有515個(gè)上調(diào)和515個(gè)下調(diào)表達(dá)基因。干旱處理下, 在A vs. Q 2 h中篩選到568個(gè)DEG, 其中有384個(gè)上調(diào)表達(dá)基因, 184個(gè)下調(diào)表達(dá)基因, 在A vs. Q 6 h中篩選到994個(gè)DEG, 其中有560個(gè)上調(diào)表達(dá)基因, 434個(gè)下調(diào)表達(dá)基因, 在A vs. Q 12 h中篩選到519個(gè)DEG, 其中有321個(gè)上調(diào)表達(dá)基因, 198個(gè)下調(diào)表達(dá)基因, 在A vs. Q 24 h中篩選到1251個(gè)DEG, 其中有471個(gè)上調(diào)表達(dá)基因, 780個(gè)下調(diào)表達(dá)基因。DNA去甲基化處理下, 在A vs. Q 2 h中篩選到1603個(gè)DEG, 其中有1033個(gè)上調(diào)表達(dá)基因, 570個(gè)下調(diào)表達(dá)基因。在A vs. Q 6 h中篩選到501個(gè)DEG, 其中有641個(gè)上調(diào)表達(dá)基因, 434個(gè)下調(diào)表達(dá)基因。在A vs. Q 12 h中篩選到1153個(gè)DEG, 其中有576個(gè)上調(diào)表達(dá)基因, 577個(gè)下調(diào)表達(dá)基因。在A vs. Q 24 h中篩選到981個(gè)DEG 其中有465個(gè)上調(diào)表達(dá)基因, 516個(gè)下調(diào)表達(dá)基因(表3)。

為了進(jìn)一步分析馬鈴薯不同品種在不同處理下特定和共同差異表達(dá)基因變化。利用R語(yǔ)言分別繪制了干旱和DNA去甲基化處理下的Upset圖(圖1-A, B)。與對(duì)照相比, 干旱處理下, A vs. Q在2、6、12和24 h時(shí)差異表達(dá)基因總數(shù)為1790, 在2、6、12和24 h中特異表達(dá)的差異基因?yàn)?70、432、175和507個(gè), 在4個(gè)時(shí)間點(diǎn)共同表達(dá)的基因有17個(gè); DNA去甲基化處理下, A vs. Q在2、6、12和24 h中差異表達(dá)基因總數(shù)為2617, 在2、6、12和24 h中特異表達(dá)的差異基因?yàn)?90、435、345和291個(gè), 在這4個(gè)時(shí)間點(diǎn)共同表達(dá)的基因有29個(gè)。同時(shí)繪制了對(duì)照、干旱和DNA甲基化抑制劑處理下A vs. Q差異表達(dá)基因的韋恩圖(圖1-C), 干旱和DNA甲基化抑制劑處理下共同表達(dá)的基因有1354個(gè), 表明這1354個(gè)基因可能受DNA甲基化調(diào)控參與馬鈴薯干旱響應(yīng), 可作為后續(xù)基因挖掘的重點(diǎn)部分。

表2 測(cè)序reads數(shù)比對(duì)到參考基因組的基本統(tǒng)計(jì)

C: 對(duì)照處理; A: DNA甲基化抑制劑處理; M: 干旱處理; At: 大西洋; Qs: 青薯9號(hào)。

C: control treatment; A: 5-azadC treatment; M: mannitol treatment; At: Atlantic; Qs: Qingshu 9.

表3 干旱、DNA去甲基化處理下不同品系間差異表達(dá)基因數(shù)量

A和B是UpSet圖, 分別表示干旱和DNA去甲基化處理下AvsQ各時(shí)間點(diǎn)比較得到的集合的重疊和DEG的數(shù)量。C是韋恩圖。

A and B are UpSet diagrams, showing the overlap of sets and the number of DEG obtained by comparison of AvsQ at each time point under drought and DNA demethylation, respectively. C is the venn diagram.

2.2 不同馬鈴薯品種干旱和DNA甲基化處理下差異表達(dá)基因GO分析

甘露醇處理下, 對(duì)A vs. Q在2、6、12和24 h中共1790個(gè)差異表達(dá)基因進(jìn)行GO功能富集分析發(fā)現(xiàn), 這些差異基因主要富集包括生物過(guò)程、分子功能、細(xì)胞組分三大分類(lèi)下的42個(gè)GO term (圖2-A)。生物學(xué)過(guò)程共富集到15個(gè)GO term, 主要是苯丙烷類(lèi)分解代謝過(guò)程(GO:0046274和GO:0009698)、木質(zhì)素代謝過(guò)程(GO:0009808)、氧化還原過(guò)程(GO:0055114)、蛋白質(zhì)磷酸化(GO:0006468)、植物型細(xì)胞壁組織或生物發(fā)生(GO:0071669)和碳水化合物代謝過(guò)程(GO:0005975)。分子功能共富集到22個(gè)GO term, 與催化活性(GO:0003824)、氧化還原酶活性(GO:0016491和GO:0016682)、碳水化合物結(jié)合(GO:0030246)、水解O-糖基化合物的水解酶活性(GO:0004553和GO:0016798)、天冬氨酸肽酶活性(GO:0070001和GO:0004190)、蛋白絲氨酸/蘇氨酸激酶活性(GO:0004674)、木糖基轉(zhuǎn)移酶活性(GO:0016762)、蛋白激酶活性(GO:0004672)、內(nèi)肽酶抑制劑活性(GO:0004866)、鈣調(diào)素結(jié)合(GO: 0005516)和磷酸轉(zhuǎn)移酶活性(GO:0016773)相關(guān)。細(xì)胞組分共富集到5個(gè)GO term, 主要是質(zhì)外體(GO:0005576)、細(xì)胞膜(GO:0005576和GO:0030312)和細(xì)胞壁(GO:0005618)。同樣地, 對(duì)DNA甲基化抑制劑處理下, A vs. Q在2、6、12和24 h共2617個(gè)差異表達(dá)基因進(jìn)行GO功能富集分析(圖2-B)發(fā)現(xiàn), 這些差異基因也顯著富集到生物過(guò)程、分子功能、細(xì)胞組分下的39個(gè)GO term。其中富集到生物學(xué)過(guò)程的共16個(gè), 主要包括氧化還原過(guò)程(GO: 0055114)、碳水化合物代謝過(guò)程(GO:0005975)、糖代謝過(guò)程(GO:0010411、GO:0006073、GO:0044042、GO: 0010410、GO:0044264、GO:0010383和GO:0005976)和蛋白水解(GO:0045861)過(guò)程。分子功能共富集到18個(gè)GO term, 與催化活性(GO:0003824)、加氧酶活性(GO:0004497)、氧化還原酶活性(GO:0016491)、血紅素結(jié)合(GO:0020037)、木糖基轉(zhuǎn)移酶活性(GO: 0016762)、輔酶因子結(jié)合(GO:0048037)、水解O-糖基化合物的水解酶活性(GO:0004553)和內(nèi)肽酶抑制劑活性(GO:0004866)相關(guān)。細(xì)胞組分共富集到5個(gè)GO term, 主要是是質(zhì)外體(GO:0005576)、細(xì)胞壁(GO: 0005618)和細(xì)胞膜(GO:0005576和GO:0030312)。

同時(shí)比較干旱處理和去甲基化處理下所富集的GO term發(fā)現(xiàn), 2種處理都共同顯著富集到催化活性(GO:0003824)、氧化還原酶活性(GO:0016491和GO:0016682)、木糖基轉(zhuǎn)移酶活性(GO:0016762)、內(nèi)肽酶抑制劑活性(GO:0004866)、碳水化合物代謝過(guò)程、質(zhì)外體(GO:0005576)、細(xì)胞壁(GO:0005618)和細(xì)胞膜(GO:0005576和GO:0030312)。

2.3 兩種處理下共同差異表達(dá)基因表達(dá)聚類(lèi)分析

聚類(lèi)分析主要是把表達(dá)模式一致的差異表達(dá)基因歸為一類(lèi)。以干旱和DNA甲基化抑制劑處理下A vs. Q共同差異表達(dá)的1354個(gè)基因的表達(dá)量為目標(biāo), 對(duì)各樣本進(jìn)行表達(dá)層次聚類(lèi)分析。由圖3可知, A vs. Q在干旱和DNA甲基化抑制劑處理下的表達(dá)模式基本一致。表達(dá)趨勢(shì)主要分為: 2種處理下持續(xù)下調(diào)表達(dá)基因, 說(shuō)明與馬鈴薯響應(yīng)干旱脅迫相關(guān)的負(fù)向調(diào)節(jié)基因受到DNA去甲基化的負(fù)調(diào)控; 以及2種處理下持續(xù)上調(diào)表達(dá)基因, 說(shuō)明與馬鈴薯響應(yīng)干旱脅迫相關(guān)的正向調(diào)節(jié)基因受到DNA去甲基化的正調(diào)控。

2.4 兩種處理下共同差異表達(dá)基因KEGG功能富集分析

為了進(jìn)一步分析DNA甲基化參與調(diào)控馬鈴薯干旱脅迫響應(yīng)的相關(guān)基因的生物學(xué)功能, 利用KEGG數(shù)據(jù)庫(kù)對(duì)2個(gè)馬鈴薯品種在甘露醇和DNA甲基化抑制劑處理下共同差異表達(dá)的基因進(jìn)行通路分析發(fā)現(xiàn), 1345個(gè)差異表達(dá)基因顯著富集到12個(gè)KEGG pathways (圖4-A), 其中與植物抗旱相關(guān)的通路有5個(gè), 分別是植物MAPK信號(hào)通路(KO04016, 富集到DEG有11個(gè))、谷胱甘肽代謝(KO00480, 富集到DEG有8個(gè))、糖酵解/糖異生(KO00010, 富集到DEG有8個(gè))、磷酸肌醇代謝(KO00562, 富集到DEG有5個(gè))、植物激素信號(hào)轉(zhuǎn)導(dǎo)通路(KO04075, 富集到DEG有14個(gè))。通過(guò)pathview繪制富集通路的基因在干旱和去甲基化處理下的表達(dá)模式和調(diào)控網(wǎng)絡(luò)發(fā)現(xiàn), 大西洋和青薯9號(hào)在響應(yīng)干旱處理和DNA去甲基化處理下, 谷胱甘肽代謝(圖4-B)和磷酸肌醇代謝(圖4-C)通路中差異基因多為上調(diào)表達(dá)基因, 且2種處理下各時(shí)間點(diǎn)表達(dá)趨勢(shì)類(lèi)似。而MAPK信號(hào)(圖4-D)、糖酵解/糖異生代謝(圖4-E)和植物激素信號(hào)轉(zhuǎn)導(dǎo)(圖4-F)通路中差異基因多為下調(diào)表達(dá)基因。這些基因的功能注釋信息如表4所示。

2.5 兩種處理下共同差異表達(dá)基因啟動(dòng)子分析

對(duì)KEGG富集到與植物抗旱相關(guān)的5個(gè)通路中的39個(gè)基因, 進(jìn)行啟動(dòng)子區(qū)CpG島預(yù)測(cè)和順式作用元件分析。首先, 從馬鈴薯基因組序列中提取39個(gè)基因1500 bp的啟動(dòng)子序列。使用在線軟件MethPrimer預(yù)測(cè)這些基因的啟動(dòng)子區(qū)域的CpG島發(fā)現(xiàn), 只有1個(gè)與植物谷胱甘肽代謝通路相關(guān)的基因 ()含有CpG島, 位于該基因啟動(dòng)子序列的-736~-866之間, 大小為131 bp (圖5)。同時(shí), 使用在線數(shù)據(jù)庫(kù)Plant Care對(duì)該基因啟動(dòng)子序列進(jìn)行順式作用元件分析發(fā)現(xiàn), 在基因的CpG島區(qū)域有ABRE、CAAT-box、A-box元件, 其中ABRE和CAAT-box與植物響應(yīng)干旱脅迫相關(guān)。表明, 馬鈴薯在響應(yīng)干旱脅迫時(shí), 谷胱甘肽轉(zhuǎn)移酶基因啟動(dòng)子區(qū)的ABRE和CAAT-box作用元件通過(guò)DNA去甲基化激活該基因的表達(dá)以應(yīng)對(duì)逆境脅迫。

2.6 qRT-PCR驗(yàn)證

從KEGG顯著富集到的5條干旱相關(guān)通路的基因中隨機(jī)挑選22個(gè)基因進(jìn)行實(shí)時(shí)熒光定量PCR驗(yàn)證。利用反轉(zhuǎn)錄cDNA試劑盒將轉(zhuǎn)錄組測(cè)序時(shí)保存的RNA樣品進(jìn)行反轉(zhuǎn)錄, 使用qRT-PCR檢測(cè)這些基因在不同處理下在大西洋和青薯9號(hào)中的相對(duì)表達(dá)量。結(jié)果表明, qRT-PCR結(jié)果和轉(zhuǎn)錄組測(cè)序結(jié)果變化趨勢(shì)基本一致(圖6)。說(shuō)明轉(zhuǎn)錄組測(cè)序數(shù)據(jù)的可靠性, 并進(jìn)一步驗(yàn)證了這些基因?qū)Ω珊岛虳NA甲基化抑制劑的響應(yīng)。

3 討論

植物對(duì)干旱脅迫反應(yīng)的分子機(jī)制包括轉(zhuǎn)錄前表觀遺傳調(diào)控和轉(zhuǎn)錄后激素、轉(zhuǎn)錄因子、脅迫反應(yīng)基因等信號(hào)分子的參與。研究發(fā)現(xiàn), 植物通過(guò)基因組DNA甲基化和去甲基化的動(dòng)態(tài)變化來(lái)迅速適應(yīng)外部干旱環(huán)境, 表明DNA甲基化調(diào)控在植物響應(yīng)干旱脅迫中起重要作用[21]。本研究利用甘露醇和DNA甲基化抑制劑分別處理2種耐旱性不同的馬鈴薯試管苗, 通過(guò)比較RNA-seq分析2個(gè)品種轉(zhuǎn)錄本在相同處理和不同處理下5個(gè)時(shí)間點(diǎn)的基因表達(dá)變化, 在分子水平探究DNA甲基化如何參與調(diào)控馬鈴薯對(duì)干旱脅迫的響應(yīng)。

植物通過(guò)升高或降低DNA在特定位點(diǎn)發(fā)生特定類(lèi)型的甲基化水平, 來(lái)打開(kāi)或關(guān)閉一系列與逆境響應(yīng)相關(guān)的基因, 進(jìn)而啟動(dòng)脅迫應(yīng)答通路, 提高響應(yīng)逆境相關(guān)酶活性, 從而防御逆境脅迫對(duì)植物功能基因表達(dá)造成的影響, 保護(hù)植物自身的生長(zhǎng)發(fā)育[15]。植物對(duì)干旱脅迫的耐受性是由復(fù)雜的多組分信號(hào)通路觸發(fā), 以恢復(fù)細(xì)胞動(dòng)態(tài)平衡并促進(jìn)植物生存[22]。信號(hào)轉(zhuǎn)導(dǎo)的一般途徑包括由應(yīng)激特異性受體(如G蛋白偶聯(lián)受體、蛋白激酶受體或磷酸肌醇)檢測(cè)和接受信號(hào), 然后由第二信使擴(kuò)增信號(hào), 導(dǎo)致參與調(diào)控響應(yīng)干旱脅迫功能基因的激活或表達(dá)調(diào)節(jié)蛋白活性的激活與抑制[23]。Molina等[24]對(duì)鷹嘴豆根部在干旱處理下的轉(zhuǎn)錄組測(cè)序結(jié)果發(fā)現(xiàn), 蛋白激酶基因和Ca2+應(yīng)答基因在感知水分脅迫信號(hào)中發(fā)揮重要作用。本研究也發(fā)現(xiàn), 馬鈴薯在干旱脅迫下, 蛋白激酶 ()、Ca2+結(jié)合蛋白(和)、肌醇加氧酶 (和)相關(guān)基因響應(yīng)了干旱脅迫, 同時(shí)這些基因也響應(yīng)了DNA去甲基化處理, 且2種處理下的表達(dá)模式基本一致, 說(shuō)明馬鈴薯在響應(yīng)干旱脅迫時(shí)蛋白激酶、Ca2+結(jié)合蛋白、肌醇加氧酶相關(guān)基因的表達(dá)受到DNA甲基化的調(diào)控。

植物激素是一類(lèi)小的內(nèi)源性信號(hào)分子, 如生長(zhǎng)素(IAA)、赤霉素(GA)、細(xì)胞分裂素(CK)、脫落酸(ABA)、乙烯(ET)、茉莉酸(JA)和水楊酸(SA)。植物激素是一種介質(zhì), 不僅作用于內(nèi)源性發(fā)育過(guò)程, 而且還向多種激素反應(yīng)途徑提供環(huán)境刺激, 以適應(yīng)不利條件[25-26]。Jain和Khurana[27]發(fā)現(xiàn), 水稻生長(zhǎng)素反應(yīng)相關(guān)基因在生殖發(fā)育和非生物脅迫過(guò)程中受到不同程度的調(diào)控。本研究也得到了類(lèi)似的結(jié)果, 在“植物激素信號(hào)轉(zhuǎn)導(dǎo)”中富集到14個(gè)基因, 其中5個(gè)是生長(zhǎng)素信號(hào)轉(zhuǎn)導(dǎo)基因, 編碼“生長(zhǎng)素響應(yīng)蛋白”和“GH3生長(zhǎng)素響應(yīng)基因”。因此, 我們推測(cè)生長(zhǎng)素可能在馬鈴薯響應(yīng)干旱脅迫中起重要作用。本研究還鑒定到與ABA、ET和SA調(diào)控相關(guān)的基因。共有9個(gè)差異表達(dá)基因, 分別編碼ERF2轉(zhuǎn)錄因子(PGSC0003 DMG400026261)、PR1蛋白(PGSC0003DMG40000 5108、PGSC0003DMG400005112、PGSC0003DMG 400005116和PGSC0003DMG400005110)和熱休克蛋白(PGSC0003DMG400023921); 同時(shí)PR1蛋白和熱休克蛋白也被富集到植物MAPK信號(hào)轉(zhuǎn)導(dǎo)通路。這些基因在2種處理下的表達(dá)模式基本一致, 這是因?yàn)樵隈R鈴薯響應(yīng)干旱脅迫時(shí), DNA甲基化參與了這些基因的表達(dá)調(diào)控。此外, 本研究還發(fā)現(xiàn)了參與植物糖酵解和糖異生途徑的相關(guān)基因, 如磷酸甘油酸變位酶()、丙酮酸脫羧酶()、己糖激酶(), 在馬鈴薯響應(yīng)干旱脅迫時(shí)受DNA甲基化的調(diào)控。

表4 共同富集KEGG與干旱相關(guān)通路基因功能注釋

(續(xù)表4)

植物谷胱甘肽代謝中有3種重要的酶, 分別是谷胱甘肽還原酶(GR)、谷胱甘肽過(guò)氧化物酶(GP)和谷胱甘肽S轉(zhuǎn)移酶(GST)。這3種酶活性的變化與植物應(yīng)對(duì)逆境脅迫密切相關(guān)[28]。其中GST是在動(dòng)植物中編碼的多功能蛋白的一個(gè)基因家族。同時(shí)GST也在參與干旱、鹽、重金屬等氧化應(yīng)激的反應(yīng)中發(fā)揮重要功能。過(guò)量的ROS誘導(dǎo)GST水平升高, 使得GST對(duì)外源物質(zhì)的解毒和應(yīng)激代謝中發(fā)揮關(guān)鍵作用[29]。Xu等[30]獲得了高表達(dá)番茄谷胱甘肽S-轉(zhuǎn)移酶的轉(zhuǎn)基因擬南芥, 命名為, 該轉(zhuǎn)基因擬南芥具有較強(qiáng)的耐鹽性和抗旱性, 同時(shí)也發(fā)現(xiàn)轉(zhuǎn)基因植株中脯氨酸、丙二醛含量以及抗氧化酶活性的變化與植物抗逆性相關(guān)。與野生型相比, 干旱處理過(guò)的擬南芥,啟動(dòng)子區(qū)重復(fù)序列的DNA甲基化水平顯著降低, 同時(shí)基因表達(dá)水平顯著升高[31]。本研究也得到類(lèi)似結(jié)果, 對(duì)抗旱性不同的馬鈴薯品種共同響應(yīng)干旱和DNA去甲基化處理下的差異表達(dá)基因進(jìn)行功能富集, 結(jié)果顯示, 有8個(gè)DEG顯著富集到植物谷胱甘肽代謝通路, 其中有6個(gè)屬于谷胱甘肽S轉(zhuǎn)移酶基因。同時(shí)差異基因啟動(dòng)子區(qū)分析表明, 在谷胱甘肽S轉(zhuǎn)移酶基因()的CpG島區(qū)域有與植物響應(yīng)干旱脅迫相關(guān)的ABRE和CAAT-box作用元件。推測(cè)馬鈴薯在響應(yīng)干旱脅迫時(shí), 可能通過(guò)DNA去甲基化來(lái)降低該基因啟動(dòng)子區(qū)的ABRE和CAAT-box作用元件的甲基化水平, 進(jìn)而激活該基因的表達(dá)以應(yīng)對(duì)干旱脅迫(圖7)。

以上結(jié)果表明, 利用比較轉(zhuǎn)錄學(xué)分析干旱和DNA去甲基化處理下差異表達(dá)基因, 可挖掘到通過(guò)DNA甲基化參與調(diào)控馬鈴薯響應(yīng)干旱脅迫的相關(guān)基因, 為研究馬鈴薯干旱脅迫響應(yīng)的表觀遺傳學(xué)機(jī)理提供新的研究思路。同時(shí), 本研究中有關(guān)DNA甲基化具體調(diào)控馬鈴薯響應(yīng)干旱脅迫的機(jī)制尚不清楚,后續(xù)研究可通過(guò)檢測(cè)差異表達(dá)基因啟動(dòng)子區(qū)DNA甲基化水平來(lái)進(jìn)一步解析其蘊(yùn)含的生物學(xué)意義。

4 結(jié)論

本研究對(duì)甘露醇和5-azadC處理下抗旱性不同的2個(gè)馬鈴薯品種進(jìn)行比較轉(zhuǎn)錄組分析發(fā)現(xiàn), 1345個(gè)DEG既響應(yīng)干旱又響應(yīng)DNA去甲基化。這些基因功能富集結(jié)果顯示, 與植物抗旱相關(guān)的通路有植物MAPK信號(hào)通路、谷胱甘肽代謝、糖酵解/糖異生、磷酸肌醇代謝和植物激素信號(hào)轉(zhuǎn)導(dǎo)通路。差異基因啟動(dòng)子分析表明, 干旱脅迫下GST基因通過(guò)DNA去甲基化來(lái)降低啟動(dòng)子區(qū)ABRE和CAAT-box作用元件的甲基化水平, 進(jìn)而激活該基因的表達(dá)以應(yīng)對(duì)干旱脅迫。本研究結(jié)果將為進(jìn)一步解析馬鈴薯響應(yīng)干旱脅迫的表觀遺傳學(xué)調(diào)控機(jī)理提供理論依據(jù)。

[1] 余林輝, 蔡曉騰, 徐萍, 向成斌. 植物抗旱節(jié)水: 從實(shí)驗(yàn)室到田間. 中國(guó)科學(xué): 生命科學(xué), 2017, 47: 145–154. Yu L H, Cai X T, Xu P, Xiang C B. Drought resistant and water-saving plants: from laboratory to field.(), 2017, 47: 145–154 (in Chinese with English abstract).

[2] Fang Y, Xiong L. General mechanisms of drought response and their application in drought resistance improvement in plants., 2015, 72: 673–689.

[3] 劉玉冰, 李新榮, 李蒙蒙, 劉丹, 張?chǎng)├? 中國(guó)干旱半干旱區(qū)荒漠植物葉片(或同化枝)表皮微形態(tài)特征. 植物生態(tài)學(xué)報(bào), 2016, 40: 1189–1207. Liu Y B, Li X R, Li M M, Liu D, Zhang W L. Leaf (or assimilation branch) epidermal micromorphology of desert plant in arid and semi-arid areas of China., 2016, 40: 1189–1207 (in Chinese with English abstract).

[4] 朱健康, 倪建平. 植物非生物脅迫信號(hào)轉(zhuǎn)導(dǎo)及應(yīng)答. 中國(guó)稻米, 2016, 22(6): 52–60. Zhu J K, Ni J P. Abiotic stress signaling and responses in plants., 2016, 22(6): 52–60 (in Chinese with English abstract).

[5] 余斌, 楊宏羽, 王麗, 劉玉匯, 白江平, 王蒂, 張俊蓮. 引進(jìn)馬鈴薯種質(zhì)資源在干旱半干旱區(qū)的表型性狀遺傳多樣性分析及綜合評(píng)價(jià). 作物學(xué)報(bào), 2018, 44: 63–74. Yu B, Yang H Y, Wang L, Liu Y H, Bai J P, Wang D, Zhang J L. Genetic diversity analysis and comprehensive assessment of phenotypic traits in introduced potato germplasm resources in arid and semi-arid area., 2018, 44: 63–74 (in Chinese with English abstract).

[6] Banerjee A, Roychoudhury A. Epigenetic regulation during salinity and drought stress in plants: histone modifications and DNA methylation., 2017, 11: 199–204.

[7] Pikaard C S, Mittelsten S O. Epigenetic regulation in plants., 2014, 6: a019315.

[8] Vanyushin B F, Ashapkin V V. DNA methylation in higher plants: past, present and future., 2011, 1809: 360–368.

[9] Gallusci P, Hodgman C, Teyssier E, Seymour G B. DNA methylation and chromatin regulation during fleshy fruit development and ripening., 2016, 7: 807.

[10] 趙云雷, 葉武威, 王俊娟, 樊保香, 宋麗艷. DNA甲基化與植物抗逆性研究進(jìn)展. 西北植物學(xué)報(bào), 2009, 29: 1479–1489. Zhao Y L, Ye W W, Wang J J, Fan B X, Song L Y. Review of DNA methylation and plant stress-tolerance., 2009, 29: 1479–1489 (in Chinese with English abstract).

[11] Fei Y, Xue Y, Du P, Yang S, Deng X. Expression analysis and promoter methylation under osmotic and salinity stress ofin wheat (L)., 2017, 254: 987–996.

[12] McLoughlin F, Arisz S A, Dekker H L, Kramer G, de Koster C G, Haring M A, Munnik T, Testerink C. Identification of novel candidate phosphatidic acid-binding proteins involved in the salt-stress response ofroots., 2013, 450: 573–581.

[13] Liang D, Zhang Z, Wu H, Huang C, Shuai P, Ye C Y, Tang S, Wang Y, Yang L, Wang J. Single-base-resolution methylomes of Populus trichocarpa reveal the association between DNA methylation and drought stress., 2014, 15: S9.

[14] Abid G, Mingeot D, Muhovski Y, Mergeai G, Aouida M, Abdelkarim S, Aroua I, El Ayed M, M’hamdi M, Sassi K. Analysis of DNA methylation patterns associated with drought stress response in faba bean (L.) using methylation-sensitive amplification polymorphism (MSAP)., 2017, 142: 34–44.

[15] 李鵬程, 畢真真, 梁文君, 孫超, 張俊蓮, 白江平. DNA甲基化參與調(diào)控馬鈴薯干旱脅迫響應(yīng). 作物學(xué)報(bào), 2019, 45: 1595–1603 Li P C, Bi Z Z, Liang W J, Sun C, Zhang J L, Bai J P. DNA methylation involved in regulating drought stress response of potato., 2019, 45: 1595–1603 (in Chinese with English abstract).

[16] Garg R, Shankar R, Thakkar B, Kudapa H, Krishnamurthy L, Mantri N, Varshney R K, Bhatia S, Jain M. Transcriptome analyses reveal genotype- and developmental stage-specific molecular responses to drought and salinity stresses in chickpea., 2016, 6: 19228.

[17] Zhou Y, Yang P, Cui F, Zhang F, Luo X, Xie J. Transcriptome analysis of salt stress responsiveness in the seedlings of Dongxiang wild rice (Griff.)., 2016, 11: e0146242.

[18] Raudvere U, Kolberg L, Kuzmin I, Arak T, Adler P, Peterson H, Vilo J. g:Profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update)., 2019, 47: W191–W198.

[19] Li L C, Dahiya R. MethPrimer: designing primers for methylation PCRs., 2002, 18: 1427–1431.

[20] Lescot M, Déhais P, Thijs G, Marchal K, Moreau Y, Van de Peer Y, Rouzé P, Rombauts S. PlantCARE, a database of plant cis-acting regulatory elements and a portal to tools for in silico analysis of promoter sequences., 2002, 30: 325–327.

[21] 李艷, 錢(qián)偉強(qiáng). 植物中DNA甲基化及去甲基化研究進(jìn)展. 生命科學(xué), 2017, 29: 302–309. Li Y, Qian W Q. Mechanisms of DNA methylation and demethylation in plants., 2017, 29: 302–309 (in Chinese with English abstract).

[22] Golldack D, Li C, Mohan H, Probst N. Tolerance to drought and salt stress in plants: unraveling the signaling networks., 2014, 5: 151.

[23] 王凱悅, 陳芳泉, 黃五星. 植物干旱脅迫響應(yīng)機(jī)制研究進(jìn)展. 中國(guó)農(nóng)業(yè)科技導(dǎo)報(bào), 2019, 21(2): 19–25. Wang K Y, Chen F Q, Huang W X. Research advance on drought stress response mechanism in plants., 2019, 21(2): 19–25 (in Chinese with English abstract).

[24] Molina C, Rotter B, Horres R, Udupa S M, Besser B, Bellarmino L, Baum M, Matsumura H, Terauchi R, Kahl G, Winter P. SuperSAGE: the drought stress-responsive transcriptome of chickpea roots., 2008, 9: 553.

[25] Ha C V, Leyva-Gonzalez M A, Osakabe Y, Tran U T, Nishiyama R, Watanabe Y, Tanaka M, Seki M, Yamaguchi S, Dong N V, Yamaguchi-Shinozaki K, Shinozaki K, Herrera-Estrella L, Tran L S. Positive regulatory role of strigolactone in plant responses to drought and salt stress., 2014, 111: 851–856.

[26] Verma V, Ravindran P, Kumar P P. Plant hormone-mediated regulation of stress responses., 2016, 16: 86.

[27] Jain M, Khurana J P. Transcript profiling reveals diverse roles of auxin-responsive genes during reproductive development and abiotic stress in rice., 2009, 276: 3148–3162.

[28] 陳坤明, 宮海軍, 王鎖民. 植物谷胱甘肽代謝與環(huán)境脅迫. 西北植物學(xué)報(bào), 2004, 24: 1119–1130. Chen K M, Gong H J, Wang S M. Glutathione metabolism and environmental stresses in plants., 2004, 24: 1119–1130 (in Chinese with English abstract).

[29] Li Z, Yu J, Peng Y, Huang B. Metabolic pathways regulated by abscisic acid, salicylic acid and gamma-aminobutyric acid in association with improved drought tolerance in creeping bentgrass ()., 2017, 159: 42–58.

[30] Xu J, Xing X J, Tian Y S, Peng R H, Xue Y, Zhao W, Yao Q H. Transgenicplants expressing tomato glutathione S-transferase showed enhanced resistance to salt and drought stress., 2015, 10: e0136960.

[31] 孟大偉, 王悅, 李沛璇, 趙宇威, 周瑤, 韓玉, 郎晨婧, 金太成, 楊麗萍. 干旱誘導(dǎo)基因DNA去甲基化. 分子植物育種, 2020, 18: 6108–6113.Meng D W, Wang Y, Li P X, Zhao Y W, Zhou Y, Han Y, Lang C J, Jin T C, Yang L P. Drought-introduced DNA demethylation ofgene., 2020, 18: 6108–6113 (in Chinese with English abstract).

Key genes mining of DNA methylation involved in regulating drought stress response in potato

LI Peng-Cheng1,2,**, BI Zhen-Zhen1,2,**, SUN Chao1,2, QIN Tian-Yuan1,2, LIANG Wen-Jun1,2, WANG Yi-Hao1,2, XU De-Rong1,2, LIU Yu-Hui1,2, ZHANG Jun-Lian1,2, and BAI Jiang-Ping1,2,*

1Gansu Provincial Key Laboratory of Aridland Crop Science / Gansu Key Laboratory of Crop Improvement & Germplasm Enhancement, Lanzhou 730070, Gansu, China;2College of Agronomy, Gansu Agricultural University, Lanzhou 730070, Gansu, China

When plants are subjected to water stress, they will make a rapid response to drought stress through DNA methylation. In order to study how DNA methylation affects the transcriptional expression of genes under drought stress in potato, comparative transcriptomic analysis was carried out on two potato varieties with different drought resistances, which were planted under mannitol simulated drought and 5-azadC treatments. The differentially expressed genes (DEGs) were identified by Fold-change > 2 and corrected< 0.01. Then DEGs were subjected to GO enrichment analysis. The results showed that these DEGs were enriched in the oxidative stress and carbohydrate metabolism, suggesting these GO term-related genes were also regulated by DNA demethylation, responded to drought stress in different drought-tolerant potatoes. The common 1345 DEGs both responding to drought stress and DNA demethylation were functionally enriched by KEGG pathway. The results showed that plant MAPK signal pathway, plant hormone signal transduction pathway, plant glutathione metabolism pathway, glycolysis and glutathione metabolism pathway and inositol phosphate metabolism pathway were related to plant drought resistance. It was suggested that the sensitivity of these pathway-related genes responding to drought stress were regulated by DNA methylation in Atlantic and Qingshu 9. The cis-acting elements and methylated CpG islands were analyzed in the 1500 bp promoter region of DEGs. The results showed that the methylation level of ABRE and CAAT-box acting elements in the promoter region of GST gene involved in plant glutathione metabolism were reduced through DNA demethylation under drought stress. Then the expression was activated in response to drought stress. Therefore, DEGs under drought and DNA demethylation treatments could be analyzed using comparative transcriptomic, and then the genes related to DNA methylation involved in regulating drought stress response in potato could be found. These results provide a new idea for further studying the epigenetic mechanism of drought stress response in potato.

potato; drought stress; DNA methylation; transcriptome; glutathione S transferase

10.3724/SP.J.1006.2021.04152

本研究由國(guó)家自然科學(xué)基金項(xiàng)目(31660432, 31960442), 國(guó)家現(xiàn)代農(nóng)業(yè)產(chǎn)業(yè)技術(shù)體系建設(shè)專(zhuān)項(xiàng)(CARS-09-P14), 甘肅省馬鈴薯產(chǎn)業(yè)體系(GARS-03-P1), 甘肅農(nóng)業(yè)大學(xué)引進(jìn)人才專(zhuān)項(xiàng)科研項(xiàng)目(2017RCZX-01)和甘肅省科技廳項(xiàng)目(18JR3RA170)資助。

This study was supported by the National Natural Science Foundation of China (31660432, 31960442), the China Agriculture Research System (CARS-09-P14), the Gansu Potato Industry System (GARS-03-P1), the Gansu Agricultural University Special Talent Research Project (2017RCZX-01), and the Gansu Science and Technology Fund (18JR3RA170).

白江平, E-mail:baijp@gsau.edu.cn

**同等貢獻(xiàn)(Contributed equally to this work)

李鵬程, E-mail: 526040572@qq.com; 畢真真, E-mail: bizhen925@sina.com

2020-07-11;

2020-10-14;

2020-11-06.

URL: https://kns.cnki.net/kcms/detail/11.1809.S.20201105.1124.004.html