刺參響應(yīng)燦爛弧菌侵染的基因組DNA甲基化水平和轉(zhuǎn)錄組差異及其關(guān)聯(lián)分析*

李欣容 廖梅杰 李 彬 榮小軍 暢孟陽 王印庚 于永翔 張 正 范瑞用 劉清兵

刺參響應(yīng)燦爛弧菌侵染的基因組DNA甲基化水平和轉(zhuǎn)錄組差異及其關(guān)聯(lián)分析*

李欣容1,2廖梅杰2,3①李 彬2,3榮小軍2,3暢孟陽2,3王印庚2,3于永翔2,3張 正2,3范瑞用4劉清兵4

(1. 上海海洋大學(xué)水產(chǎn)與生命學(xué)院 上海 201306；2. 中國水產(chǎn)科學(xué)研究院黃海水產(chǎn)研究所 農(nóng)業(yè)農(nóng)村部海洋漁業(yè)可持續(xù)發(fā)展重點(diǎn)實(shí)驗(yàn)室 山東 青島 266071；3. 青島海洋科學(xué)與技術(shù)試點(diǎn)國家實(shí)驗(yàn)室海洋漁業(yè)科學(xué)與食物產(chǎn)出過程功能實(shí)驗(yàn)室 山東 青島 266071；4. 青島瑞滋集團(tuán)有限公司 山東 青島 266408)

為探討病原菌脅迫下刺參()基因組DNA甲基化水平和轉(zhuǎn)錄組表達(dá)的差異，本研究采用人工攻毒侵染脅迫，獲得刺參化皮體壁組織，并以未攻毒組的健康體壁組織為對照，對刺參2種體壁組織進(jìn)行全基因組甲基化測序(WGBS)和轉(zhuǎn)錄組高通量測序，解析刺參體壁基因組DNA甲基化差異，篩選響應(yīng)病原脅迫的差異甲基化區(qū)域和差異表達(dá)基因。同時(shí)，通過基因組甲基化和轉(zhuǎn)錄組聯(lián)合分析，篩選負(fù)相關(guān)關(guān)聯(lián)基因，為解析刺參響應(yīng)燦爛弧菌()侵染的分子機(jī)理提供基礎(chǔ)數(shù)據(jù)。結(jié)果顯示，刺參對照組和侵染組基因組總甲基化水平分別為(3.59±0.04)%和(3.87±0.27)%，侵染組甲基化水平顯著升高；在發(fā)生甲基化的位點(diǎn)中，攻毒組和對照組的mCpG占比分別為83.06%和81.91%，mCpG為最主要的甲基化形式。本研究共篩選出差異甲基化區(qū)域(DMRs) 626677個(gè)，注釋到23706個(gè)功能基因。轉(zhuǎn)錄組測序共檢測到29290個(gè)基因，篩選到差異表達(dá)基因496個(gè)，其中，上調(diào)基因214個(gè)，下調(diào)基因282個(gè)。在基因組甲基化與轉(zhuǎn)錄組的關(guān)聯(lián)分析中，篩選到180個(gè)負(fù)相關(guān)關(guān)聯(lián)基因，其中，差異甲基化區(qū)域位于啟動(dòng)子區(qū)域的負(fù)相關(guān)基因?yàn)?0個(gè)。對負(fù)相關(guān)關(guān)聯(lián)基因的GO和KEGG富集分析，篩選到相關(guān)通路和、和等關(guān)鍵基因。本研究將為解析刺參抗病的表觀遺傳調(diào)控機(jī)制提供數(shù)據(jù)，也為刺參抗病性狀選育提供科學(xué)參考。

刺參；DNA甲基化；轉(zhuǎn)錄組；關(guān)聯(lián)分析

仿刺參又稱刺參()，具有重要的營養(yǎng)價(jià)值和藥用價(jià)值，是我國第五次海水養(yǎng)殖浪潮的代表性物種之一，其養(yǎng)殖業(yè)為我國沿海經(jīng)濟(jì)結(jié)構(gòu)調(diào)整和漁民增收開辟了一條新途徑(王印庚等, 2014)。在人工養(yǎng)殖刺參過程中，刺參的病害問題日益凸顯，給刺參養(yǎng)殖業(yè)造成了嚴(yán)重的經(jīng)濟(jì)損失，嚴(yán)重制約了該產(chǎn)業(yè)的穩(wěn)定發(fā)展。開展刺參抗病機(jī)制解析對刺參良種選育和病害防控具有重要意義。DNA甲基化是表觀遺傳修飾的主要方式，在不改變DNA序列的前提下，改變遺傳表現(xiàn)在調(diào)控基因表達(dá)和染色質(zhì)構(gòu)象等方面發(fā)揮著重要作用(Moore, 2013)。研究表明，刺參在腸道組織再生、夏眠、化皮及倍性等過程或現(xiàn)象中，相關(guān)組織均有顯著的基因組甲基化水平差異變化，表明甲基化在刺參相關(guān)生理過程中發(fā)揮重要的調(diào)控作用(Yang, 2020; Zhao, 2015; Sun, 2020; Han, 2021; 高杉等, 2017; 李曉妮, 2018; 趙業(yè), 2015)。本研究以刺參腐皮綜合征重要致病原——燦爛弧菌()為病原對刺參苗種進(jìn)行人工攻毒侵染實(shí)驗(yàn)，以攻毒后發(fā)病個(gè)體和未攻毒的健康個(gè)體的體壁組織為樣本，采用全基因組甲基化測序技術(shù)(WGBS)和轉(zhuǎn)錄組高通量測序技術(shù)，分析燦爛弧菌脅迫下刺參體壁DNA甲基化水平變化及基因表達(dá)差異變化，并進(jìn)一步對甲基化位點(diǎn)和轉(zhuǎn)錄組差異進(jìn)行關(guān)聯(lián)分析，篩選出關(guān)鍵位點(diǎn)及關(guān)聯(lián)基因，為解析刺參響應(yīng)燦爛弧菌侵染的分子機(jī)制提供基礎(chǔ)數(shù)據(jù)，也有助于為刺參抗病性狀選育提供參考。

