李慧，張雨峰,2，李曉剛，王中華，藺經(jīng)，常有宏

全基因組DNA甲基化和轉(zhuǎn)錄組聯(lián)合分析鑒定杜梨耐鹽相關(guān)轉(zhuǎn)錄因子

李慧1，張雨峰1,2，李曉剛1，王中華1，藺經(jīng)1，常有宏1

1江蘇省農(nóng)業(yè)科學(xué)院果樹研究所/江蘇省高效園藝作物遺傳改良重點(diǎn)實(shí)驗(yàn)室，南京 210014；2南京林業(yè)大學(xué)生物與環(huán)境學(xué)院，南京 210037

【】鑒定不同杜梨株系根中響應(yīng)鹽脅迫信號的相關(guān)轉(zhuǎn)錄因子，分析鹽脅迫下基因序列DNA甲基化變化與基因表達(dá)改變之間的關(guān)系，探討參與調(diào)控不同杜梨株系耐鹽能力的轉(zhuǎn)錄因子成員?！尽恳远爬婺望}株系和普通株系為試材，在苗期使用200 mmol?L-1NaCl對90日齡組培生根苗進(jìn)行水培處理，以Hoagland營養(yǎng)液為對照。利用火焰石墨爐原子吸收光譜儀測定鈉離子含量；利用全基因組DNA甲基化和轉(zhuǎn)錄組測序技術(shù)從表觀遺傳修飾和轉(zhuǎn)錄調(diào)控水平對鹽脅迫下轉(zhuǎn)錄因子進(jìn)行生物信息學(xué)分析；最后用McrBC-PCR和qPCR對差異轉(zhuǎn)錄因子進(jìn)行驗(yàn)證。【】外源NaCl處理24 h后，杜梨植株中鈉離子含量顯著增加，其中耐鹽株系的增加幅度比普通株系小，為普通株系鈉含量的73.1%，但根中積累的鈉是普通株系含量的1.1倍；杜梨根中檢測到69類共2 682個(gè)轉(zhuǎn)錄因子的表達(dá)，鹽脅迫后243個(gè)轉(zhuǎn)錄因子在兩個(gè)株系中都發(fā)生了差異表達(dá)，包括AP2/ERF（37個(gè)）、bHLH（19個(gè)）、bZIP（7個(gè)）、HD-Zip（10個(gè)）、MYB（30個(gè)）、NAC（18個(gè)）、WRKY（8個(gè)）和ZFP（23個(gè)）等家族成員；鹽脅迫后，耐鹽株系基因組中轉(zhuǎn)錄因子甲基化水平下降，而普通株系轉(zhuǎn)錄因子甲基化水平上升，其發(fā)生DNA差異甲基化區(qū)域主要在基因啟動子位置，差異甲基化類型主要為mCHH，占mCG、mCHG、mCHH三種類型總和的93%以上。AP2/ERF、bHLH、DREB、GRAS、GT因子、HB Zip、MYB、NAC、Trihelix和Zinc finger ZFP家族的23個(gè)轉(zhuǎn)錄因子響應(yīng)鹽脅迫表達(dá)量上調(diào)而甲基化水平降低，可能參與調(diào)節(jié)鈉在根中的吸收和積累。對部分候選基因的表達(dá)模式和啟動子區(qū)域分別進(jìn)行實(shí)時(shí)熒光定量（qPCR）和甲基化依賴型限制性內(nèi)切酶PCR（McrBC-PCR），驗(yàn)證了生物信息學(xué)分析結(jié)果。【】鹽脅迫后在兩個(gè)杜梨株系根中均差異表達(dá)的轉(zhuǎn)錄因子數(shù)目為243個(gè)，其中8個(gè)轉(zhuǎn)錄因子（、、、、、、、）DNA序列的甲基化改變與基因轉(zhuǎn)錄水平變化呈負(fù)相關(guān)，研究結(jié)果為揭示轉(zhuǎn)錄因子參與不同杜梨株系耐鹽能力調(diào)控的分子機(jī)制提供了依據(jù)。

杜梨；鹽脅迫；轉(zhuǎn)錄因子；DNA甲基化；轉(zhuǎn)錄組

0 引言

