鎘脅迫下桑樹幼苗的轉(zhuǎn)錄組分析

王曉紅，韓世玉，孫運(yùn)朋，張芳，羅澤虎，陳彪

（貴州省農(nóng)業(yè)科學(xué)院蠶業(yè)研究所，貴陽 550025）

0 引言

來自自然、農(nóng)業(yè)和工業(yè)來源的重金屬對耕地的污染已成為世界性的公共衛(wèi)生問題[1]。鎘(Cd)是一種重要的重金屬污染物，被認(rèn)為是最具毒性的重金屬之一。鎘污染土壤的修復(fù)是當(dāng)前農(nóng)業(yè)面臨的一個(gè)重要問題。用于這一目的的不同的物理和化學(xué)方法都存在成本高、勞動(dòng)強(qiáng)度大、土壤性質(zhì)改變和土壤原生微生物區(qū)系干擾等嚴(yán)重的局限性。相比之下，植物修復(fù)技術(shù)是一種環(huán)境友好、經(jīng)濟(jì)高效、人們接受度較高的技術(shù)[2]。植物可以提取或穩(wěn)定土壤中的Cd，在污染土壤修復(fù)中得到廣泛應(yīng)用，一些Cd超富集植物可以在不影響其正常生長的情況下從污染土壤中吸收和提取Cd。然而，多數(shù)Cd超富集植物是草本植物，其生物量小、根系淺、生長緩慢。相比之下，人們逐漸提出使用具有大生物量和多分枝生根系統(tǒng)的木本植物來修復(fù)Cd污染土壤[3]。

桑樹是一個(gè)具有較高經(jīng)濟(jì)和生態(tài)價(jià)值的木本植物，具有適應(yīng)性強(qiáng)、生長快、根系深、地上生物量大、耐剪伐等優(yōu)點(diǎn)，應(yīng)用于石漠化、鹽堿地、戈壁沙漠等多種植物修復(fù)系統(tǒng)中[4]。同時(shí)，桑樹還被證明能夠耐受Cd、Cr、Pb、Co、Cu等重金屬，具有修復(fù)重金屬污染土壤的極好潛力[5-7]。此外，利用桑樹進(jìn)行Cd污染土壤修復(fù)也是一種經(jīng)濟(jì)可行的做法，因?yàn)樯Ｈ~是家蠶的飼料，可用于蠶業(yè)發(fā)展。目前，多數(shù)研究集中在桑樹對Cd脅迫生理反應(yīng)或部分耐Cd 基因的功能研究[8-9]。全面闡明桑樹中Cd 吸收、積累、轉(zhuǎn)運(yùn)和解毒的機(jī)制對進(jìn)一步提高桑樹在Cd污染土壤中的植物修復(fù)價(jià)值具有重要意義。

Cd 通過抑制光合作用、呼吸和氮代謝，擾亂水和礦物質(zhì)養(yǎng)分的吸收和分配，阻礙植物的生長和發(fā)育[10-11]。為減輕Cd 的毒害作用，植物已經(jīng)進(jìn)化出一系列的解毒策略，包括螯合作用、細(xì)胞壁固定、區(qū)室化作用和抗氧化系統(tǒng)解除重金屬毒害[12]。這些解毒機(jī)制的闡明需要廣泛的轉(zhuǎn)錄、蛋白質(zhì)組學(xué)和代謝研究。轉(zhuǎn)錄組測序技術(shù)是從分子水平闡明其機(jī)制的有效方法，該技術(shù)具有效率高、通量大、運(yùn)行時(shí)間快、準(zhǔn)確性高的特點(diǎn)，常用于發(fā)現(xiàn)和表征基因、分析功能基因變異、識別和量化轉(zhuǎn)錄本[13]。利用這種方法已鑒定出柳樹[3]、白菜[14]和蘿卜[15]等的Cd脅迫相關(guān)基因。此外，比較轉(zhuǎn)錄組學(xué)分析已用于揭示植物中Cd積累的機(jī)制[14,16]。這些研究極大地提高了人們對Cd耐受和積累的分子調(diào)控機(jī)制的理解。本研究將桑樹幼苗在Cd 處理12、24、48 h后，采用Illumina HiSeq 2000平臺進(jìn)行測序及分析，旨在了解桑樹對Cd的吸收、轉(zhuǎn)運(yùn)、隔離和耐受機(jī)制，并在植物修復(fù)研究中發(fā)揮作用。