1 材料與方法

1.1 實(shí)驗(yàn)材料

本研究所用健康苗種取自山東青島瑞滋集團(tuán)有限公司，選取活力良好、健康的個(gè)體，苗種規(guī)格為(50.0±2.0) g/只。苗種運(yùn)回實(shí)驗(yàn)室后暫養(yǎng)3 d，暫養(yǎng)水溫為(13.0±0.5)℃，待苗種狀態(tài)穩(wěn)定后用于后續(xù)實(shí)驗(yàn)。

攻毒用菌株為本團(tuán)隊(duì)病原庫中保存的分離自患腐皮綜合征的刺參病樣的燦爛弧菌菌株(AJ-Vb1801)。對該菌株用胰蛋白胨大豆肉湯培養(yǎng)基(TSB)固體培養(yǎng)基復(fù)蘇，然后用TSB液體培養(yǎng)基擴(kuò)大培養(yǎng)。

1.2 實(shí)驗(yàn)樣品的采集

燦爛弧菌人工攻毒侵染與取樣：根據(jù)實(shí)驗(yàn)?zāi)康膶?shí)驗(yàn)分為2組，每組設(shè)3個(gè)平行，每個(gè)平行實(shí)驗(yàn)苗種數(shù)量為30頭，實(shí)驗(yàn)水槽容積為50 L。對照組(PT10H)為在水質(zhì)良好的自然海水中養(yǎng)殖的健康苗種，侵染組(PT16S)按照水體體積在養(yǎng)殖水槽中添加培養(yǎng)的燦爛弧菌菌液至終濃度1×106CFU/mL (該濃度為燦爛弧菌對刺參苗種的半致死濃度)。養(yǎng)殖條件：溫度(13.0±0.5)℃，鹽度(28.0±0.5)，每天換1/3水，換水后及時(shí)添加燦爛弧菌菌液，使其維持在1×106CFU/mL。每天換水后定時(shí)投喂刺參配合飼料，投喂量為刺參苗種體重的2%，養(yǎng)殖期間持續(xù)充氣。實(shí)驗(yàn)期間，每天觀察刺參苗種的生理狀態(tài)，及時(shí)收集侵染組發(fā)生化皮的瀕死個(gè)體，分別標(biāo)記為PT16S1、PT16S2和PT16S3，在相同時(shí)間節(jié)點(diǎn)采集對照組健康個(gè)體，分別標(biāo)記為PT10H1、PT10H2和PT10H3。剖取所采集樣品的體壁組織，置于2 mL凍存管后迅速放入液氮罐中，送回實(shí)驗(yàn)室，–80℃保存，用于后期DNA和RNA的提取。

1.3 DNA和RNA的提取

以對照組和侵染組的體壁組織為材料，分別采用Omega公司的Mollusc DNA kit和QIAGEN RNeasy mini kit提取各樣品體壁基因組DNA和RNA，1.5%瓊脂糖凝膠電泳檢測DNA和RNA完整性，NanoDrop核酸測定儀測定所提取核酸樣品的純度，檢測合格的樣品于–80℃保存?zhèn)溆谩?/p>

1.4 全基因組DNA甲基化測序

將所提取的6個(gè)DNA樣品送至杭州聯(lián)川生物技術(shù)股份有限公司分別構(gòu)建甲基化文庫，并利用HiSeq4000測序平臺(tái)150PE上機(jī)策略進(jìn)行測序。

1.5 數(shù)據(jù)處理與比對

測序完成后對原始數(shù)據(jù)進(jìn)行預(yù)處理，去除可能含有測序接頭序列和低質(zhì)量的測序數(shù)據(jù)(含有N的比例大于5%，以及質(zhì)量值Q≤10的堿基數(shù)占整個(gè)read的20%以上的reads)，得到有效數(shù)據(jù)，用于后續(xù)的生物信息學(xué)分析。

以本團(tuán)隊(duì)組裝的刺參全基因組序列(數(shù)據(jù)待發(fā)表)為基礎(chǔ)，對預(yù)處理后的有效測序數(shù)據(jù)使用Bismark (Krueger, 2011)進(jìn)行參考基因組序列比對，統(tǒng)計(jì)基因組比對結(jié)果。

1.6 刺參體壁基因組DNA甲基化分布特征分析

根據(jù)基因組比對分析結(jié)果，統(tǒng)計(jì)所測樣品各類型甲基化C位點(diǎn)(mCpG、mCHG和mCHH)的數(shù)目，及其在全部mC的位點(diǎn)中所占的比例，解析物種的全基因組DNA甲基化修飾特征。

1.7 差異甲基化區(qū)域分析

使用R語言程序包中的methylKit篩選差異甲基化區(qū)域(differentially methylated regions, DMRs)，默認(rèn)選擇1000 bp windows、500 bp overlap、FDR<0.05為差異篩選閾值，進(jìn)行DMR分析。統(tǒng)計(jì)差異表達(dá)分析中侵染組的甲基化水平顯著性升高或下降的基因數(shù)目。對篩選到的具有DMR的基因進(jìn)行注釋，按照基因結(jié)構(gòu)繪制柱狀統(tǒng)計(jì)圖。

1.8 轉(zhuǎn)錄組中差異表達(dá)基因的篩選

將所采集的6個(gè)樣品的RNA送至杭州聯(lián)川生物技術(shù)股份有限公司分別構(gòu)建cDNA文庫，采用Illumina NovaseqTM6000進(jìn)行測序，測序讀長為雙端2×150 bp (PE150)。使用String Tie軟件計(jì)算不同樣品間各基因表達(dá)定量的結(jié)果，以FPKM (fragments per kilobase of exon model per million mapped reads)為單位對已知基因統(tǒng)計(jì)其在不同樣本中的表達(dá)豐度，篩選差異表達(dá)基因(different expression genes, DEGs)。