【研究意義】由于極端氣候變化和不合理的農(nóng)業(yè)活動，全球約7%面積的土地（超過9億hm2）受到不同程度的鹽漬化危害[1]。耐鹽新品種的選育是減少土壤鹽漬化對農(nóng)業(yè)生產(chǎn)不良影響的有效途徑[2]，研究不同物種應(yīng)對鹽脅迫的具體機(jī)制是促進(jìn)耐鹽作物育種、合理利用鹽漬化土地的首要前提。杜梨是我國第三大果樹梨的常用抗逆砧木，用其作砧木嫁接，可以提高梨品種的耐鹽能力。發(fā)掘杜梨耐鹽相關(guān)基因，研究該物種的耐鹽機(jī)制對促進(jìn)梨生產(chǎn)及鹽漬土地的開發(fā)利用具有重要意義。【前人研究進(jìn)展】鹽脅迫主要產(chǎn)生滲透脅迫、離子脅迫和次生脅迫，植物在基因水平上進(jìn)行轉(zhuǎn)錄調(diào)控是應(yīng)對上述脅迫的關(guān)鍵步驟，轉(zhuǎn)錄因子在應(yīng)激反應(yīng)基因的表達(dá)中起著關(guān)鍵的調(diào)控作用[3-4]。利用微陣列技術(shù)分析鹽脅迫條件下擬南芥轉(zhuǎn)錄組的變化，發(fā)現(xiàn)有289個(gè)轉(zhuǎn)錄因子受鹽脅迫誘導(dǎo)上調(diào)表達(dá)，139個(gè)轉(zhuǎn)錄因子的表達(dá)被抑制[5]。陸地棉根中有123個(gè)轉(zhuǎn)錄因子基因受鹽脅迫誘導(dǎo)表達(dá)[6]，苜蓿中受鹽脅迫誘導(dǎo)而上調(diào)表達(dá)的基因中20%為轉(zhuǎn)錄因子[7]。同一物種中不同耐鹽能力的品種之間轉(zhuǎn)錄因子對鹽脅迫的轉(zhuǎn)錄響應(yīng)也存在差異，如鹽脅迫下兩個(gè)甜瓜品種差異表達(dá)的轉(zhuǎn)錄因子既表現(xiàn)特異性，也存在部分重疊[8]。除了基因轉(zhuǎn)錄調(diào)控外，植物還改變?nèi)蚪MDNA甲基化狀態(tài)以應(yīng)對鹽脅迫[2]，植株基因型是影響鹽脅迫下基因組DNA甲基化多態(tài)性改變的主要因素[9]。例如：水稻鹽耐受品種通過去甲基化來響應(yīng)鹽脅迫，其甲基化水平比鹽敏感品種更低[10]。而棉花的鹽耐受品種比敏感品種的DNA甲基化水平更高[11]。植物在高鹽逆境下DNA甲基化狀況的自我調(diào)節(jié)，是應(yīng)對鹽脅迫、調(diào)控基因時(shí)空表達(dá)所必需的表觀遺傳適應(yīng)機(jī)制。已有研究發(fā)現(xiàn)轉(zhuǎn)錄因子通過改變自身基因序列的DNA甲基化狀況來響應(yīng)鹽脅迫，例如水稻14個(gè)鋅指蛋白基因的甲基化狀態(tài)受鹽脅迫誘導(dǎo)改變，從而激活或抑制該類基因的轉(zhuǎn)錄，參與調(diào)控植株應(yīng)對高鹽逆境[12]。水稻基因組中的基因型特異性表觀遺傳變化可能是通過改變鹽脅迫響應(yīng)基因的表達(dá)網(wǎng)絡(luò)來感知和響應(yīng)鹽脅迫的重要替代調(diào)控機(jī)制[13]。而某些轉(zhuǎn)錄因子的DNA甲基化對基因表達(dá)的激活或抑制在大豆鹽脅迫耐受過程起關(guān)鍵作用[14]。【本研究切入點(diǎn)】目前杜梨耐鹽研究集中于生理評價(jià)和基因表達(dá)分析[15-18]，該物種中不同轉(zhuǎn)錄因子對鹽脅迫存在表達(dá)響應(yīng)[19-20]，但缺乏鹽脅迫下所有轉(zhuǎn)錄因子表達(dá)情況的全面分析，其基因DNA序列甲基化狀況對轉(zhuǎn)錄水平的影響也不清楚?！緮M解決的關(guān)鍵問題】以耐鹽能力不同的2個(gè)典型杜梨株系為材料，擬從轉(zhuǎn)錄組和全基因組甲基化聯(lián)合分析入手，系統(tǒng)研究杜梨鹽脅迫下轉(zhuǎn)錄因子表達(dá)調(diào)控情況，結(jié)合DNA甲基化狀態(tài)變化情況，鑒定杜梨耐鹽相關(guān)轉(zhuǎn)錄因子，為耐鹽品種培育和鹽堿地高效利用提供技術(shù)參考。

1 材料與方法

試驗(yàn)于2019—2022年在江蘇省農(nóng)業(yè)科學(xué)院果樹研究所進(jìn)行。

1.1 試驗(yàn)材料處理與鈉離子測定

以杜梨耐鹽株系（收集于連云港花果山，可耐6‰ NaCl）和普通株系（收集于南京紫金山，僅耐3‰ NaCl）組培苗為試材，經(jīng)生根培養(yǎng)基誘導(dǎo)生根30 d后，轉(zhuǎn)入溫室土培（營養(yǎng)土﹕蛭石3﹕1）90 d，選擇生長一致的植株，輕輕洗凈根部基質(zhì)，兩種株系各60 株分為兩組，一組加入額外添加200 mmol?L-1NaCl的Hoagland營養(yǎng)液進(jìn)行脅迫處理，一組加入Hoagland營養(yǎng)液作為對照[21]。處理24 h后分別收集植株全株或單獨(dú)的根系，去離子水充分沖洗后，全株分為地上部和根系兩部分烘干，用于鈉離子測定；而新鮮根系經(jīng)液氮急凍后保存于-80℃超低溫冰箱待用。

杜梨地上部和根系經(jīng)80℃烘干6 h，每0.5 g樣品加入25 mL醋酸（1 mol·L-1），90℃振蕩消化2 h，過濾定容后于火焰石墨爐原子吸收光譜儀（ZEEnit?700P，德國耶拿）測定鈉離子含量，全株的鈉離子含量由地上部和根系兩部分含量相加。鈉標(biāo)準(zhǔn)溶液購自美國默克公司。

1.2 轉(zhuǎn)錄組測序與數(shù)據(jù)分析

使用TRIzol 試劑（Invitrogen）分離和純化杜梨根總RNA。建庫RNA滿足濃度＞100 ng·μL-1、RIN number＞7.0、OD260/280＞1.8、總含量＞20 μg。NEBNext?Ultra? RNA Library Prep Kit for Illumina?（NEB）生成測序文庫及Illumina HiSeqTM2500雙端測序均在南京集思慧遠(yuǎn)生物技術(shù)有限公司進(jìn)行，每個(gè)處理的樣品均進(jìn)行3次生物學(xué)重復(fù)。

使用FASTQC（http://www.bioinformatics.babraham. ac.uk/projects/fastqc/）讀取質(zhì)量控制指標(biāo)，NGS QC工具包用于去除接頭序列和低質(zhì)量序列[22]。使用HISAT（v2.1.0）[23]將過濾后的干凈序列比對至杜梨參考基因組[24]，參數(shù)設(shè)置為默認(rèn)標(biāo)準(zhǔn)。獲得的比對序列組成bam文件計(jì)算FPKM（fragments per kilobases per million reads）值來代表基因表達(dá)水平。R包DESeq2[25]用于鑒定差異表達(dá)基因，標(biāo)準(zhǔn)為|log2Fold change|＞1且FDR＜0.05。使用iTAK（1.7a）軟件從杜梨基因組的所有基因中鑒定轉(zhuǎn)錄因子[26]。

