鎘脅迫后旱柳轉(zhuǎn)錄組變化分析

曹繼敏 李雙財 何德

（西南林業(yè)大學(xué)生命科學(xué)學(xué)院，昆明 650224）

隨著工業(yè)的發(fā)展，鎘（Cd）成為全球污染最主要的重金屬之一，如工業(yè)廢氣會將Cd帶到氣流層中，并以大氣沉降、降雨以及下雪等方式進入土壤中，造成土壤Cd污染；同時，農(nóng)業(yè)上的污水灌溉以及含Cd藥劑的使用都會使Cd進入土壤，從而加重土壤的Cd污染。此外土壤中的Cd主要通過根的吸收進入植物體內(nèi)，然后與根細胞的細胞壁結(jié)合，后經(jīng)木質(zhì)部運輸?shù)角o、葉、果實等器官中積累，抑制植物的生長發(fā)育［1］。Cd可以通過食物鏈進入動物和人體內(nèi)，對機體造成嚴(yán)重傷害。過量Cd積累導(dǎo)致肺紋理增多、紊亂而模糊，白細胞和中性粒細胞數(shù)量增多，引起肺泡Ⅱ型上皮細胞毛細血管內(nèi)皮細胞受損，使肺泡間隔增厚，間質(zhì)纖維增生，也會直接作用于骨骼，使有機體骨骼礦化發(fā)生障礙、骨鈣溶出的增加以及骨膠原和骨的固化作用異常等［2-3］。20世紀(jì)60年代發(fā)生在日本神通川流域的“骨痛病”，原因就是當(dāng)?shù)鼐用袷秤酶缓亟饘貱d的大米造成的。在我國也發(fā)生了由于土壤Cd污染日益嚴(yán)重出現(xiàn)了例如“鎘大米”，農(nóng)作物減產(chǎn)等一系列惡性事件［4］。

考慮到傳統(tǒng)的物理、化學(xué)修復(fù)重金屬污染的方法中存在不可避免的缺陷［5］，植物修復(fù)以其成本低、安全可靠、對環(huán)境干擾小且改善生態(tài)環(huán)境等優(yōu)點而備受關(guān)注［6］，木本植物因為其具有龐大的根系、巨大的生物量、發(fā)達的維管組織、堅固的木質(zhì)組織、以及高效的蒸騰作用，并且可以有吸收多種重金屬的特點而成為環(huán)境重金屬污染治理的重要方向［7］。旱柳（Saliz matsudana）作為一種常見的木本植物，具有生長快、易繁殖、生物量大、根系發(fā)達、再生能力強和對多種重金屬耐受的特點［8］，此外旱柳的葉與根等新陳代謝較快的器官可以大量積累重金屬Cd，而其它的主要營養(yǎng)儲存器官（果實、籽粒、塊莖等）對于重金屬積累較少，因此旱柳具有巨大的治理重金屬污染的潛力［9］。

當(dāng)前對旱柳重金屬耐受性的研究主要集中在對重金屬的植物體內(nèi)富集情況的調(diào)查、影響因素、生理變化等方面，對重金屬富集后植物耐受的分子機理研究較少，尤其在轉(zhuǎn)錄水平上研究很少。Claudia等［10］對柳樹耐Cd無性系葉片中Cd的沉積研究發(fā)現(xiàn)Cd主要分布于嫩葉的葉尖以及老葉的葉頸，水培時間較短的旱柳主要以嫩葉為主。目前我國土壤一般Cd污染濃度為2.5 mg/L，徐愛春［11］研究表明Cd濃度為50 mg/L時對旱柳具有較大影響效應(yīng)，本研究以旱柳葉片為實驗材料從轉(zhuǎn)錄組水平著手，探究旱柳在Cd脅迫條件下的耐受性反應(yīng)情況，為旱柳修復(fù)土壤Cd污染提供理論指導(dǎo)。

1 材料與方法

1.1 材料