1.9 DNA甲基化與轉(zhuǎn)錄組聯(lián)合分析

對轉(zhuǎn)錄組測序和基因甲基化測序數(shù)據(jù)進(jìn)行聯(lián)合分析，統(tǒng)計(jì)DMRs和DEGs注釋的共有基因，檢測這些共有基因的DMR甲基化水平與DEG表達(dá)水平之間的相關(guān)性，統(tǒng)計(jì)DMR甲基化水平與DEG表達(dá)水平存在負(fù)相關(guān)的基因，分別采用GOseq和KOBAS軟件對篩選到的差異甲基化基因進(jìn)行GO和KEGG pathway富集分析，采用ggplot2軟件繪制差異甲基化基因的KEGG散點(diǎn)圖。

2 結(jié)果

2.1 甲基化測序結(jié)果統(tǒng)計(jì)

對實(shí)驗(yàn)獲得的對照組(PT10H)和侵染組(PT16S)的6個(gè)樣品的WGBS測序和數(shù)據(jù)統(tǒng)計(jì)分析結(jié)果見表1。高通量測序平均每個(gè)樣品產(chǎn)出165510349條的原始數(shù)據(jù)，對原始數(shù)據(jù)進(jìn)行預(yù)處理后，截去測序數(shù)據(jù)的測序接頭和去除低質(zhì)量的數(shù)據(jù)，得到均值為23.20 Gb的有效數(shù)據(jù)，有效序列占比在99.18%以上，Q30以上的數(shù)據(jù)占比大于91.96%，有37.31%～39.81%的mapped reads可比對到刺參的參考基因組，說明測序數(shù)據(jù)的質(zhì)量較好。

2.2 基因組中甲基化位點(diǎn)統(tǒng)計(jì)

根據(jù)基因組比對分析結(jié)果，表2統(tǒng)計(jì)了對照組和侵染組刺參體壁組織中C位點(diǎn)以及不同類型(CG、CHG和CHH)的平均甲基化水平。侵染組樣品基因組甲基化C位點(diǎn)占全部C位點(diǎn)數(shù)的(3.87±0.27)%，對照組為(3.59±0.04)%，說明在病原菌侵染脅迫下，基因組甲基化C位點(diǎn)占比顯著升高(<0.05)。不同類型(CpG、CHG和CHH)的平均甲基化水平統(tǒng)計(jì)結(jié)果顯示，對照組mCpG、mCHG和mCHH占CpG、CHG和CHH的比例分別為24.14%、0.56%和0.47%，而侵染組的mCpG、mCHG和mCHH占CpG、CHG和CHH的比例分別為24.22%、0.98%和0.76%。由此可見，對于刺參基因組來說，CpG位點(diǎn)發(fā)生甲基化的比例顯著高于CHG和CHH。從對照組和侵染組不同類型位點(diǎn)甲基化對比結(jié)果可以看出，侵染組的mCHG和mCHH比例顯著高于對照組。

表1 WGBS測序數(shù)據(jù)統(tǒng)計(jì)

Tab.1 Statistics of the WGBS sequencing data

表2 不同類型C位點(diǎn)甲基化水平

Tab.2 Methylation level of different cytosine contexts

注：*代表有顯著差異(<0.05)。下同

Note: * represents a significant difference (<0.05). The same as below

對照組和侵染組基因組DNA不同形式甲基化C位點(diǎn)占總甲基化C位點(diǎn)的比例見圖1。如圖1所示，在總的mC位點(diǎn)中，侵染組和對照組mCpG占比分別為83.06%和81.91%，mCHH占比分別為13.41%和14.28%，mCHG占比分別為3.53%和3.81%，說明甲基化位點(diǎn)主要集中在CpG二核苷酸序列的胞嘧啶上，CpG甲基化是主要的甲基化形式。

圖1 對照組(PT10H)和侵染組(PT16S) mCpG、 mCHG和mCHH占總mC的百分比

2.3 差異甲基化區(qū)域分析

使用R語言程序包中的methylKit分析自對照組和侵染組樣品測序結(jié)果中篩選到的626677個(gè)DMRs。根據(jù)分組的差異比較結(jié)果，對發(fā)生顯著甲基化的基因數(shù)量的統(tǒng)計(jì)結(jié)果見圖2。從對照組和侵染組共篩選到的啟動(dòng)子區(qū)域DMRs共有40329個(gè)(6.44%)，其中，甲基化水平升高的有22613個(gè)，下降的有17716個(gè)；外顯子區(qū)域存在58522個(gè)DMRs (9.34%)，其中，升高的有29731個(gè)，下降的有28791個(gè)；內(nèi)含子區(qū)域DMRs共有109635個(gè)(17.49%)，其中，升高的有56603個(gè)，下降的有53032個(gè)；基因間區(qū)DMRs共有401174個(gè)(64.02%)，其中，升高的有216152個(gè)，下降的有185022個(gè)；CGI.CG島DMRs共有17017個(gè)(2.72%)，其中，升高的有8877個(gè)，下降的有8140個(gè)。刺參基因間區(qū)的差異甲基化基因比例最高，CGI.CG島區(qū)域的差異甲基化比例最低。

圖2 不同基因功能區(qū)域的DMR數(shù)量

2.4 轉(zhuǎn)錄組顯著性差異基因表達(dá)分析

由對照組和侵染組轉(zhuǎn)錄組測序分析共檢測出29290個(gè)基因，對所檢測的基因表達(dá)進(jìn)行差異顯著性分析，共篩選到496個(gè)顯著性差異表達(dá)基因(DEGs)，其中，上調(diào)性基因214個(gè)，下調(diào)性基因有282個(gè)(圖3)。從整體上來看，下調(diào)基因數(shù)目多于上調(diào)基因數(shù)目。