1 材料與方法

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

試驗(yàn)于2021 年3—8 月在貴州省蠶業(yè)研究所植物生理實(shí)驗(yàn)室進(jìn)行。試驗(yàn)桑樹品種為‘桂桑優(yōu)12’，購自廣西蠶業(yè)推廣總站。挑選完整、飽滿的種子，75%酒精滅菌60 s，無菌水沖洗5 次。將滅菌后的桑樹種子接種到MS 培養(yǎng)基上，16 h/8 h 光周期、24℃條件下萌發(fā)培養(yǎng)。待胚根長至2 cm左右時(shí)取出幼苗，接種到添加50 μmol/L CdCl2的MS培養(yǎng)基中，每瓶接種10棵幼苗；同時(shí)，以未添加CdCl2的幼苗為對照。分別在處理12、24、48 h取樣，迅速用無菌水將幼苗沖洗干凈后放入液氮中速凍，保存于-80℃?zhèn)溆谩?/p>

1.2 RNA提取、文庫構(gòu)建及轉(zhuǎn)錄組測序

分別收集Cd 處理和對照組12、24、48 h 的桑樹幼苗提取RNA用于轉(zhuǎn)錄組測序?？俁NA提取、cDNA文庫構(gòu)建和轉(zhuǎn)錄組測序由北京百邁克生物公司完成。主要流程如下：用帶有Oligo(dT)的磁珠富集真核生物mRNA；加入Fragmentation Buffer 將mRNA 進(jìn)行隨機(jī)打斷；以mRNA 為模板，用六堿基隨機(jī)引物合成第一條cDNA 鏈，然后加入緩沖液、dNTPs、RNase H、DNA polymerase I 合成第二條cDNA 鏈，利用AMPure XP beads純化cDNA；純化的雙鏈cDNA再進(jìn)行末端修復(fù)、加A 尾并連接測序接頭，然后用AMPure XP beads 進(jìn)行片段大小選擇；最后通過PCR 富集得到cDNA 文庫。對檢測合格的文庫進(jìn)行測序。

1.3 轉(zhuǎn)錄組從頭組裝和功能注釋

fastq格式的原始讀取首先通過內(nèi)部的perl腳本進(jìn)行處理，從原始數(shù)據(jù)中刪除包含接頭的序列，超過10%未知核苷酸的序列以及低質(zhì)量的序列，以獲得干凈序列。使用Trinity軟件對干凈讀數(shù)進(jìn)行組裝。隨后，基于數(shù)據(jù)庫NR(NCBI non-redundant protein sequences)、Swiss-Prot、GO(Gene Ontology)、KOG/COG/eggNOG(Clusters of Orthologous Groups of proteins)和KEGG(Kyoto Encyclopedia of Genes and genome)，使用BLASTx對轉(zhuǎn)錄本進(jìn)行注釋，期望值E-value≤10-5。使用HMMER 將預(yù)測的氨基酸序列與Pfam 數(shù)據(jù)庫進(jìn)行比較，期望值E-value≤10-10，獲得注釋信息。

1.4 Unigenes的差異表達(dá)分析

為了確定Cd 脅迫下的差異表達(dá)基因，用Bowtie軟件將干凈序列映射到組裝的轉(zhuǎn)錄組上。使用RSEM估計(jì)基因表達(dá)水平，unigenes 的表達(dá)水平用FPKM 歸一化。使用DEGSeq 鑒定Cd 處理樣本和對照組之間的差異基因。2 組之間的錯(cuò)誤發(fā)現(xiàn)率(FDR)≤0.01 和|log2 Fold Change|≥1 被作為差異表達(dá)的閾值。隨后，使用topGO 和KOBAS 分別對差異表達(dá)的基因進(jìn)行GO富集分析和KEGG通路富集分析。

1.5 qRT-PCR分析