1.3 DNA甲基化測序與數(shù)據(jù)分析

使用 DNeasy Plant Mini Kit（Qiagen）提取杜梨根DNA，MinElute PCR Cleanup（Qiagen）對DNA樣品進(jìn)行純化并測定濃度。1 μg基因組DNA經(jīng)超聲波處理機(jī)械打斷、末端修復(fù)、3'端加A堿基、連接甲基化接頭后，使用EZ DNA Methylation-Gold?試劑盒（Zymo Research）將DNA轉(zhuǎn)化為亞硫酸氫鹽，全基因組重亞硫酸鹽測序（whole-genome bisulfite sequencing，WGBS）文庫的構(gòu)建以及Illumina HiSeqTM2500雙端測序均在南京集思慧遠(yuǎn)生物技術(shù)有限公司進(jìn)行，每組樣品均進(jìn)行3次生物學(xué)重復(fù)。

使用Trim fastp（v0.20.0）[27]對原始測序讀數(shù)進(jìn)行質(zhì)量過濾和接頭去除。隨后使用Bismark（v0.22.3）- bowtie2（2.3.5.1）[28]將高質(zhì)量修剪后的讀段與杜梨參考基因組[24]進(jìn)行比對。使用Bismark（命令為—no_ overlap--ignore 5--ignore_r2 5）識別比對讀數(shù)中的甲基化胞嘧啶，5mC檢測的篩選條件為：覆蓋度≥4X、FDR＜0.05。使用MethylKit（v1.12.0）[29]在R包（v3.6.0）中識別差異甲基化區(qū)域（differentially methylated region，DMR），參數(shù)為span=1 000、FDR＜0.05、C＞0.2、CG＞0.25、CHG＞0.20、CHH＞0.1。使用Bedtools（v2.21.0）[30]對DMR進(jìn)行基因特征注釋。

1.4 轉(zhuǎn)錄組和甲基化聯(lián)合分析

使用Bedtools（v2.21.0）在R包（v3.6.0）中注釋來自WGBS分析的DMR區(qū)域[30]。使用Perl編程語言獲得目標(biāo)基因進(jìn)行分析，它們是在指定位置與差異甲基化區(qū)域（DMR）相交的轉(zhuǎn)錄組結(jié)果的差異甲基化基因（differentially methylated gene，DMG）。最后，結(jié)合iTAK(1.7a)軟件鑒定出的轉(zhuǎn)錄因子ID號信息確定DMG中包含的轉(zhuǎn)錄因子成員[26]。

1.5 McrBC-PCR和qPCR

qPCR 基因特異性引物使用Primer Premier 5.0設(shè)計(jì)，上海生工合成，引物序列見表1。所有引物均用PCR擴(kuò)增、電泳和溶解曲線進(jìn)行測試以保證引物特異性。使用LightCycler?480 II (Roche)進(jìn)行qPCR，反應(yīng)體系按Genious 2×SYBR Green Fast qPCR Mix (AB clonal)說明書進(jìn)行。以為內(nèi)參，用2-ΔΔCT公式計(jì)算相對表達(dá)量[16,31]。

表1 實(shí)時(shí)熒光定量PCR和McrBC-qPCR引物

用核酸內(nèi)切酶McrBC（M0272S，New England Biolabs）按照說明書的方法對1 μg的杜梨根基因組DNA進(jìn)行酶切消化。未加GTP的酶切消化反應(yīng)作為陰性對照用于標(biāo)準(zhǔn)化分析。甲基化的DNA可以被McrBC消化，PCR結(jié)果表現(xiàn)為qPCR信號水平與甲基化水平呈負(fù)相關(guān)[32]。McrBC qPCR引物設(shè)計(jì)和合成同qPCR分析，引物序列見表1。

1.6 數(shù)據(jù)分析

試驗(yàn)數(shù)據(jù)采用SPSS20進(jìn)行ANOVA分析，不同處理之間差異采用Duncan檢測，數(shù)值表示為3個(gè)生物學(xué)重復(fù)的平均值±標(biāo)準(zhǔn)誤，不同小寫字母表示顯著性差異（＜0.05）。

2 結(jié)果

2.1 不同杜梨株系鈉離子含量

在正常生長條件下，耐鹽株系和普通株系全株鈉離子含量無顯著差異，它們根系中鈉離子含量占全株的比例相差不大；200 mmol·L-1NaCl處理24 h后，耐鹽株系全株鈉離子含量是未處理前的3.41倍，根系中鈉離子含量占全株的47.09%；普通株系全株鈉離子含量是未處理前的4.52倍，根系中鈉離子含量僅占全株的35.79%（表2），表明鹽脅迫條件下耐鹽株系能夠?qū)⑤^多的鈉離子儲存在根中，限制其向上運(yùn)輸，從而達(dá)到鹽耐受的目的。

表2 不同杜梨株系鹽脅迫后鈉離子含量

不同小寫字母表示不同處理間差異顯著（＜0.05）

Different lowercase letters indicate significant difference (＜0.05)

2.2 不同杜梨株系表達(dá)差異轉(zhuǎn)錄因子