圖3 PT16S與PT10H差異表達(dá)基因火山圖

2.5 基因組甲基化與轉(zhuǎn)錄組關(guān)聯(lián)分析

對DMRs功能基因注釋篩選到的23706個(gè)功能基因和DEGs篩選到的496個(gè)功能基因進(jìn)行比對，篩選到DMRs和DEGs注釋的共有基因有261個(gè)(圖4)。對這261個(gè)基因的DMRs差異和基因表達(dá)差異進(jìn)行關(guān)聯(lián)分析，結(jié)果見表3。共篩選到180個(gè)負(fù)相關(guān)關(guān)聯(lián)基因，其中，甲基化水平高而基因表達(dá)下調(diào)的負(fù)相關(guān)關(guān)聯(lián)基因有58個(gè)，甲基化水平低而基因表達(dá)上調(diào)的負(fù)相關(guān)關(guān)聯(lián)基因有122個(gè)。對所篩選的負(fù)相關(guān)關(guān)聯(lián)基因差異甲基化區(qū)域進(jìn)行統(tǒng)計(jì)，差異甲基化區(qū)域位于啟動(dòng)子區(qū)域的負(fù)相關(guān)基因有60個(gè)。

圖4 差異表達(dá)基因和差異甲基化區(qū)域共注釋基因的韋恩圖

對所篩選的甲基化和轉(zhuǎn)錄組負(fù)相關(guān)的180個(gè)基因進(jìn)行GO富集分析，結(jié)果顯示(圖5)，甲基化水平與表達(dá)水平呈負(fù)相關(guān)的180個(gè)基因中，有124個(gè)基因注釋到896個(gè)GO terms，其中，50個(gè)GO terms顯著富集(<0.05)。富集到生物過程的GO terms主要包括細(xì)胞粘附、DNA轉(zhuǎn)錄模板、先天性免疫應(yīng)答、氧化還原過程、RNA聚合酶Ⅱ?qū)D(zhuǎn)錄的負(fù)調(diào)控等；富集到細(xì)胞成分的GO terms主要包括膜的組成部分、細(xì)胞質(zhì)、細(xì)胞膜、細(xì)胞核等；富集到分子功能的GO terms主要包括蛋白結(jié)合、金屬離子結(jié)合、鈣離子結(jié)合、ATP結(jié)合、鋅離子結(jié)合、DNA結(jié)合、轉(zhuǎn)運(yùn)蛋白活性等。篩選到的功能基因包括免疫球蛋白、富含亮氨酸的重復(fù)序列(LRR)、有機(jī)陰離子轉(zhuǎn)運(yùn)蛋白多肽(OATP)、半胱天冬酶募集結(jié)構(gòu)域(CARD)、黃素結(jié)合單加氧酶(almA)、CD36、血紅素結(jié)合蛋白等。

利用KEGG數(shù)據(jù)庫對差異表達(dá)基因進(jìn)行pathway分析，結(jié)果見圖6，甲基化水平與基因表達(dá)水平呈負(fù)相關(guān)的有180個(gè)基因，其中，39個(gè)富集基因注釋到112個(gè)KEGG信號(hào)通路中，有20個(gè)通路顯著富集(<0.05)，內(nèi)質(zhì)網(wǎng)中的蛋白質(zhì)加工、破骨細(xì)胞分化、乙型肝炎、維生素消化吸收這4個(gè)通路注釋基因最多，均有3個(gè)負(fù)相關(guān)關(guān)聯(lián)基因富集在相應(yīng)通路上，篩選到的功能基因包括20和等。

3 討論

3.1 甲基化檢測和基因組測序?qū)z測效果的影響

DNA甲基化是表觀遺傳學(xué)研究的重要內(nèi)容，是基因組DNA的一種重要修飾方式，參與生物體或細(xì)胞的生物調(diào)控過程(Moore, 2013)。水生生物在響應(yīng)溫度、病原、飼料、倍性及雜種優(yōu)勢過程中DNA甲基化均發(fā)生了改變(吳彪等, 2016; 高杉等, 2017; 卓梅琴等, 2019; Han, 2021; 李炎璐等, 2019)。

目前，水產(chǎn)動(dòng)物研究中常用的DNA甲基化檢測方法有基于限制性酶切預(yù)處理的甲基化敏感擴(kuò)增多態(tài)性技術(shù)(MASP)和基于高通量測序的亞硫酸氫鹽處理的全基因組DNA甲基化測序技術(shù)(WGBS)。MASP技術(shù)具有簡便、高效和可靠的特點(diǎn)。研究人員利用該技術(shù)探究了刺參不同群體(左之良等, 2016; 左閃等, 2021)、不同組織(郭婷婷, 2013)、高溫脅迫(李尚俊等, 2017)、夏眠(趙業(yè), 2015)和化皮(高杉等, 2017)的甲基化水平差異。但MASP無法完全準(zhǔn)確地顯示基因組中甲基化情況的全貌及具體的甲基化位點(diǎn)的序列等信息。

表3 基因組甲基化DMRs與轉(zhuǎn)錄組DEGs關(guān)聯(lián)分析篩選基因統(tǒng)計(jì)

Tab.3 Statistics of gene number detected by association analysis between DMRs and DEGs

圖5 甲基化與基因表達(dá)負(fù)相關(guān)基因GO富集

圖6 甲基化與基因表達(dá)負(fù)相關(guān)基因KEGG富集