本研究所用的旱柳枝條取樣自昆明市北京路延長線的一株旱柳。一年生旱柳枝條被截取長度為25 cm繼而用改良后霍格蘭德營養(yǎng)液［12］進行水培。水培后的旱柳隨機分成3組（每組5盆，每盆3棵），1組為對照組，其余2組為試驗組。Cd（NO3）2作為唯一的Cd來源對試驗組進行脅迫處理，使?fàn)I養(yǎng)液中Cd濃度分別為2.5 mg/L和 50 mg/L，營養(yǎng)液每3 d換一次。研究材料為3組的第0、1、7與30天并處于同一位置的葉片，采樣方法使用隨機混采法；所有試驗葉片采集完成后液氮急速冷凍后放置于-80℃冰箱內(nèi)保存?zhèn)溆谩φ战M用CK表示，2.5 mg/L處理組用L表示，50 mg/L處理組用H表示；樣品編號如下:第0天為 CK0，第1 天分別為CK1、L1、H1，第7天分別為CK7、L7、H7，第30天分別為CK30、L30、H30。

1.2 方法

1.2.1 總RNA提取及樣品檢測 總RNA提取使用OMEGA公司的E.Z.N.ATM Plant RNA Kit試劑盒進行提取，提取方法嚴(yán)格按照試劑盒說明書進行。并將獲得的總RNA用瓊脂糖凝膠電泳分析RNA的降解程度及是否存在污染，用液相色譜Agilent 2100和OD260/280的比值檢測RNA的純度。旱柳總RNA的提取、質(zhì)量檢測及無參考基因組序列的轉(zhuǎn)錄組高通量測序分析均由深圳華大基因生物科技有限公司進行。

1.2.2 文庫構(gòu)建及轉(zhuǎn)錄組測序 樣品檢測合格后，加入Fragmentation buffer將mRNA打斷成短片段，用六堿基隨機引物（Random hexamers）以mRNA為模板合成一鏈cDNA，然后加入緩沖液、dNTPs和DNA polymerase I合成二鏈cDNA，雙鏈cDNA再利用AMPure XP beads純化。純化后雙鏈cDNA末端修復(fù)、加A尾并連接測序接頭，片段大小選擇是用AMPure XP bead，最終的cDNA文庫通過PCR富集得到。cDNA 文庫質(zhì)量用Agilent 2100檢測，再通過Illumina Hi-Seq2500平臺完成RNA-Seq測序分析。

1.2.3De nove組裝 獲得的原始數(shù)據(jù)Raw reads經(jīng)過濾得到Clean reads。使用Trinity對clean reads組裝得到contigs，之后再用Tgicl對轉(zhuǎn)錄組樣本進行兩次聚類去冗余得到最終Unigene。Unigene分為兩部分，一部分是clusters，是進一步去冗余后的結(jié)果；其余的是singletons，指沒有聚類、單獨的 Unigene。每個樣品測序產(chǎn)出不少于 6 Gb Clean data，Q30堿基百分比達到85%。

1.2.4 生物信息學(xué)分析 使用BLAST 對Unigene序 列 與 Nr、Nt、KOG、KEGG、GO、Swiss-Prot及InterPro數(shù)據(jù)庫進行比對（Evalue＜1e-5），獲取與旱柳Unigene 具有最高序列相似性的蛋白，得到Unigene 的蛋白功能注釋信息。

本研究從兩個角度分析旱柳在Cd脅迫后的差異基因表達情況:（1）同一時間不同濃度Cd脅迫的差異基因表達情況的對比分析，命名為L1-CK1（即L1和CK1兩個樣品的基因表達對比，下同）、H1-CK1、L7-H7、L30-CK30、H30-CK30；（2）同一濃度不同時間的差異基因表達情況對比分析，分為2.5 mg/L組和50 mg/L組，命名為L1-CK0、H1-CK0、L7-L1、H7-H1、L30-L7、H30-H7。

根據(jù)差異表達分析結(jié)果，使用 Blast2GO和WEGO軟件對差異表達的Unigene進行GO功能分類［13］。將差異表達的Unigene注釋到KEGG數(shù)據(jù)庫中得到注釋信息，之后進行pathway分析［14］。得到2.5 mg/L和50 mg/L濃度Cd脅迫后旱柳的生理生化響應(yīng)結(jié)果。