采用實(shí)時(shí)熒光定量PCR(qRT-PCR)技術(shù)對轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行驗(yàn)證，隨機(jī)選擇4個(gè)DEGs進(jìn)行qRT-PCR。以桑樹MaACT基因作為內(nèi)參基因[17]，根據(jù)各基因序列設(shè)計(jì)引物（表1），在Bio-Rad實(shí)時(shí)熒光定量PCR儀上檢測差異基因表達(dá)水平。根據(jù)2×RealStar Green Fast Mixture（GenStar，北京）試劑盒的說明書配制qRT-PCR反應(yīng)體系：2×RealStar Green Fast Mixture 10 μL、模板cDNA 1 μL、正/反引物各0.5 μL(10 μmol/L)，ddH2O 8 μL。PCR程序：94℃預(yù)變性2 min；94℃變性15 s，60℃退火30 s，循環(huán)40次。采用2ΔΔCt方法計(jì)算基因相對表達(dá)量。

表1 研究所用qRT-PCR引物

2 結(jié)果與分析

2.1 桑樹幼苗轉(zhuǎn)錄組測序、組裝與注釋結(jié)果統(tǒng)計(jì)

為全面分析Cd脅迫下桑樹幼苗的轉(zhuǎn)錄組，將對照組T1N (12 h)、T2N (24 h)、T3N (48 h)和Cd 處理組T1P(12 h)、T2P(24 h)、T3P(48 h)的6個(gè)樣品進(jìn)行轉(zhuǎn)錄組測序，6個(gè)樣本分別產(chǎn)生了37398420～49460056條原始序列，去除低質(zhì)量序列后，獲得了36139366～47208418 條干凈序列，占原始序列的94.99%以上，Q30在92.15%～93.44%之間（表2），說明轉(zhuǎn)錄組測序質(zhì)量較高，可以進(jìn)行下一步的數(shù)據(jù)整理和分析。使用Trinity軟件對干凈序列進(jìn)行從頭組裝，結(jié)果得到39758個(gè)unigenes，長度在201～13312 bp之間。這些unigenes的平均長度和N50 分別為1134.51、1968 bp，平均GC含量為41.72%（表3）。

表2 Illumina轉(zhuǎn)錄組測序總結(jié)

表3 非冗余轉(zhuǎn)錄本統(tǒng)計(jì)

2.2 Unigenes功能注釋與分類

組裝后，通過查詢7個(gè)公共數(shù)據(jù)庫預(yù)測基因功能，39758 個(gè)unigenes 得到注釋（表4）。其中，有26280 個(gè)unigenes(66.10%)在Nr數(shù)據(jù)庫中獲得了注釋，18589個(gè)unigenes(46.76%)在Swiss-prot 數(shù)據(jù)庫中獲得了注釋，18606個(gè)unigenes(46.8%)在GO數(shù)據(jù)庫中獲得了注釋，只有7201個(gè)unigenes(18.11%)在KEGG數(shù)據(jù)庫中獲得了注釋。

表4 桑樹unigenes在公共數(shù)據(jù)庫中的注釋結(jié)果

在KOG分析中，有13734個(gè)單基因被分成25組用于功能預(yù)測和分類（圖1）。其中，僅一般功能預(yù)測(3831)包含的基因數(shù)量最多，其次是翻譯后修飾、蛋白質(zhì)轉(zhuǎn)換、伴侶(1465)、信號轉(zhuǎn)導(dǎo)機(jī)制(1344)、碳水化合物運(yùn)輸和代謝(833)和轉(zhuǎn)錄(762)。此外，18606 個(gè)單基因被注釋為一個(gè)或多個(gè)GO 術(shù)語，并被分為3 個(gè)類別和57 個(gè)組（圖2）。對于生物過程(biological process，BP)類別，細(xì)胞過程(12445)代表最大的組，其次是代謝過程(11271)、生物調(diào)控(5578)和刺激反應(yīng)(4620)。在分子功能(molecular function，MF)類別中，結(jié)合(12335)、催化活性(10218)、轉(zhuǎn)運(yùn)活性(1229)和轉(zhuǎn)錄調(diào)節(jié)活性(1190)是最常見的組。在細(xì)胞成分(cellular component，CC)類別中，前5 名分別是細(xì)胞組分(14509)、細(xì)胞器(9053)、膜(5491)、細(xì)胞器的部分(4674)和膜部分(4125)。

圖2 桑樹轉(zhuǎn)錄本GO注釋分類

2.3 差異表達(dá)基因分析