隨著高通量測序技術(shù)的發(fā)展，WGBS成為檢測DNA甲基化的最可靠技術(shù)之一，該方法可以直觀地反映全基因組DNA甲基化的詳細(xì)信息，可以直接提供篩選出的位點(diǎn)及其對應(yīng)的信息序列。Yang等(2020)運(yùn)用該技術(shù)篩選到夏眠期刺參的DMR區(qū)域和關(guān)鍵功能基因；溫爭爭等(2021)利用WGBS技術(shù)完成了刺參高溫脅迫下基因組DNA甲基化水平及模式的變化研究；Sun等(2020)利用該技術(shù)篩選到了高濃度燦爛弧菌侵染下刺參化皮組織DMR區(qū)域。全基因組序列信息是開展WBGS分析的關(guān)鍵，參考基因組拼接的完整性對DMR區(qū)域的篩選和注釋效果產(chǎn)生直接影響。Zhang等(2017)和Li等(2018)獲得NCBI數(shù)據(jù)庫中的刺參基因組序列信息，共包括3821個(gè)Scaffolds (N50為786 kb)。Sun等(2020)利用高濃度燦爛弧菌(5×109CFU/mL)攻毒制備化皮實(shí)驗(yàn)材料，利用NCBI中的全基因組數(shù)據(jù)比對篩選得到116522個(gè)DMRs。本團(tuán)隊(duì)于2020年采用納米孔測序技術(shù)和染色體構(gòu)象捕獲技術(shù)(Hi-C)開展了刺參基因組精細(xì)圖譜繪制，拼接獲得911 Mb的染色體水平的刺參基因組，共定位到23條染色體，N50為39.61 Mb (相關(guān)數(shù)據(jù)未發(fā)表)，基因組拼接質(zhì)量得到較大提升。

本研究利用燦爛弧菌半致死濃度(1×106CFU/mL)制備化皮實(shí)驗(yàn)材料，利用染色體水平基因組序列篩選得到的DMRs數(shù)目為626677個(gè)，顯著高于Sun等(2020)篩選到的116522個(gè)DMRs，所篩選的DMRs的差異數(shù)除了與實(shí)驗(yàn)條件差異有關(guān)外，對比時(shí)所選用的基因組拼接質(zhì)量是產(chǎn)生此差異的主要原因之一。

3.2 刺參基因組甲基化特征及其響應(yīng)病原脅迫的變化

高通量測序獲得的所有樣品的基因組甲基化水平檢測結(jié)果表明，本研究所測定的刺參體壁組織甲基化水平為3.52～4.17%。Yang等(2020)測定的夏眠期刺參基因組甲基化水平在3.5%左右；Sun等(2020)檢測刺參甲基化水平在4%左右?？梢钥闯觯虆⒒蚪M整體甲基化水平較低，這與Tweedie等(1997)發(fā)現(xiàn)的無脊椎動(dòng)物基因組甲基化水平通常低于脊椎動(dòng)物的結(jié)果相一致。

對刺參的不同類型C位點(diǎn)甲基化水平的統(tǒng)計(jì)結(jié)果可以看出，CG類型位點(diǎn)的甲基化比例較高，而CHG和CHH位點(diǎn)甲基化比例很低，與目前已報(bào)道的其他海洋無脊椎動(dòng)物CG類型甲基化水平范圍(20%～ 50%)相一致(曹哲明等, 2009; 于濤等, 2010; Sun,2014)。進(jìn)一步對發(fā)生甲基化C位點(diǎn)的分型比例統(tǒng)計(jì)結(jié)果可以看出，mCpG占甲基化位點(diǎn)的比例達(dá)到80%以上，表明mCpG是最主要的甲基化形式。與李玉強(qiáng)等(2018)利用methylRAD-seq技術(shù)對仿刺參的甲基化圖譜研究結(jié)果相一致。這些結(jié)果表明海洋動(dòng)物的DNA甲基化模式都有一定的相似性和保守性。

對病原侵染刺激下基因組甲基化水平的變化可以看出，侵染組樣本的甲基化水平顯著高于對照組。雖然CG位點(diǎn)甲基化比例未發(fā)生顯著變化，但CHH和CHG的甲基化比例顯著升高，證明病原脅迫會(huì)引起基因組甲基化水平的升高。高杉等(2017)、Yang等(2020)、Sun等(2020)和溫爭爭等(2021)研究結(jié)果也表明，溫度和高濃度病原刺激等外界脅迫會(huì)引起基因組甲基化水平的升高，尤其是CHH和CHG的甲基化比例的升高。雖然mCpG是刺參最主要的甲基化形式，但在不同刺激影響下的占比仍存在差異。本研究在病原半致死濃度刺激下，mCpG的占比在82%左右。Sun等(2020)測定的在致死濃度病原刺激下，mCpG的占比為87%～89%，但在溫度刺激或夏眠狀態(tài)下，mCpG的占比達(dá)到94%以上。相關(guān)對比分析證明，不同的刺激形式造成的基因組甲基化存在差異。

3.3 刺參響應(yīng)燦爛弧菌侵染差異甲基化位點(diǎn)和差異表達(dá)基因分析

本研究利用染色體水平基因組為基礎(chǔ)，篩選得到626677個(gè)DMRs，且在侵染組中篩選到的高甲基化DMR數(shù)量大于低甲基化DMR數(shù)量，這與侵染組全基因組甲基化水平增高相一致。根據(jù)這些DMR所在的功能區(qū)域的位置統(tǒng)計(jì)結(jié)果可以看出，基因間區(qū)的差異甲基化位點(diǎn)數(shù)量最高，這與李玉強(qiáng)(2018)利用methylRAD-seq技術(shù)分析結(jié)果解釋的甲基化標(biāo)簽較大比例分布于基因間區(qū)的結(jié)果相一致。對篩選的DMR區(qū)域進(jìn)行功能基因注釋得到23706個(gè)功能基因，所注釋的功能基因數(shù)目也遠(yuǎn)高于Sun等(2020)，這可能也與實(shí)驗(yàn)條件的差異和比對所選用的基因組序列不同有關(guān)。轉(zhuǎn)錄組測序結(jié)果表明，侵染組篩選到的下調(diào)基因數(shù)目多于上調(diào)基因數(shù)目，這個(gè)趨勢與侵染組篩選到高甲基化位點(diǎn)數(shù)目較多正好相反，也從側(cè)面驗(yàn)證了Choy等(2010)提出的甲基化水平的升高會(huì)抑制基因表達(dá)，二者存在負(fù)調(diào)控關(guān)系。在所篩選到的496個(gè)差異表達(dá)基因中發(fā)現(xiàn)DNA甲基化轉(zhuǎn)移酶的表達(dá)顯著上調(diào)，相關(guān)結(jié)果與Yang等(2020)和李尚俊等(2017)檢測到的相關(guān)基因上調(diào)表達(dá)相一致，也解釋了該基因上調(diào)表達(dá)與基因組甲基化水平的升高直接相關(guān)。