1.2.5 差異表達基因q-PCR驗證 利用 primer premier 5.0 對Cd脅迫響應(yīng)的相關(guān)基因和內(nèi)參基因肌動蛋白基因（ACT）進行 q-PCR 引物設(shè)計。在BIORAD公司的C1000TM熒光定量PCR儀進行real-time PCR反應(yīng)，每個樣品重復(fù)3次，反應(yīng)體系20 μL。以內(nèi)參基因ACT為對照，利用公式2-△△Ct計算其相對表達量。

2 結(jié)果

2.1 總RNA質(zhì)檢結(jié)果

10個樣品中除CK7外9個樣品的RNA質(zhì)量較好，能滿足建庫測序要求，可用于RNA轉(zhuǎn)錄組測序，且總量滿足1次但不足2次建庫需要。CK7因RNA提取質(zhì)量不佳，所以沒有進行轉(zhuǎn)錄組測序。

2.2 轉(zhuǎn)錄組測序與De nove組裝

9個樣品的raw reads在72.22 Mb到85.59 Mb之間，過濾后總clean reads在 65.17 Mb到67.15 Mb之間；Q20、Q30以及過濾前后reads的比率分別是95.04%到 97.38%，86.40%到 90.91%，72.75%到89.4%（表1），這表明低質(zhì)量的堿基比率低，測序質(zhì)量好，可以用于de novo組裝。用Trinity對其進行組裝，所有樣品的總reads均超過10萬，總長度超過7千萬，平均長度超過600，N50超過1 000，GC含量42%左右。之后使用 Tgicl 對轉(zhuǎn)錄本進行聚類去冗余得到Unigene，聚類后的 Unigene質(zhì)量（表2），10個樣品的Unigene數(shù)為53 705到62 933個，總長度從43 247 382 nt 到50 853 969 nt，N50 值都在1 200之上，這表明組裝質(zhì)量高，可用于后續(xù)分析。

表1 過濾前后的reads質(zhì)量統(tǒng)計表

表2 Unigene的質(zhì)量指標(biāo)

2.3 對Unigene的功能注釋

得到的Unigene有83.93%在7個數(shù)據(jù)庫（Nr、Nt、KOG、KEGG、GO、Swiss-Prot及 InterPro）中任意一個注釋到。Nr數(shù)據(jù)庫進行比對篩選，篩選條件為E值小于10-5，結(jié)果顯示:有74 919條Unigenes在Nr數(shù)據(jù)庫中找到同源蛋白序列，占總Unigenes數(shù)的73.02%。Nt數(shù)據(jù)庫注釋到的Unigenes最多，占總Unigenes的81.12%，但被全部數(shù)據(jù)庫注釋到的Unigenes僅為29.04%，說明有許多旱柳Unigenes的功能還沒有完全明確。

2.4 Unigene GO功能分類

根據(jù) Nr 數(shù)據(jù)庫注釋的信息結(jié)果顯示，有51 494條Unigene 映射到GO不同的功節(jié)點（Term）上，根據(jù)生物過程（Biological process）、細胞組成（Cellular component）、分子功能（Molecular function）進行分類。由于經(jīng)常有同一個轉(zhuǎn)錄本映射到不同節(jié)點現(xiàn)象，所以在生物過程中有106 517條，占42.69%，參與代 謝 過 程（Metabolic process）（22 955，21.54%）、細 胞 過 程（Cellular process）（21 034，19.75%）的Unigene最多；在細胞組成中有92 900條，占37.24%， 參 與 膜（Membrane）（15 182，16.34%）、細胞器（Organelle）（13 010，14.01%）、細胞（Cell）（16 820，18.11%）和細胞部分（Cell part）（16 671，17.95%）的Unigene最多；在分子功能中有50 054條，占20.06%，參與催化活性（Catalytic activity）（20 318，40.59%）和結(jié)合蛋白（Binding）（20 533，41.02%）的Unigene最多（圖1）。

圖1 Unigene的GO分類圖

2.5 Unigene 的KEGG注釋和代謝通路分析