分別對Cd 脅迫和對照處理12 h (T1N_vs_T1P)、24 h (T2N_vs_T2P)、48 h (T3N_vs_T3P) 樣本的unigenes 進(jìn)行差異表達(dá)分析（圖3）。在T1N_vs_T1P組，共有6319 個(gè)差異表達(dá)基因，其中，4411 個(gè)基因上調(diào)表達(dá)，1908 個(gè)基因下調(diào)表達(dá)；7737 個(gè)基因在T2N_vs_T2P 中存在差異，其中，2129 個(gè)基因上調(diào)表達(dá)，5608 個(gè)基因下調(diào)表達(dá)；而在T3N_vs_T3P 中，共有6561 個(gè)差異表達(dá)基因，其中，3389 個(gè)基因上調(diào)表達(dá)，3172 個(gè)基因下調(diào)表達(dá)（圖3A），表明桑樹幼苗在Cd 脅迫不同時(shí)期轉(zhuǎn)錄水平的動(dòng)態(tài)變化過程。值得注意的是，有534 個(gè)基因在不同時(shí)期均有差異（圖3B），其中有4 個(gè)基因(TRINITY_DN10246_c0_g4、TRINITY_DN9070_c0_g1、TRINITY_DN12871_c0_g1、TRINITY_DN7347_c1_g2)注釋為ABC 轉(zhuǎn)運(yùn)蛋白家族成員，其在Cd 脅迫過程中上調(diào)表達(dá)，表明它們可能在整個(gè)Cd脅迫過程中發(fā)揮作用。

圖3 桑樹各處理組差異表達(dá)基因

2.4 差異表達(dá)基因GO富集分析

為了進(jìn)一步研究差異表達(dá)基因在Cd 脅迫下的生物學(xué)功能，進(jìn)行了GO 富集分析。Cd 脅迫12、24、48 h后，每組DEGs可分為生物過程、細(xì)胞成分和分子功能3類，選取富集顯著程度最高的10個(gè)GO條目分析（圖4）。在生物過程類別，Cd 脅迫不同時(shí)間后，差異表達(dá)基因主要集中于能量代謝、DNA代謝和應(yīng)激反應(yīng)相關(guān)的組中，如ATP代謝過程(ATP metabolic process)、ADP代謝過程(ADP metabolic process)、DNA 重組(DNA recombination)、DNA 代謝過程(DNA metabolic process)和急性炎癥反應(yīng)(acute inflammatory response)等；在細(xì)胞組分類別，MHC 蛋白復(fù)合體(MHC protein complex)、NADPH 脫氫酶復(fù)合體[NAD(P)H dehydrogenase complex (plastoquinone)]、細(xì)胞(cell body)和細(xì)胞壁(cell wall)等是幾個(gè)富集的組；在分子功能類別，主要富集的組與能量轉(zhuǎn)換、DNA 復(fù)制與轉(zhuǎn)錄過程酶活性有關(guān)，例如，NADH 脫氫酶活性(NADH dehydrogenase activity)、RNA 聚合酶II 轉(zhuǎn)錄因子活性(RNA polymerase II transcription factor activity)、DNA結(jié)合轉(zhuǎn)錄因子活性(DNA binding transcription factor activity)和DNA 聚合酶活性(DNA polymerase activity)等；其他富集的組還包括?；d體蛋白質(zhì)合成酶活性(3-oxoacyl-[acyl-carrier-protein] synthase activity)和7-脫氧馬錢酸葡糖基轉(zhuǎn)移酶活性(7-deoxyloganetic acid glucosyltransferase activity)等。這些GO 富集結(jié)果表明，Cd脅迫引起桑樹幼苗能量代謝、酶催化活性，以及DNA復(fù)制和轉(zhuǎn)錄相關(guān)途徑基因的差異表達(dá)，這些過程可能在桑樹幼苗抵抗Cd脅迫中起到重要的作用。

2.5 差異表達(dá)基因KEGG通路富集分析