3.4 刺參響應(yīng)病原脅迫的基因組甲基化與基因表達(dá)的關(guān)聯(lián)分析

基因組DNA甲基化作為基因的表達(dá)調(diào)控的重要方式，在通常情況下，DNA低甲基化會(huì)激活基因的轉(zhuǎn)錄，而高甲基化則抑制基因的轉(zhuǎn)錄(Bird, 2002)。如早期的海膽()甲基化研究報(bào)道了DNA甲基化對于海膽發(fā)育具有重要意義(Tosi, 1995)。本研究自甲基化和轉(zhuǎn)錄組關(guān)聯(lián)分析中篩選到261個(gè)共同注釋的基因，其中180個(gè)為負(fù)相關(guān)基因，占比為68.97%，說明所篩選到的基因其DNA甲基化對基因表達(dá)發(fā)揮了重要作用。這些負(fù)相關(guān)基因是后期解析刺參表觀調(diào)控需要重點(diǎn)關(guān)注的對象。

本研究對所篩選的180個(gè)負(fù)相關(guān)基因進(jìn)行GO和KEGG分析。在甲基化和轉(zhuǎn)錄組負(fù)相關(guān)的180個(gè)基因的GO富集分析中，篩選到生物進(jìn)程、細(xì)胞成分和分子功能多個(gè)GO terms中的富含亮氨酸的重復(fù)序列基因，該基因包含許多抗病基因，能參與固有免疫(Gou, 2010)。富集到的CARD蛋白參與先天性免疫應(yīng)答、細(xì)胞質(zhì)、ATP結(jié)合、乙型肝炎、RIG-I樣受體信號(hào)通路、麻疹、甲型流感等通路，進(jìn)而參與包括炎癥反應(yīng)、免疫、分化、細(xì)胞生長、腫瘤發(fā)生和凋亡等許多生物過程(Beg, 1994)。此外，篩選到的免疫球蛋白能夠避免免疫系統(tǒng)的過分激活，降低炎癥反應(yīng)。雖然在KEGG通路中富集的負(fù)相關(guān)基因數(shù)量僅為39個(gè)，但篩選到的和基因參與免疫調(diào)控、細(xì)胞粘附等多種生理和病理過程(Li, 2015; 魏書磊等, 2013)。這些重要的功能基因?qū)⒊蔀榻馕龃虆㈨憫?yīng)病原侵染分子機(jī)理的重要關(guān)注點(diǎn)。

綜上所述，本研究以染色體水平的基因組序列為基礎(chǔ)，通過對病原侵染后化皮和健康樣品體壁組織的WGBS和轉(zhuǎn)錄組測序分析，篩選到相應(yīng)的差異甲基化位點(diǎn)和差異表達(dá)基因，進(jìn)一步通過關(guān)聯(lián)分析篩選到刺參響應(yīng)病原侵染的關(guān)鍵基因及其甲基化特征，相關(guān)結(jié)果將為解析刺參抗病調(diào)控機(jī)理奠定基礎(chǔ)。

BEG A A, BALDWIN A S. Activation of multiple NF-kappa B/Rel DNA-binding complexes by tumor necrosis factor. Oncogene, 1994, 9(5): 1487–1492

BIRD A. DNA methylation patterns and epigenetic memory. Genes and Development, 2002, 16(1): 6–21

CAO Z M, YANG J. Analysis of the methylation in genome DNA from different tissues of. Ecology and Environmental Sciences, 2009, 18(6): 2011–2016 [曹哲明, 楊健. 背角無齒蚌不同組織的基因組DNA甲基化分析. 生態(tài)環(huán)境學(xué)報(bào), 2009, 18(6): 2011–2016]

CHOY M K, MOVASSAGH M, GOH H,. Genome-wide conserved consensus transcription factor binding motifs are hyper-methylated. BMC Genomics, 2010, 11(1): 519

GAO S, YANG A F, DONG Y,. MSAP analysis of DNA methylation in the body wall of. Acta Hydrobiologica Sinica, 2017, 41(3): 637–642 [高杉, 楊愛馥, 董穎, 等. 仿刺參“化皮病”體壁組織DNA甲基化的MSAP分析. 水生生物學(xué)報(bào), 2017, 41(3): 637–642]

GOU X P, HE K, YANG H,. Genome-wide cloning and sequence analysis of leucine-rich repeat receptor-like protein kinase genes in. BMC Genomics, 2010, 11: 19

GUO T T. The method to assess genome DNA methylation ofby HPLC and MSAP analysis in different tissues. Master′s Thesis of Shanghai Ocean University, 2013 [郭婷婷. 刺參不同組織基因組DNA甲基化狀態(tài)MSAP分析及HPLC方法的建立. 上海海洋大學(xué)碩士研究生學(xué)位論文, 2013]

HAN L S, SUN Y, CAO Y,. Analysis of the gene transcription patterns and DNA methylation characteristics of triploid sea cucumbers (). Scientific Reports, 2021, 11(1): 7564

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

LI J K, WU X W, TAN J,. Molecular cloning of the heat shock protein 20 gene from Paphia textile and its expression in response to heat shock. Chinese Journal of Oceanology and Limnology, 2015, 33(4): 919–927