從12個(gè)杜梨根系樣品中分別獲得2 069.3—3 025.4萬條不同數(shù)目的轉(zhuǎn)錄組有效序列，Q30≥90.8%，各樣品數(shù)據(jù)量≥6.2 G，共獲得96.4 Gb數(shù)據(jù)供后續(xù)分析（NCBI number: PRJNA812627），每個(gè)處理的3次生物學(xué)重復(fù)樣品間的一致性為91.1%—98.9%。轉(zhuǎn)錄組數(shù)據(jù)與杜梨參考基因組(http://bigd.big.ac.cn/gwh/ Assembly/505/show)的比對率為84.8%—86.9%。借助iTAK（1.7a）軟件從杜梨基因組中鑒定出69類共2 682個(gè)轉(zhuǎn)錄因子基因（圖1）。將200 mmol·L-1NaCl處理24 h后的杜梨耐鹽株系或普通株系根轉(zhuǎn)錄組數(shù)據(jù)分別和其對照根轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行比對分析，分別獲得702個(gè)和435個(gè)差異表達(dá)轉(zhuǎn)錄因子基因（圖2）；其中耐鹽株系處理24 h后，共計(jì)342個(gè)轉(zhuǎn)錄因子基因表達(dá)上調(diào)，360個(gè)轉(zhuǎn)錄因子基因表達(dá)下調(diào)；普通株系處理24 h后，共計(jì)192個(gè)轉(zhuǎn)錄因子基因表達(dá)上調(diào)，243個(gè)轉(zhuǎn)錄因子基因表達(dá)下調(diào)；耐鹽株系和普通株系鹽脅迫前后共同的表達(dá)差異轉(zhuǎn)錄因子基因有243個(gè)，包括AP2/ERF（37個(gè)）、bHLH（19個(gè)）、bZIP（7個(gè)）、HD-Zip（10個(gè)）、MYB（30個(gè)）、NAC（18個(gè)）、WRKY（8個(gè)）和ZFP（23個(gè)）等基因家族成員；其中表達(dá)趨勢相同的轉(zhuǎn)錄因子基因?yàn)?31個(gè)（115個(gè)上調(diào)，116個(gè)下調(diào)）；表達(dá)趨勢不同的差異轉(zhuǎn)錄因子基因有12個(gè)。

2.3 不同杜梨株系轉(zhuǎn)錄因子差異甲基化區(qū)域

12個(gè)杜梨根系樣品的WGBS測序原始數(shù)據(jù)過濾后得到6 411.3—8 305.5萬條有效序列，Q30≥91.1%，分別獲得19.2—24.9 Gb數(shù)據(jù)（NCBI number: PRJNA812739），是杜梨參考基因組[24]大小（532.7 Mb，http://bigd.big.ac.cn/gwh/Assembly/505/show）的36.1— 46.8倍，WGBS數(shù)據(jù)與杜梨參考基因組的比對率為64.4%—72.8%。不同樣品之間轉(zhuǎn)錄因子相同基因區(qū)域的甲基化水平類似，以轉(zhuǎn)錄起始位點(diǎn)為界限，基因編碼區(qū)的甲基化程度最低，基因上游（2 kb）及下游（2 kb）區(qū)域甲基化程度較高；鹽脅迫后，耐鹽株系基因組中轉(zhuǎn)錄因子甲基化水平下降，而普通株系轉(zhuǎn)錄因子甲基化水平上升；兩個(gè)株系發(fā)生DNA差異甲基化的區(qū)域均主要為基因啟動子序列（圖3）。比較杜梨耐鹽株系和普通株系鹽脅迫前后WGBS數(shù)據(jù)，分別獲得轉(zhuǎn)錄因子基因序列中差異甲基化區(qū)域（包括mC、mCG、mCHG、mCHH 4種類型）889和857個(gè)，其中mCHH類型的差異甲基化區(qū)域占所有差異甲基化區(qū)域數(shù)目的87.7%和86.0%。耐鹽株系鹽脅迫處理24 h后，轉(zhuǎn)錄因子基因序列中386個(gè)差異甲基化區(qū)域甲基化水平升高，503個(gè)差異甲基化區(qū)域的甲基化水平降低；普通株系鹽脅迫處理24 h后，轉(zhuǎn)錄因子基因序列中569個(gè)差異甲基化區(qū)域的甲基化水平升高，288個(gè)差異甲基化區(qū)域的甲基化水平降低（圖4、圖5）。

圖1 杜梨根中檢測到的轉(zhuǎn)錄因子種類及比例

DR：耐鹽株系對照根；DNR：耐鹽株系200 mmol·L-1 NaCl處理24 h的根；UR：普通株系對照根；UNR：普通株系200 mmol·L-1 NaCl處理24 h的根。下同

圖3 杜梨根中轉(zhuǎn)錄因子不同基因區(qū)域的甲基化水平

2.4 聯(lián)合分析與驗(yàn)證

將WGBS和RNA-Seq數(shù)據(jù)進(jìn)行聯(lián)合分析，發(fā)現(xiàn)鹽脅迫后耐鹽株系中同時(shí)發(fā)生表達(dá)差異和甲基化差異的轉(zhuǎn)錄因子基因共203個(gè)（包括275個(gè)差異甲基化區(qū)域），其中負(fù)相關(guān)的轉(zhuǎn)錄因子基因?yàn)?8個(gè)，而普通株系中表達(dá)差異和甲基化差異相關(guān)的轉(zhuǎn)錄因

圖5 鹽脅迫后杜梨不同株系根中在轉(zhuǎn)錄因子上的差異甲基化區(qū)域熱圖

子基因共121個(gè)（包括178個(gè)差異甲基化區(qū)域），其中負(fù)相關(guān)的轉(zhuǎn)錄因子基因?yàn)?0個(gè)（表3）。聯(lián)合分析獲得的鹽脅迫相關(guān)差異轉(zhuǎn)錄因子基因中，耐鹽株系和普通株系共有的mC類甲基化差異轉(zhuǎn)錄因子基因只有1個(gè)（GWHGAAYT031490），無相同的mCG類或mCHG類甲基化差異轉(zhuǎn)錄因子基因，而它們共有的mCHH類甲基化差異轉(zhuǎn)錄因子基因?yàn)?3個(gè)。

表3 聯(lián)合分析杜梨不同株系鹽脅迫后差異轉(zhuǎn)錄因子數(shù)量

*：基因區(qū)域中差異甲基化區(qū)域和差異表達(dá)基因的交集數(shù)量，括號表示基因數(shù)量；**：啟動子區(qū)域中差異甲基化區(qū)域和差異表達(dá)基因的交集數(shù)量，括號表示基因數(shù)量

*: The number of intersections between differentially methylated regions and differentially expressed genes in the gene region, the number of genes was showed in the brackets; **: The number of intersections of differentially methylated regions and differentially expressed genes in the promoter region, the number of genes was showed in the brackets