將旱柳差異表達基因在KEGG數(shù)據(jù)庫中進行比對分析顯示，53 966個Unigene分別注釋到KEGG數(shù)據(jù)庫的6個一級層級和21個二級層級中，共包含137條通路。其中在同一個pathway通路中10個以上差異表達基因主要涉及細胞過程、環(huán)境信息處理、遺傳信息處理、人類疾病、代謝以及有機系統(tǒng)6個方面。這6個方面中代謝類的涉及最多，包括碳水化合物、氨基酸、多糖生物合成、能量、次生物質(zhì)等；涉及人類疾病的最少，只有2個，分別是內(nèi)毒素和代謝疾病（Endocine and metabolic diseases）和耐藥性∶抗菌素（Drug resistance∶Antimicrobial）。在旱柳Cd脅迫后KEGG分析中發(fā)現(xiàn)涉及人類疾病的基因變化，這也表明Cd對旱柳中有關(guān)人疾病的內(nèi)毒素和抗菌性的代謝有影響。同時，pathway富集結(jié)果也表明，前述的25個Cd脅迫響應(yīng)候選基因主要參與代謝途徑（Metabolic pathway）（10541，19.53%）、次生代謝產(chǎn)物的合成（Synthesis of secondary metabolites）（5473，10.14%）、植物MAPK信號通路（Plant MAPK signaling pathway）（1592，2.95%）、 植 物 病原互作（Plant pathogen interaction）（1761，3.26%）、植物激素信號轉(zhuǎn)導(dǎo)等通路（Plant hormone signal transduction）（1862，3.45%）（表3）。由此推測旱柳可能主要是通過這些代謝通路調(diào)控Cd抗逆性過程。

表3 富集KEGG通路分析

2.6 差異表達基因分析

FDR（False discovery rate）值越小，表明基因表達差異越大，差異越顯著。本研究在差異基因時設(shè)置的閾值為Fold Change≥2.00，F(xiàn)DR≤0.001，Q-value＜0.05。

由表4可知，相同時間不同濃度的差異表達基因（Differentially expressed genes，DEG）分析顯示:第1天組（L1-CK1和H1-CK1）中的差異基因數(shù)為8 550個（3 660個上調(diào)和4 890個下調(diào)）；第7天組（L7-H7）的差異基因數(shù)為3 444個（1 307個上調(diào)和2 137個下調(diào)）；第30 天組（L30-CK30和H30-CK30）中的差異基因數(shù)為11 428個（5 055個上調(diào)和6 373個下調(diào)）。相同濃度不同時間的DEG分析顯示:2.5 mg/L組（L1-CK0，L7-L1，L30-L7）中的差異基因數(shù)為26 623個（14 531個上調(diào)和12 092個下調(diào)；50 mg/L組（H1-CK0，H7-H1，H30-H7）中的總差異基因數(shù)為32 154個（18 318個上調(diào)和13 836個下調(diào)），這說明50 mg/L濃度的Cd脅迫比2.5 mg/L的對旱柳的影響更加明顯。此外，在L1-CK1、H1-CK1、L30-CK30和H30-CK30這4組中有共同的170個DEG（其中46個共同上調(diào)和29個共同下調(diào)），它們的基因表達量差異達2倍（|Log2ration|≥1）以上。170個DEG中又有25個其表達量差異是大于1 000倍（|Log2ration|≥10）以上的，其中19個是基因表達上調(diào)的，6個是基因表達下調(diào)的（表5）。這25個基因表達變化如此之大，說明它們在旱柳受Cd脅迫耐受過程中等起到了重要作用，可以作為旱柳Cd脅迫響應(yīng)的重要候選基因。

表4 差異表達基因統(tǒng)計

表5 抗逆相關(guān)基因詳細情況

2.7 差異基因的q-PCR