LI S J, SUN G H, LI X Y,. Characteristics of epigenetic regulation of related genes under high temperature stress in sea cucumber. Journal of Fishery Sciences of China, 2017, 24(3): 470–476 [李尚俊, 孫國華, 李雪燕, 等. 高溫脅迫下仿刺參表觀遺傳調(diào)控相關(guān)基因的表達(dá)特征. 中國水產(chǎn)科學(xué), 2017, 24(3): 470–476]

LI X N. Analysis of DNA methylation regulation and function of candidated genes during intestinal regeneration in the sea cucumber. Doctoral Dissertation of Institute of Oceanology, Chinese Academy of Sciences, University of Chinese Academy of Sciences, 2018 [李曉妮. 刺參腸道再生的DNA甲基化調(diào)控解析及再生候選基因的功能分析. 中國科學(xué)院大學(xué)(中國科學(xué)院海洋研究所)博士研究生學(xué)位論文, 2018]

LI Y L, CHEN C, CHEN J G,. DNA methylation analysis of,and their F1hybrid. Progress in Fishery Sciences, 2019, 40(6): 98–104 [李炎璐, 陳超, 陳建國, 等. 云紋石斑魚、鞍帶石斑魚及其雜交F1的DNA甲基化分析. 漁業(yè)科學(xué)進(jìn)展, 2019, 40(6): 98–104]

LI Y Q, WANG R J, LI Y L,. Genome-wide profiling of DNA methylation inbased on methylRAD-Seq. Periodical of Ocean University of China (Natural Science), 2018, 48(9): 41–50 [李玉強(qiáng), 王睿甲, 李語麗, 等. 基于MethylRAD-Seq技術(shù)對仿刺參DNA甲基化圖譜的研究. 中國海洋大學(xué)學(xué)報(bào)(自然科學(xué)版), 2018, 48(9): 41–50]

LI Y, WANG R, XUN X,. Sea cucumber genome provides insights into saponin biosynthesis and aestivation regulation. Cell Discovery, 2018, 4(1): 29

MOORE L D, LE T, FAN G P. DNA methylation and its basic function. Neuropsychopharmacology, 2013, 38(1): 23–38

SUN H J, ZHOU Z C, DONG Y,. Insights into the DNA methylation of sea cucumberin response to skin ulceration syndrome infection. Fish and Shellfish Immunology, 2020, 104: 155–164

SUN Y, HOU R, FU X,Genome-wide analysis of DNA methylation in five tissues of Zhikong scallop,. PloS One, 2014, 9(1): e86232

TOSI L, ANIELLO F, GERACI G,. DNA methyltransferase activity in the early stages of a sea urchin embryo-evidence of differential control. FEBS Letters, 1995, 361(1): 115–117

TWEEDIE S, CHARLTON J, CLARK V,. Methylation of genomes and genes at the invertebrate-vertebrate boundary. Molecular and Cellular Biology, 1997, 17(3): 1469–1475

WANG Y G, RONG X J, LIAO M J,. Sea cucumber culture and disease control technology. Beijing: China Agriculture Press, 2014 [王印庚, 榮小軍, 廖梅杰, 等. 刺參健康養(yǎng)殖與病害防控技術(shù)叢解. 北京: 中國農(nóng)業(yè)出版社, 2014]

WEI S L. Expression and functional analysis of CD36 gene in zebrafish. Master′s Thesis of China Ocean University, 2013 [魏書磊. 斑馬魚()重組CD36基因的原核表達(dá)及功能研究. 中國海洋大學(xué)碩士研究生學(xué)位論文, 2013]

WEN Z Z, ZUO S, CHEN M,. DNA methylation level of genomic DNA ofat different temperatures. Progress in Fishery Sciences, 2021, 42(3): 46–54 [溫爭爭, 左閃, 陳夢, 等. 刺參基因組DNA甲基化水平及模式對溫度變化的響應(yīng). 漁業(yè)科學(xué)進(jìn)展, 2021, 42(3): 46–54]

WU B, YANG A G, SUN X J,. Effects of acute temperature stress on genome-wide DNA methylation profiles inProgress in Fishery Sciences, 2016, 37(5): 140–146 [吳彪, 楊愛國, 孫秀俊, 等. 急性溫度脅迫對蝦夷扇貝()基因組DNA甲基化的影響. 漁業(yè)科學(xué)進(jìn)展, 2016, 37(5): 140–146]

YANG Y J, ZHENG Y Q, SUN L N,. Genome-wide DNA methylation signatures of sea cucumberduring environmental induced aestivation. Genes, 2020, 11(9): 1020

YU T, YANG A G, WU B,. Analysis of,and their offspring using methylation-sensitive amplification polymorphism (MSAP). Journal of Fisheries of China, 2010, 34(9): 1335–1342 [于濤, 楊愛國, 吳彪, 等. 櫛孔扇貝、蝦夷扇貝及其雜交子代的MSAP分析. 水產(chǎn)學(xué)報(bào), 2010, 34(9): 1335–1342]

ZHANG X J, SUN L, YUAN J B,. The sea cucumber genome provides insights into morphological evolution and visceral regeneration. PLoS Biology, 2017, 15(10): e2003790

ZHAO Y, CHEN M Y, STOREY K B,. DNA methylation levels analysis in four tissues of sea cucumberbased on fluorescence-labeled methylation- sensitive amplified polymorphism (F-MSAP) during aestivation. Comparative Biochemistry and Physiology Part B: Biochemistry and Molecular Biology, 2015, 181: 26–32

ZHAO Y. Study on gene expression patterns and basic characteristics of DNA methylation in(Selenka) during aestivation. Doctoral Dissertation of Institute of Oceanology, Chinese Academy of Sciences, 2015 [趙業(yè). 刺參(Selenka)夏眠期間基因表達(dá)模式及DNA甲基化基礎(chǔ)特征研究. 中國科學(xué)院研究生院(海洋研究所)博士研究生學(xué)位論文, 2015]