選擇鹽脅迫后基因轉(zhuǎn)錄變化與甲基化狀態(tài)相關(guān)的12個(gè)轉(zhuǎn)錄因子（IGV軟件展示差異甲基化區(qū)域，附圖1），對它們在兩個(gè)株系中的表達(dá)水平和差異甲基化區(qū)域進(jìn)行qPCR和McrBC-qPCR驗(yàn)證。結(jié)果表明，鹽脅迫后(GWHGAAYT045012)、(GWHGAAYT011887)、(GWHGAAYT041697)、(GWHGAAYT007070)、(GWHGAAYT038287)在耐鹽單株和普通單株中表達(dá)水平上升，并且檢測區(qū)域的甲基化水平下降，其中耐鹽單株的變化幅度遠(yuǎn)大于普通單株；而（GWHGAAYT020116）、（GWHGAAYT056118）、（GWHGAAYT 009070）在耐鹽單株和普通單株中表達(dá)水平上升，并且檢測區(qū)域的甲基化水平下降，但是它們在普通單株中表達(dá)水平和甲基化水平變化不顯著；（GWHGAAYT023975）、（GWHGAAYT016435）在耐鹽單株和普通單株中表達(dá)水平上升，而檢測區(qū)域的甲基化水平在耐鹽單株中下降、普通單株中上升；（GWHGAAYT040185）、（GWHGAAYT 025897）在耐鹽單株和普通單株中表達(dá)水平上升，而檢測區(qū)域的甲基化水平在耐鹽單株中無顯著變化、在普通單株中下降（圖6、圖7）。上述結(jié)果說明鹽脅迫誘導(dǎo)的DNA甲基化變異在不同的耐鹽單株中均可調(diào)控部分相關(guān)轉(zhuǎn)錄因子基因（、、、、）的表達(dá)變化，而有些響應(yīng)鹽脅迫的差異表達(dá)基因（、、）可能并不直接受DNA甲基化所調(diào)控，表現(xiàn)為它們的轉(zhuǎn)錄水平上升與甲基化變化并不是負(fù)相關(guān)。在鹽脅迫條件下，DNA甲基化如何精確調(diào)控特定基因表達(dá)，影響杜梨不同株系耐鹽能力差異的具體調(diào)控機(jī)制仍需進(jìn)一步探究。

3 討 論

3.1 轉(zhuǎn)錄因子參與杜梨耐鹽調(diào)控過程

前人研究發(fā)現(xiàn)，全基因組水平上發(fā)生轉(zhuǎn)錄改變是植物受到鹽脅迫后的普遍反應(yīng)，而轉(zhuǎn)錄因子作為RNA轉(zhuǎn)錄網(wǎng)絡(luò)調(diào)控的樞紐，它們在高鹽逆境條件下表達(dá)水平的變化尤為常見[3-4]。例如擬南芥、陸地棉和苜蓿在鹽脅迫條件下，都檢測到大量轉(zhuǎn)錄因子的表達(dá)量提高或降低[5-7]。本研究結(jié)果表明，杜梨根系在鹽脅迫后16.2%—26.2%的轉(zhuǎn)錄因子基因受誘導(dǎo)上調(diào)或下調(diào)表達(dá)。除了基因轉(zhuǎn)錄調(diào)控之外，全基因組DNA甲基化狀態(tài)改變也是應(yīng)對鹽脅迫的重要手段，植物可通過DNA低甲基化或高甲基化來調(diào)節(jié)基因表達(dá)以適應(yīng)高鹽逆境[2,33-35]，NaCl處理后，杜梨兩個(gè)株系根的全基因組DNA甲基化水平均發(fā)生了變化[21]，表明該物種也存在應(yīng)對鹽脅迫逆境的表觀遺傳調(diào)控機(jī)制。已有研究發(fā)現(xiàn)轉(zhuǎn)錄因子通過改變自身基因序列的DNA甲基化狀況來響應(yīng)鹽脅迫，比如水稻、大豆和苜蓿等植物中AP2/ERF、NAC、ZFP轉(zhuǎn)錄因子家族成員可改變自身序列的DNA甲基化水平來激活或抑制基因表達(dá)，從而參與鹽脅迫應(yīng)答調(diào)控，在鹽耐受過程起關(guān)鍵作用[12-14,36]。本研究發(fā)現(xiàn)，杜梨根系中轉(zhuǎn)錄因子基因DNA序列中差異甲基化區(qū)域（包括mC、mCG、mCHG、mCHH四種類型）為857—889個(gè)，同時(shí)發(fā)生表達(dá)差異和甲基化差異的轉(zhuǎn)錄因子基因?yàn)?21—203個(gè)，包括AP2/ERF、bHLH、NAC、MYB、WRKY等家族成員，并且它們的轉(zhuǎn)錄水平受鹽脅迫影響，其中23個(gè)轉(zhuǎn)錄因子響應(yīng)鹽脅迫表達(dá)量上調(diào)而甲基化水平降低，可能參與調(diào)節(jié)鈉在根中的吸收和積累，表明杜梨與其他植物類似，其根系中的轉(zhuǎn)錄因子也可通過改變自身基因序列DNA甲基化水平來調(diào)節(jié)轉(zhuǎn)錄情況，從而參與植株的耐鹽調(diào)控過程。

圖7 鹽脅迫下杜梨不同株系根中轉(zhuǎn)錄因子DNA甲基化水平和差異甲基化區(qū)域驗(yàn)證

3.2 杜梨響應(yīng)鹽脅迫存在基因型特異性