黃酮類化合物生物合成和油菜素內(nèi)酯合成的關(guān)鍵基因:FLS，F(xiàn)3H和ROT3，這3個基因在旱柳Cd脅迫響應(yīng)代謝通路中表達量顯著上升，鈣依賴性蛋白激酶2基因表達量顯著下調(diào)，因此我們選擇這4個基因進行q-PCR驗證轉(zhuǎn)錄組測序分析的準(zhǔn)確性和可靠性，肌動蛋白基因（ACT）因其具有穩(wěn)定性而被選擇為內(nèi)參基因［15］。圖2-6的q-PCR結(jié)果顯示:5個差異表達基因在轉(zhuǎn)錄組測序中的表達變化與熒光定量結(jié)果基本一致，說明轉(zhuǎn)錄組測序分析的結(jié)果真實可靠。

圖2 CL11387.Contig1基因q-PCR結(jié)果

圖3 Unigene3229_All基因q-PCR結(jié)果

圖4 CL13372.Contig3基因q-PCR結(jié)果

圖5 CL1644.Contig4基因q-PCR結(jié)果

3 討論

圖6 CL10874.Contig2_All基因q-PCR結(jié)果

旱柳是中國土生土長的高大喬木，具有易扦插、易成活、生長快、枝葉發(fā)達和抗逆性強等特點，是一種非常適合用修復(fù)土壤重金屬污染的物種。為研究旱柳響應(yīng)不同濃度Cd脅迫的基因表達情況，選擇Cd脅迫濃度分別為2.5 mg/L和50 mg/L的旱柳葉片作為轉(zhuǎn)錄組測序樣品，通過旱柳基因表達的差異來研究旱柳響應(yīng)Cd脅迫的分子機制。轉(zhuǎn)錄組數(shù)據(jù)進行過濾、拼接、組裝后共獲得102 595個基因。

近年來對植物響應(yīng)重金屬脅迫后的基因表達研究已取得一系列成果，鑒定和推測出不少關(guān)于重金屬響應(yīng)的基因，包括ABC轉(zhuǎn)運家族、鋅和錳轉(zhuǎn)運蛋白家族、重金屬ATP酶（HMAs）家族、陽離子擴散促進者（CDF）家族、植物螯合肽合成基因（PCS）、金屬硫蛋白基因（MT）和金屬忍耐蛋白基因（MTP）等［16］。本研究發(fā)現(xiàn)在旱柳響應(yīng)Cd脅迫的過程中，金屬硫蛋白基因、ABC轉(zhuǎn)運家族、鋅和錳轉(zhuǎn)運蛋白家族的基因表達會隨著Cd濃度的增高和脅迫時間的增加而有較大的變化，L7-L1、L30-L7、H7-H1、H30-H7組中金屬硫蛋白共9個基因表達上調(diào)、11個下調(diào)；ABC轉(zhuǎn)運家族共90個上調(diào)、88個下調(diào)，其中H30-H7中32個上調(diào)，50個下調(diào)；鋅和錳轉(zhuǎn)運蛋白家族有28個下調(diào)，8個上調(diào)，上面列的其它重金屬響應(yīng)基因的差異表達在本研究沒有發(fā)現(xiàn)。這些結(jié)果表明旱柳在應(yīng)對Cd脅迫的過程中金屬硫蛋白基因、ABC轉(zhuǎn)運家族及鋅和錳轉(zhuǎn)運蛋白起了很大的作用，也說明Cd脅迫會影響旱柳對其他金屬離子的吸收，鋅和錳轉(zhuǎn)運蛋白表達量的變化可能是因為Cd在生命體內(nèi)會置換鋅，導(dǎo)致生命體需鋅的蛋白質(zhì)發(fā)生“饑餓”，為了生命正?；顒拥倪M行，旱柳加快了鋅的轉(zhuǎn)運過程，從而使鋅、錳蛋白轉(zhuǎn)運蛋白家族的表達發(fā)生變化。ABC轉(zhuǎn)運家族基因表達量發(fā)生明顯變化，說明旱柳在Cd脅迫響應(yīng)中細胞物質(zhì)跨膜轉(zhuǎn)運主要是通過ABC轉(zhuǎn)運家族蛋白來實現(xiàn)的，也說明ABC轉(zhuǎn)運家族在旱柳細胞運送重金屬相關(guān)物質(zhì)中起重要作用。