ZHUO M Q, YANG S B, LING S C,. Effects of dietary lipid on lipid metabolism, methylation and expression of PI3KCa in the ovary of yellow catfish (). Journal of Fisheries of China, 2019, 43(10): 2186–2196 [卓梅琴, 楊水波, 凌仕誠, 等. 飼料脂肪對黃顙魚卵巢脂類代謝以及PI3KCa甲基化和轉(zhuǎn)錄水平的影響. 水產(chǎn)學(xué)報(bào), 2019, 43(10): 2186–2196]

ZUO S, WEN Z Z, ZHOU H X,. Evaluation of epigenetic and genome sequence diversity in sea cucumberselected population based on MSAP technology. Progress in Fishery Sciences, 2021, 42(3): 38–45 [左閃, 溫爭爭, 周紅學(xué), 等. 基于MSAP技術(shù)的刺參選育群體基因組表觀與序列遺傳多樣性分析. 漁業(yè)科學(xué)進(jìn)展, 2021, 42(3): 38–45]

ZUO Z L, TAN J, WU B,. MSAP analysis of genomic DNA in the tissues ofand whiteProgress in Fishery Sciences, 2016, 37(3): 93–100 [左之良, 譚杰, 吳彪, 等. 普通刺參()和白刺參不同組織基因組DNA的MSAP研究. 漁業(yè)科學(xué)進(jìn)展, 2016, 37(3): 93–100]

Genomic DNA Methylation Levels and Transcriptome Differences ofin Response toInfection and Their Association Analysis

LI Xinrong1,2, LIAO Meijie2,3①, LI Bin2,3, RONG Xiaojun2,3, CHANG Mengyang2,3, WANG Yingeng2,3, YU Yongxiang2,3, ZHANG Zheng2,3, FAN Ruiyong4, LIU Qingbing4

(1. College of Fisheries and Life Science, Shanghai Ocean University, Shanghai 201306, China; 2. Key Laboratory of Sustainable and Development of Marine Fisheries, Ministry of Agriculture and Rural Affairs, Yellow Sea Fisheries Research Institute, Chinese Academy of Fishery Sciences, Qingdao, Shandong 266071, China;3. Pilot National Laboratory for Marine Science and Technology (Qingdao), Laboratory for Marine Fisheries Science and Food Production Processes, Qingdao, Shandong 266071, China;4. Qingdao Ruizi Company, Qingdao, Shandong 266408, China)

DNA methylation is an important epigenetic modification that plays a key role in gene expression regulation. In this study, two groups of sea cucumbers () were prepared. One group had skin ulceration syndrome body wall (PT16S) under stress frominfection at a concentration of 1×106CFU/mL (LD50); the other group had a healthy body wall (PT10H). Genomic DNA methylation levels and gene expression differences between the two groups were detected using whole-genome bisulfite sequencing (WGBS) and transcriptome sequencing. The key gene ontology (GO) terms and KEGG terms engaged in the immune response were selected using association analysis. The results showed that the total methylation levels of thegenome of PT10H and PT16S were (3.59±0.04)% and (3.87±0.27)%, respectively. The methylation levels of thegenome under pathogen challenge significantly increased; mCpG accounted for 83.06% and 81.91% of all the methylated sites in PT16S and PT10H, respectively, indicating that mCpG was the most important methylation form in the sea cucumber. A total of 626677 differentially methylated regions (DMRs) were screened and annotated into 23706 functional genes. A total of 496 differentially expressed genes were screened, including 214 up-regulated and 282 down-regulated genes. A total of 180 negatively correlated genes were isolated using association analysis between genomic methylation and transcriptome, of which 60 genes had DMRs located in promoter regions. Based on GO and KEGG enrichment of the 180 negatively correlated genes, key genes such as,, andwere selected to play a critical role in the immune response. The results would provide primary data for the epigenetic regulatory mechanisms ofand provide a theoretical reference forbreeding.

; DNA methylation; Transcriptome; Association analysis

S968.9

A

2095-9869(2022)03-0176-10

10.19663/j.issn2095-9869.20210424001

http://www.yykxjz.cn/

李欣容, 廖梅杰, 李彬, 榮小軍, 暢孟陽, 王印庚, 于永翔, 張正, 范瑞用, 劉清兵. 刺參響應(yīng)燦爛弧菌侵染的基因組DNA甲基化水平和轉(zhuǎn)錄組差異及其關(guān)聯(lián)分析. 漁業(yè)科學(xué)進(jìn)展, 2022, 43(3): 176–185

LI X R, LIAO M J, LI B, RONG X J, CHANG M Y, WANG Y G, YU Y X, ZHANG Z, FAN R Y, LIU Q B. Genomic DNA methylation levels and transcriptome differences ofin response toinfection and their association analysis. Progress in Fishery Sciences, 2022, 43(3): 176–185

LIAO Meijie, E-mail: liaomj@ysfri.ac.cn

* 國家重點(diǎn)研發(fā)計(jì)劃(2018YFD0901603)、山東省農(nóng)業(yè)良種工程課題(2020LZGC015)、中國水產(chǎn)科學(xué)研究院中央級(jí)公益性科研院所基本科研業(yè)務(wù)費(fèi)專項(xiàng)資金(2020TD40; 2021GH05)共同資助 [This work was supported by National Key R&D Program of China (2018YFD0901603), Agriculture Seed Improvement Project of Shandong Province (2020LZGC015), and Central Public-Interest Scientific Institution Basal Research Fund, CAFS (2020TD40; 2021GH05)]. 李欣容，E-mail: 1123295662@qq.com

廖梅杰，研究員，E-mail: liaomj@ysfri.ac.cn

2021-04-24,

2021-04-27

(編輯 馮小花)