KEGG數(shù)據(jù)庫是系統(tǒng)分析基因產(chǎn)物在細(xì)胞中的代謝途徑以及這些基因產(chǎn)物功能的數(shù)據(jù)庫。以q＜0.05為標(biāo)準(zhǔn)，對KEGG中每個(gè)通路應(yīng)用超幾何檢驗(yàn)進(jìn)行富集分析，找出差異表達(dá)基因中顯著性富集的通路。結(jié)果如圖5 所示。在Cd 脅迫不同時(shí)間，谷胱甘肽代謝(glutathione metabolism)、倍半萜和三萜生物合成(sesquiterpenoid and triterpenoid biosynthesis)、光合作用-天線蛋白(photosynthesis-antenna proteins)、苯丙素的生物合成(phenylpropanoid biosynthesis)、光合作用(photosynthesis)、戊糖和葡萄糖醛酸轉(zhuǎn)換(pentose and glucuronate interconversions)、異黃酮生物合成(isoflavonoid biosynthesis)、半乳糖代謝(galactose metabolism)、氧化磷酸化(oxidative phosphorylation)、類黃酮生物合成(flavonoid biosynthesis)和DNA 復(fù)制(DNA replication)等11 條代謝通路差異表達(dá)基因富集程度較高，這些代謝通路可能與桑樹中Cd的代謝相關(guān)。

圖5 差異表達(dá)基因KEGG富集分析

2.6 差異轉(zhuǎn)錄因子分析

許多類型的轉(zhuǎn)錄因子廣泛參與非生物和生物脅迫反應(yīng)，并在各種植物發(fā)育過程和脅迫反應(yīng)中發(fā)揮重要的調(diào)節(jié)作用。與對照相比，在Cd 脅迫不同時(shí)期，桑樹幼苗中共鑒定了148 個(gè)差異表達(dá)的轉(zhuǎn)錄因子，代表了36個(gè)不同的家族（圖6）。其中，bHLH、MYB、B3、NAC、MYB-related和WRKY家族是擁有差異基因較多的轉(zhuǎn)錄因子家族。

圖6 差異表達(dá)轉(zhuǎn)錄因子分類

2.7 qRT-PCR驗(yàn)證轉(zhuǎn)錄組數(shù)據(jù)

為進(jìn)一步驗(yàn)證轉(zhuǎn)錄組的比較結(jié)果，隨機(jī)選擇4 個(gè)差異表達(dá)基因進(jìn)行qRT-PCR 檢測（圖7），這些基因的表達(dá)模式與轉(zhuǎn)錄組檢測到的基因表達(dá)變化是一致的，表明轉(zhuǎn)錄組譜數(shù)據(jù)是可靠的。

3 結(jié)論與討論