黃酮類化合物是植物體內(nèi)一類在抗菌、抗逆性上具有重要的作用的化合物，是一種良好的抗氧化劑，能夠清除生物膜周圍和細胞內(nèi)的H2O2［17-18］。有

研究表明，黃酮化合物的代謝通路在水稻響應(yīng)Cu2+和Cd2+的脅迫以及組培旱柳苗響應(yīng)Cd2+中起了積極作用［19-20］。孫愛清等［21］利用 Solexa高通量測序?qū)ㄉ珊淀憫?yīng)基因表達基因分析，發(fā)現(xiàn)有9個類黃酮代謝相關(guān)基因在干旱脅迫下顯著表達；宋中邦［22］通過煙草基因組數(shù)據(jù)庫中搜索獲得兩個類黃酮合成酶基因，即NtFLS1和NtFLS2，二者表達模式一致，均在葉片中高水平表達。本研究通過KEGG代謝通路圖進行分析發(fā)現(xiàn)，旱柳黃酮類化合物的生物合成通路在樣品中均有大量的差異基因上調(diào)，L1-CK1、H1-CK1、L30-CK30、H30-CK30組 中 一 共 170個差異基因，其中有46個共同上調(diào)基因，29個共同下調(diào)基因，而這46個共同上調(diào)基因中包含了FLS、F3H基因。這說明旱柳在應(yīng)對Cd脅迫時，黃酮類化合物起了重要作用。

同時，本研究發(fā)現(xiàn)油菜素內(nèi)酯合成通路的 3，6-脫氧油菜素淄酮酶（ROT3）在除L1-CK1外的其他樣品中都是上調(diào)。肖瑞雪等［23］研究發(fā)現(xiàn)油菜素內(nèi)酯是一種被譽為第六激素的天然激素，其生理活性甚至超過現(xiàn)在5種激素，能提高植物的抗逆性。靳開川等［24］結(jié)合遺傳學(xué)、基因與蛋白質(zhì)組學(xué)、細胞生物學(xué)等多學(xué)科方法和手段，發(fā)現(xiàn)油菜素內(nèi)酯在植物的抗逆性（抗干旱、高鹽、高溫、低溫、重金屬）過程中的起重要作用。王喆等［25］研究表明油菜素內(nèi)酯可能會與膜蛋白結(jié)合，然后通過減輕重金屬脅迫來提高代謝活性，從而降低植物對重金屬的攝取能力［26］。本研究在 H1-CK1、L30-CK30、H30-CK30中油菜素內(nèi)酯合成通路的 3，6-脫氧油菜素淄酮酶（ROT3）上調(diào)是由于Cd2+在旱柳葉片中積累，旱柳為減輕Cd2+帶來的損害，而促進了油菜素內(nèi)酯的合成。

本研究對旱柳受Cd脅迫的轉(zhuǎn)錄組進行了差異基因的分析，找到了受脅迫的相關(guān)基因以及與之相關(guān)的代謝通路，但這些基因是如何互相作用調(diào)節(jié)旱柳Cd脅迫，還需要在分子方面進行更加深入的分析。

4 結(jié)論

旱柳在受重金屬Cd脅迫后展現(xiàn)出一定程度的脅迫抗性，同時基因表達與代謝過程發(fā)生了重大的變化。（1）鎘脅迫后不同濃度不同時間的旱柳轉(zhuǎn)錄組測序共獲得102 595個Unigenes；（2）篩選得到與Cd脅迫響應(yīng)密切相關(guān)的基因25個，其中金屬硫蛋白、ABC轉(zhuǎn)運蛋白、鋅和錳轉(zhuǎn)運蛋白基因在旱柳植物抗逆性過程中發(fā)揮重要作用；（3）Cd脅迫的GO條目主要集中在代謝過程、細胞過程、膜、細胞器、細胞、細胞部分、催化活化和結(jié)合蛋白上；（4）油菜素內(nèi)酯合成通路的3，6-脫氧油菜素淄酮酶（ROT3）和黃酮類化合物合成通路的黃酮醇合成酶（FLS）、黃烷酮-3-羥化酶（F3H）均明顯上調(diào)。