植物應(yīng)對高鹽脅迫的復(fù)雜過程涉及脅迫相關(guān)基因的高效有序表達(dá)，離不開體內(nèi)各類轉(zhuǎn)錄因子的有效調(diào)控[5-8]，同一物種中不同耐鹽能力的品種之間對鹽害的耐受能力主要與基因?qū)γ{迫的響應(yīng)程度相關(guān)[37]，例如鹽脅迫下兩個(gè)甜瓜品種差異表達(dá)的轉(zhuǎn)錄因子既表現(xiàn)特異性，也存在部分重疊[8]。杜梨耐鹽株系和普通株系鹽脅迫前后差異表達(dá)的轉(zhuǎn)錄因子也存在相同特征，各有特異表達(dá)的成員，也有部分重疊的成員，某一具體轉(zhuǎn)錄因子基因在杜梨鹽脅迫過程所起的作用需要進(jìn)一步試驗(yàn)。此外，轉(zhuǎn)錄因子在高鹽逆境條件下表達(dá)水平的變化往往伴隨著基因DNA序列的甲基化模式改變[5-8]，水稻、大豆和苜蓿等植物中AP2/ERF、NAC、ZFP轉(zhuǎn)錄因子家族成員在鹽脅迫下DNA甲基化改變和轉(zhuǎn)錄水平變化與品種特性密切相關(guān)，同一物種不同基因型之間存在特異性的表觀遺傳變化，從而在不同程度上改變鹽脅迫響應(yīng)基因的表達(dá)網(wǎng)絡(luò)來感知和響應(yīng)鹽脅迫，最后表現(xiàn)為植株對鹽脅迫的不同耐受能力[12-14,36]。NaCl處理后，杜梨兩個(gè)株系根的全基因組DNA甲基化水平均發(fā)生了變化，普通單株整體甲基化水平上升，耐鹽株系通過去甲基化來響應(yīng)鹽脅迫，表現(xiàn)為全基因組甲基化水平下降[21]。本研究發(fā)現(xiàn)杜梨耐鹽單株根系比普通單株根系DNA甲基化程度低，可能與其能在根系中積累更多的鈉離子，從而更好地適應(yīng)高鹽逆境相關(guān)。此外，DNA甲基化對基因表達(dá)的影響因組織類型、甲基化序列背景、基因間區(qū)和基因體內(nèi)甲基化區(qū)域而異，啟動子區(qū)域的DNA甲基化通常會導(dǎo)致基因表達(dá)量下降[38]。本研究中杜梨耐鹽株系基因組中轉(zhuǎn)錄因子DNA甲基化水平下降，而普通株系轉(zhuǎn)錄因子DNA甲基化水平上升；兩個(gè)株系發(fā)生DNA差異甲基化的區(qū)域主要為轉(zhuǎn)錄因子基因的啟動子序列，而前者的轉(zhuǎn)錄因子表達(dá)水平普遍高于后者，與上述結(jié)論相符。杜梨耐鹽單株和普通單株中篩選獲得共有的鹽脅迫相關(guān)差異轉(zhuǎn)錄因子基因64個(gè)，包括AP2/ERF、NAC、ZFP等家族成員，這些基因的表達(dá)水平和啟動子序列DNA甲基化水平密切相關(guān)，但是它們在兩個(gè)株系中的表現(xiàn)有所差異。表明杜梨在應(yīng)對鹽脅迫時(shí)不同株系的根系中轉(zhuǎn)錄因子基因的表達(dá)情況和DNA甲基化模式不盡相同，該物種的逆境表觀遺傳調(diào)控機(jī)制存在基因型特征，而某一具體的轉(zhuǎn)錄因子成員如何通過DNA甲基化來調(diào)控基因轉(zhuǎn)錄的詳細(xì)過程，進(jìn)而參與植株適應(yīng)高鹽逆境的表觀遺傳機(jī)制仍需進(jìn)一步驗(yàn)證。

4 結(jié)論

鹽脅迫后，杜梨不同耐鹽株系基因組中轉(zhuǎn)錄因子總體DNA甲基化變化具有基因型特異性，但是兩個(gè)株系DNA差異甲基化區(qū)域均主要在基因啟動子位置，差異甲基化類型主要為mCHH。鹽脅迫下杜梨根中差異表達(dá)基因涉及243個(gè)轉(zhuǎn)錄因子，其中8個(gè)轉(zhuǎn)錄因子DNA序列的甲基化改變與轉(zhuǎn)錄變化呈負(fù)相關(guān)，可能參與調(diào)控杜梨不同株系耐鹽能力差異的分子機(jī)制。

[1] LI J G, PU L J, HAN M F, ZHU M, ZHANG R S, XIANG Y Z. Soil salinization research in China: Advances and prospects. Journal of Geographical Sciences, 2014, 24(5): 943-960.

[2] FANG S M, HOU X, LIANG X L. Response mechanisms of plants under saline-alkali stress. Frontiers in Plant Science, 2021, 12: 667458.

[3] YANG Y Q, GUO Y. Elucidating the molecular mechanisms mediating plant salt-stress responses. The New Phytologist, 2018, 217(2): 523-539.

[4] GOLLDACK D, LüKING I, YANG O. Plant tolerance to drought and salinity: Stress regulating transcription factors and their functional significance in the cellular transcriptional network. Plant Cell Reports, 2011, 30(8): 1383-1391.

[5] JIANG Y Q, DEYHOLOS M K. Comprehensive transcriptional profiling of NaCl-stressedroots reveals novel classes of responsive genes. BMC Plant Biology, 2006, 6: 25.

[6] GUO J Y, SHI G Y, GUO X Y, ZHANG L W, XU W Y, WANG Y M, SU Z, HUA J P. Transcriptome analysis reveals that distinct metabolic pathways operate in salt-tolerant and salt-sensitive upland cotton varieties subjected to salinity stress. Plant Science, 2015, 238: 33-45.

[7] GRUBER V, BLANCHET S, DIET A, ZAHAF O, BOUALEM A, KAKAR K, ALUNNI B, UDVARDI M, FRUGIER F, CRESPI M. Identification of transcription factors involved in root apex responses to salt stress in. Molecular Genetics and Genomics, 2009, 281(1): 55-66.

[8] 陳嘉貝, 張芙蓉, 黃丹楓, 張利達(dá), 張屹東. 鹽脅迫下兩個(gè)甜瓜品種轉(zhuǎn)錄因子的轉(zhuǎn)錄組分析. 植物生理學(xué)報(bào), 2014, 50(2): 150-158.

CHEN J B, ZHANG F R, HUANG D F, ZHANG L D, ZHANG Y D. Transcriptome analysis of transcription factors in two melon (L.) cultivars under salt stress. Plant Physiology Journal, 2014, 50(2): 150-158. (in Chinese)

[9] ZHANG H M, LANG Z B, ZHU J K. Dynamics and function of DNA methylation in plants. Nature Reviews Molecular Cell Biology, 2018, 19(8): 489-506.

[10] FERREIRA L J, AZEVEDO V, MAROCO J, OLIVEIRA M M, SANTOS A P. Salt tolerant and sensitive rice varieties display differential methylome flexibility under salt stress. PLoS ONE, 2015, 10(5): e0124060.

[11] WANG B H, FU R, ZHANG M, DING Z Q, CHANG L, ZHU X Y, WANG Y F, FAN B X, YE W W, YUAN Y L. Analysis of methylation-sensitive amplified polymorphism in different cotton accessions under salt stress based on capillary electrophoresis. Genes & Genomics, 2015, 37(8): 713-724.

[12] AHMAD F, FARMAN K, WASEEM M, RANA R M, NAWAZ M A, REHMAN H M, ABBAS T, BALOCH F S, AKREM A, HUANG J, ZHANG H S. Genome-wide identification, classification, expression profiling and DNA methylation (5mC) analysis of stress-responsive ZFP transcription factors in rice (L.). Gene, 2019, 718: 144018.

[13] KARAN R, DELEON T, BIRADAR H, SUBUDHI P K. Salt stress induced variation in DNA methylation pattern and its influence on gene expression in contrasting rice genotypes. PLoS ONE, 2012, 7(6): e40203.

[14] SONG Y G, JI D D, LI S, WANG P, LI Q, XIANG F N. The dynamic changes of DNA methylation and histone modifications of salt responsive transcription factor genes in soybean. PLoS ONE, 2012, 7(7): e41274.

[15] LI H, LIN J, YANG Q S, LI X G, CHANG Y H. Comprehensive analysis of differentially expressed genes under salt stress in pear () using RNA-Seq. Plant Growth Regulation, 2017, 82(3): 409-420.

[16] LI H, LIU W, YANG Q S, LIN J, CHANG Y H. Isolation and comparative analysis of two Na+/H+antiportergenes from. Plant Molecular Biology Reporter, 2018, 36(3): 439-450.

[17] WU Q S, ZOU Y N. Adaptive responses of birch-leaved pear () seedlings to salinity stress. Notulae Botanicae Horti Agrobotanici Cluj-Napoca, 2009, 37(1): 133-138.

[18] YU X J, SHI P J, HUI C, MIAO L F, LIU C L, ZHANG Q Y, FENG C N. Effects of salt stress on the leaf shape and scaling ofbunge. Symmetry, 2019, 11(8): 991.

[19] LIN L K, YUAN K L, HUANG Y D, DONG H Z, QIAO Q H, XING C H, HUANG X S, ZHANG S L. A WRKY transcription factorfromfunctions positively in salt tolerance and modulating organic acid accumulation by regulatingexpression. Environmental and Experimental Botany, 2022, 196: 104782.

[20] LIU C X, XU X Y, KAN J L, CHENG Z M, CHANG Y H, LIN J, LI H. Genome-wide analysis of the C3H zinc finger family reveals its functions in salt stress responses of. PeerJ, 2020, 8: e9328.

[21] LI H, ZHANG Y F, ZHOU X Y, LIN J, LIU C X, LI X G, CHANG Y H. Single-base resolution methylome of different ecotype fromreveals epigenomic changes in response to salt stress. Scientia Horticulturae, 2022, 306: 111437.

[22] PATEL R K, JAIN M. NGS QC Toolkit: A toolkit for quality control of next generation sequencing data. PLoS ONE, 2012, 7(2): e30619.

[23] KIM D, LANGMEAD B, SALZBERG S L. HISAT: A fast spliced aligner with low memory requirements. Nature Methods, 2015, 12(4): 357-360.

[24] DONG X G, WANG Z, TIAN L M, ZHANG Y, QI D, HUO H L, XU J Y, LI Z, LIAO R, SHI M, ALI WAHOCHO S, LIU C, ZHANG S M, TIAN Z X, CAO Y F.assembly of a wild pear() genome. Plant Biotechnology Journal, 2020, 18(2): 581-595.

[25] LOVE M I, HUBER W, ANDERS S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 2014, 15(12): 550.

[26] ZHENG Y, JIAO C, SUN H H, ROSLI H G, POMBO M A, ZHANG P F, BANF M, DAI X B, MARTIN G B, GIOVANNONI J J, ZHAO P X, RHEE S Y, FEI Z J. iTAK: A program for genome-wide prediction and classification of plant transcription factors, transcriptional regulators, and protein kinases. Molecular Plant, 2016, 9(12): 1667-1670.

[27] CHEN S F, ZHOU Y Q, CHEN Y R, GU J. Fastp: An ultra-fast all-in-one FASTQ preprocessor. Bioinformatics, 2018, 34(17): i884-i890.

[28] KRUEGER F, ANDREWS S R. Bismark: A flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics, 2011, 27(11): 1571-1572.

[29] AKALIN A, KORMAKSSON M, LI S, GARRETT-BAKELMAN F E, FIGUEROA M E, MELNICK A, MASON C E. methylKit: A comprehensive R package for the analysis of genome-wide DNA methylation profiles. Genome Biology, 2012, 13(10): R87.

[30] QUINLAN A R, HALL I M. BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics, 2010, 26(6): 841-842.

[31] LIVAK K J, SCHMITTGEN T D. Analysis of relative gene expression data using real-time quantitative PCR and the 2-ΔΔCTmethod. Methods, 2001, 25(4): 402-408.

[32] MAO H D, WANG H W, LIU S X, LI Z G, YANG X H, YAN J B, LI J S, TRAN L S P, QIN F. A transposable element in agene is associated with drought tolerance in maize seedlings. Nature Communications, 2015, 6(1): 8326.

[33] VIGGIANO L, DE PINTO M C. Dynamic DNA methylation patterns in stress response//RAJEWSKY N, JURGA S, BARCISZEWSKI J. Plant Epigenetics. Cham, Switzerland: Springer, 2017: 281-302.

[34] FINNEGAN E J, PEACOCK W J, DENNIS E S. Reduced DNA methylation inresults in abnormal plant development. Proceedings of the National Academy of Sciences of the United States of America, 1996, 93(16): 8449-8454.

[35] AL-HARRASI I, AL-YAHYAI R, YAISH M W. Differential DNA methylation and transcription profiles in date palm roots exposed to salinity. PLoS ONE, 2018, 13(1): e0191492.

[36] YAISH M W, AL-LAWATI A, AL-HARRASI I, PATANKAR H V. Genome-wide DNA Methylation analysis in response to salinity in the model plant caliph medic (). BMC Genomics, 2018, 19(1): 78.

[37] 潘教文, 李臻, 王慶國, 管延安, 李小波, 戴紹軍, 丁國華, 劉煒. NaCl處理谷子萌發(fā)期種子的轉(zhuǎn)錄組學(xué)分析. 中國農(nóng)業(yè)科學(xué), 2019, 52(22): 3964-3976.doi: 10.3864/j.issn.0578-1752.2019.22. 003.

PAN J W, LI Z, WANG Q G, GUAN Y A, LI X B, DAI S J, DING G H, LIU W. Transcriptomics analysis of Na Cl response in foxtail millet (L.) seeds at germination stage. Scientia Agricultura Sinica, 2019, 52(22): 3964-3976. doi: 10.3864/j.issn.0578-1752.2019. 22.003. (in Chinese)

[38] WANG J, MAROWSKY N C, FAN C Z. Divergence of gene body DNA methylation and evolution of plant duplicate genes. PLoS ONE, 2014, 9(10): e110357.

Identification of Salt-Tolerant Transcription Factors in the Roots ofby the Association Analysis of Genome-Wide DNA Methylation and Transcriptome

1Institute of Pomology, Jiangsu Academy of Agricultural Sciences/Jiangsu Key Laboratory for Horticultural Crop Genetic Improvement, Nanjing 210014;2College of Biology and the Environment, Nanjing Forestry University, Nanjing 210037

【】Here, two ecotypes offrom Huaguo Mountain, Lianyungang (the salt-tolerant ecotype, D) and Purple Mountain, Nanjing (the common ecotype, U) were collected for this research. The purpose of this study was to analyze the role of transcription factor genes in the roots of two ecotypes ofdiffering in terms of salt stress. Transcription factors involving in the regulation of the salt tolerance of differentecotypes were identified on the grounds of differential expression under salt stress and the relationship between the methylation status and the relative expression level of relevant tolerance genes after exposure to salt stress was investigated. 【】The 90-day-oldseedlings were grown hydroponically in Hoagland’s nutrient solution supplemented with 200 mmol?L-1NaCl, with seedlings grown in Hoagland’s nutrient solution as the control. The sodium ion content in the tissues was determined by flame graphite furnace atomic absorption spectrometry. Whole-genome DNA methylation analysis and transcriptome sequencing were performed on three replicates for the following four root samples: ecotype D and ecotype U, each grown in the presence or absence of salt stress. Bioinformatics analysis of transcription factor gene expression under salt stress at the levels of transcriptional regulation and epigenetic methylation were carried out using transcriptome sequencing data and whole-genome DNA methylation results, respectively. Then, McrBC-PCR and real-time fluorescence quantitative PCR (qPCR)were used to confirm the levels of methylation and transcription of differential transcription factor genes. 【】After exogenous NaCl treatment for 24 h, the concentration of sodium ions inroots increased significantly, with the increase in sodium ion concentration in the salt-tolerant ecotype being significantly less than that in the common ecotype. In the whole seedling, the final salt concentration of tolerant ecotype was only 73.1% of that of the common ecotype. Whereas, in the roots, the sodium content of the salt-tolerant ecotype was 1.1 times of that in the common ecotype. These results indicated that the salt-tolerant ecotype could store more sodium ions in roots and limit their upward transport after salt stress. A total of 2 682 transcription factor (TF) genes from 69 gene families were detected in roots. Among them, 243 TF genes displayed differential expression in response to salt stress, including 37 AP2/ERF, 19 bHLH, 7 bZIP, 10 HD-Zip, 30 MYB, 18 NAC, 8 WRKY, and 23 ZFP family genes. The global methylation level of transcription factor genes in the genome of the salt-tolerant rootstock ecotype decreased, whereas the overall methylation level of these genes in the common ecotype increased after exposure to 200 mmol?L-1NaCl. The differentially methylated regions in both ecotypes were mainly in the position of gene promoters, with the type of differentially methylated sequences being mostly mCHH, constituting more than 93% of the sum of all three types of methylated sequences. The expression levels of twenty-three transcription factor genes, which belonged to the AP2/ERF, bHLH, DREB, GRAS, GT factor, HB Zip, MYB, NAC, Trihelix, and zinc-finger ZFP gene families, were upregulated, and their methylation levels were downregulated in both two ecotypes in response to salt stress. These genes may be involved in the regulation of sodium uptake and accumulation in roots under salt stress. The expression patterns and promoter methylation of representative candidate genes identified by bioinformatics analysis were confirmed by qPCR and McrBC-qPCR.【】The differentially expressed genes in roots ofunder salt stress included 243 transcription factor genes in both ecotypes. The methylation changes in DNA sequences in eight transcription factor genes (,,,,,,, and) were correlated with their transcriptional activity. Our results provided preliminary experimental evidence for supporting a relationship between promoter DNA methylation and expression of TF genes inin response to salt stress as part of the molecular role of TFs involved in the regulation of salt tolerance among differentecotypes, which would increase our understanding of the role of epigenetics in the response of woody trees to abiotic stress.

;salt stress; transcription factors; DNA methylation; transcriptome

2022-05-05；

2023-02-01

江蘇省自然科學(xué)基金（BK20191238）、江蘇現(xiàn)代農(nóng)業(yè)（梨）產(chǎn)業(yè)技術(shù)體系項(xiàng)目[ATS（2021）436]、國家自然科學(xué)基金面上項(xiàng)目（31772287）

通信作者李慧，E-mail：lihui7904@163.com

（責(zé)任編輯 趙伶俐）