Cd是一種常見的環(huán)境污染物，對人類健康構(gòu)成潛在的威脅。如何有效地治理Cd 污染，降低Cd 對生物體的潛在風(fēng)險(xiǎn)，已成為一個(gè)重要的社會(huì)問題[18]。雖然已經(jīng)開展了各種研究來探討植物對Cd反應(yīng)的生理、遺傳和分子基礎(chǔ)，但桑樹中Cd 吸收、積累和運(yùn)輸?shù)姆肿诱{(diào)控機(jī)制在很大程度上仍未被認(rèn)識。作為一種重要的生態(tài)經(jīng)濟(jì)樹種，桑樹被認(rèn)為是解析重金屬調(diào)控網(wǎng)絡(luò)的最重要的木本植物之一[9,19]。

最近，RNA-Seq方法已被廣泛應(yīng)用于許多重要的植物的從頭轉(zhuǎn)錄組測序組裝，如蘿卜[15]、玉米[20]、柳樹[3]和水稻[21]。在本研究中，通過從頭轉(zhuǎn)錄組測序和組裝，從Cd 脅迫下的桑樹幼苗中獲得了39758 個(gè)非冗余轉(zhuǎn)錄本?；赑fam、Nr、GO 和KEGG 數(shù)據(jù)庫，通過BLAST和功能生物信息學(xué)分析，成功地注釋了全部轉(zhuǎn)錄本序列。本研究中的轉(zhuǎn)錄組數(shù)據(jù)將為進(jìn)一步研究桑樹中Cd脅迫響應(yīng)基因和功能提供有價(jià)值的信息。

轉(zhuǎn)錄組測序分析使人們能夠全面了解桑樹響應(yīng)Cd 脅迫有關(guān)的基因和過程。在Cd 脅迫下，桑樹共檢測到15882 個(gè)差異表達(dá)基因。其中，分別有6319、7737、6561 個(gè)差異表達(dá)基因在脅迫12、24、48 h 被鑒定，表明Cd脅迫導(dǎo)致桑樹大量基因進(jìn)行了重編程。在這些差異表達(dá)基因中有4 個(gè)轉(zhuǎn)錄本被注釋為ABC 轉(zhuǎn)運(yùn)蛋白家族成員，它們可能對桑樹Cd耐受性具有重要作用。ABC轉(zhuǎn)運(yùn)蛋白廣泛存在于動(dòng)植物中，負(fù)責(zé)各類底物在生物膜上的運(yùn)輸，包括重金屬、激素和羧酸鹽等[22]。已有的報(bào)道表明，ABC轉(zhuǎn)運(yùn)蛋白對Cd吸收轉(zhuǎn)運(yùn)具有重要的調(diào)節(jié)作用，AtABCC3缺陷型擬南芥幼苗對不同Cd 濃度的敏感性增加，過度表達(dá)AtABCC3基因的幼苗對Cd的耐受性增加[23]。

Cd 通過多種毒害作用抑制植物生長。已有研究表明，Cd會(huì)損害光合機(jī)構(gòu)，改變光合過程，最終導(dǎo)致植物生長抑制[24]。Cd 還會(huì)改變氧化酶的活性來誘導(dǎo)氧化應(yīng)激，從而對植物產(chǎn)生不利影響[25]。相反，植物會(huì)通過增強(qiáng)防御反應(yīng)，如金屬螯合、液泡隔離、轉(zhuǎn)運(yùn)蛋白調(diào)節(jié)以及強(qiáng)化抗氧化機(jī)制來抵抗重金屬威脅[26]。本研究中，KEGG富集分析表明，具有差異表達(dá)的基因在谷胱甘肽代謝、次生代謝、光合作用、能量代謝和DNA復(fù)制等相關(guān)代謝通路富集程度較高。這些結(jié)果與構(gòu)樹類似，在Cd 脅迫下，構(gòu)樹差異表達(dá)基因在多個(gè)KEGG 途徑富集，包括苯丙烷生物合成、氧化磷酸化和類黃酮生物合成等[16]。表明植物在抵抗Cd 脅迫時(shí)具有相似的機(jī)制。

已有的研究表明，許多轉(zhuǎn)錄因子家族在Cd脅迫中起著重要作用。Hu 等[27]報(bào)道水稻的一個(gè)轉(zhuǎn)錄因子OsMYB45參與水稻的Cd脅迫反應(yīng)，OsMYB45突變導(dǎo)致對Cd 處理的超敏反應(yīng)，而OsMYB45過度表達(dá)補(bǔ)充了突變體表型。Yang等[28]報(bào)道ThWRKY7可能通過激活ThVHAc1基因表達(dá)來提高轉(zhuǎn)基因擬南芥對Cd的抗性。此外，蘿卜中bZIP、MYB、ERF、C2H2、NAC 家族被證明在介導(dǎo)Cd 脅迫反應(yīng)中起著重要的調(diào)節(jié)作用[15]。bHLH、GRAS、WRKY、bZIP 家族在構(gòu)樹Cd 脅迫中具有復(fù)雜的轉(zhuǎn)錄調(diào)控作用[16]。在本研究中，大量轉(zhuǎn)錄因子在Cd 脅迫下差異表達(dá)，包括bHLH、MYB、NAC、WRKY等轉(zhuǎn)錄因子家族。表明不同轉(zhuǎn)錄因子家族的調(diào)控可能會(huì)觸發(fā)植物對Cd 脅迫的響應(yīng)，也說明Cd信號調(diào)控的復(fù)雜性。

總之，本研究從桑樹幼苗轉(zhuǎn)錄組中分離到39758個(gè)非冗余轉(zhuǎn)錄本，其中15882個(gè)基因在Cd脅迫下受到差異調(diào)控。KEGG 富集分析表明，這些差異表達(dá)基因參與了許多代謝、生物合成和非生物脅迫響應(yīng)途徑的過程。同時(shí)，還鑒定了一些ABC轉(zhuǎn)運(yùn)蛋白和大量轉(zhuǎn)錄因子，這些發(fā)現(xiàn)將為研究桑樹Cd脅迫反應(yīng)的分子機(jī)制提供有價(jià)值的信息。