水稻芽期響應鎘脅迫的轉錄組分析

王燕寧，黃 濤，吳光亮，黃詩穎，鐘 奇，王 鵬，楊朦朦，李才敬，程 琴，賀浩華，2，金 偉，郭 凌，邊建民，2

(1.作物生理生態與遺傳育種教育部重點實驗室，江西 南昌 330045；2.江西省水稻高水平工程研究中心，江西 南昌 330045；3.南昌市農業技術推廣中心，江西 南昌 330200)

水稻是我國的主要糧食作物之一，具有富集鎘的特性[1-2]。水稻根、莖、葉、籽粒容易富集鎘，嚴重影響水稻的光合作用、蒸騰作用和其他生理過程，而且會產生過量的氧自由基影響水稻的正常生長發育[3]，鎘還會通過稻米進入人體，危害人類健康[4-5]。近年來，水稻鎘危害引起人們廣泛關注，因此，了解水稻鎘脅迫的遺傳機制，培育耐鎘水稻品種對于降低水稻鎘危害、保障稻米安全具有重要作用。

水稻鎘脅迫相關性狀屬于復雜性狀，受到多個遺傳位點共同調控，近年來研究人員通過外源鎘處理的方法，定位了多個與根長、株高、干質量和鎘含量相關的QTLs，大多數效應較小。Yan等[6]通過一個重組自交系群體(RILs)定位到與鎘含量相關的5個主效QTLs(scc10、gcc3、gcc9、gcc11和sgr5)，其中sgr5是一個新的潛在QTL，對鎘從地上部到籽粒的分布影響較大。Ishikawa等[7]利用Sasanishiki和Habataki構建的回交自交系(BILs)定位群體在7號染色體上鑒定到一個新型主效QTLqGCd7，qGCd7是提高水稻籽粒鎘含量的特異性QTL。駱旭添[8]利用Lemont(美國)和Dular(印度)雜交構建的RILs，以葉綠素含量、根長、株高、葉長為指標共檢測到16個加性QTL，涉及水稻第1，2，3，4，6，7，11號等7條染色體。Pan等[9]對生長在不同程度鎘污染土壤上的338份不同稻種資源的鎘積累能力進行了評估，通過全基因組關聯分析(GWAS)研究發現了35個與鎘積累顯著相關的QTL。Li等[10]利用CH891和02428構建的BILs定位了鎘脅迫下種子發芽率相關的6個加性QTL和16對上位性QTL。

鎘相關基因的克隆也取得了一定的進展，在水稻12條染色體上均有分布。Ueno等[11]利用鎘高積累品種Anjana Dhan和日本晴構建的定位群體，克隆了水稻低鎘積累基因OsHMA3，該基因主要影響鎘從地下向地上部的轉運效率以及地上部的鎘積累，此后發現該基因在品種間存在多個自然變異類型。Luo等[12]利用鎘積累差異顯著的2個品種臺中1號和春江06，構建雜交群體克隆到了一個特異性調控水稻葉片鎘積累的主效QTL—CAL1基因。該基因編碼植物防御素類似蛋白，定向調控鎘在葉片等營養器官的積累。鋅鐵調控轉運蛋白基因OsZIP3能降低水稻根中的鎘含量[13]。Chen等[14]分離出一個水稻耐鎘突變體，并認為OsCADT1是造成突變體耐鎘的基因。LCD是一個水稻低鎘基因，定位于細胞質和細胞核，被認為是低鎘水稻育種中的潛在可利用基因。鎘脅迫下，其功能缺失型突變體植株的籽粒鎘含量比野生型植株減少1/2[15]。OsCd1是Yan等[16]利用127份自然群體材料通過GWAS定位到的QTL中分離到的水稻籽粒鎘積累基因，OsCd1屬于MFS超家族轉運蛋白，能轉運鎘離子，參與根系鎘吸收。

多年來，國內外學者已經對水稻鎘脅迫的遺傳機理進行了許多研究，但主要集中在成熟期鎘的吸收和轉運，而水稻幼芽期對鎘脅迫較為敏感[17-18]，因此，如何降低鎘脅迫對水稻幼芽的危害也是目前抗逆育種急需解決的問題。研究還表明，不同水稻品種對鎘的吸收和積累能力存在明顯差異[19]，秈稻和粳稻在鎘脅迫下的抗逆性也有顯著差異，但是造成這種差異的機制還不清楚[20]。

本研究以秈稻品種CH891和粳稻品種02428為供試材料，利用RNA-Seq技術對2個品種在正常條件和鎘處理3 d后幼苗的轉錄本進行測序，并對轉錄組測序結果進行GO富集分析和KEGG通路富集分析，旨在為闡明秈、粳稻響應鎘脅迫的遺傳差異提供理論依據。

1 材料和方法

1.1 試驗材料

供試材料為鎘敏感的秈稻恢復系CH891與耐鎘的粳稻品種02428。

1.2 水稻種子鎘處理與樣品采集

挑選健康飽滿的谷粒用1%的次氯酸鈉溶液浸泡消毒15 min，用蒸餾水沖洗3～4次。將消毒過的谷粒浸種24 h，25 ℃催芽1～2 d。分別選取CH891和02428大小、露白基本一致的30粒種子放置于鋪有單層濾紙的培養皿中，處理組用濃度為200 mg/L的CdCl2溶液處理，對照組用相同體積的蒸餾水處理，處理組和對照組均設置3次生物學重復。處理組和對照組放入同一個光照培養箱內進行培養。白天溫度設置為(27.0±0.5)℃，夜間溫度設置為(25 .0±0.5)℃，光期時長設為14 h，暗期時長設為10 h。第3天時將2個品種水稻鎘處理和正常條件下的幼苗分別取下，用錫箔紙包裹置于液氮中速凍，放置在-80 ℃冰箱保存備用。

1.3 總RNA提取

利用TRIzol試劑提取第3天鎘處理和對照組CH891和02428的RNA，利用1%的瓊脂糖凝膠電泳檢測RNA質量，利用2100 Bioanalyzer檢測RNA純度和濃度，RNA質量合格后送去聯川生物公司。

1.4 cDNA文庫構建和轉錄組測序

使用mRNA片段作為模板構建cDNA文庫，cDNA文庫構建由聯川生物公司完成，最后在Illumina HiSeq 4000平臺上進行雙末端測序[21]。

1.5 生物信息學分析

將原始讀長數據(Raw read)進行過濾，將過濾掉低質量reads的高質量序列(Clean read)比對到參考基因組，對所有比對上的基因進行表達水平分析，并用FPKM值表示各基因表達水平，再進一步進行基因表達水平差異研究和差異基因功能分析[22-23]。

1.6 差異基因表達分析

利用edgeR對StringTie組裝和定量完畢的基因進行差異分析(顯著差異的閾值為|Log2Fold Change|≥1，P<0.05)。Fold Change是基因的變化倍數值，一般認為數值在2以上的表示該基因的表達量在2個比較組之間存在顯著差異。P-value是各個基因差異表達概率，其數值一般控制在0.05水平內，即<0.05表示在95%的置信區間內該基因是差異表達基因(可過濾假陽性差異表達基因)。此外，通過繪制火山圖可以了解差異表達基因的整體分布情況，根據樣品基因表達譜的相近程度，將差異基因進行聚類分析，直觀地展示基因在不同樣品(或是不同處理)中的表達情況，由此獲取生物學相關信息[24]。

1.7 差異表達基因功能注釋和分類富集分析

所有差異表達基因均在Gene Ontology(GO)和KEGG(Kyoto Encyclopedia of Genes and Genomes)中分別進行GO和KEGG注釋，網址分別為http：//www.geneontology.org/page/download-go-annotatios和www.genome.jp/kegg。利用WEGO和OMICSHARE進一步對注釋得到的差異表達基因進行GO和KEGG分類富集分析，網址分別為http：//wego.genomics.org.cn和http：//www.omicshare.com/tools/index.php/。

1.8 qRT-PCR驗證

為了驗證RNA-Seq結果，通過定量實時RT-PCR(qRT-PCR)證實了幾個基因的不同表達模式。按照TaKaRa反轉錄試劑盒說明書對本試驗樣品RNA進行反轉錄。按照cDNA 1 μL，Premix Ex TaqTMⅡ 10 μL，ddH2O 7 μL，Rox Reference Dye 0.4 μL，上下引物10 μmol/L各0.8 μL的比例配置qRT-PCR反應體系(表1)。反應程序設定為：95 ℃ 1 min；95 ℃ 15 s，55 ℃ 30 s，72 ℃ 60 s，65 ℃ 5 s，95 ℃ 5 s，40個循環，每個基因設置3次生物學重復和3次技術性重復。利用qRT-PCR所得Ct值計算樣品中基因的相對表達量(2-ΔΔCt)[25-26]。

表1 qRT-PCR引物序列Tab.1 The sequences for qRT-PCR primers

2 結果與分析

2.1 鎘脅迫下CH891和02428的表型差異

正常情況下秈稻CH891幼芽顯著長于粳稻02428，鎘處理后，CH891和02428的生長均受到明顯抑制，其中CH891受到抑制的程度大于02428，芽長顯著變短，但鎘處理3 d后CH891和02428的芽長沒有顯著差異，表明CH891幼芽生長對鎘的敏感性高于02428(圖1)。

2.2 樣品RNA質量檢測與測序序列統計和質量評估

使用高通量 Illumina HiSeq 4000(CA，美國)對CH891和02428在鎘處理3 d和正常生長的幼芽進行RNA-Seq，總共產生545 010 200個原始測序讀長數(Raw read)，平均每個樣本產生45 417 517個原始測序讀長數。過濾掉不合格序列后共得到539 524 490個有效讀長(Valid read)，占原始讀長的99.01%以上。所有樣品測序讀長與參考基因組比對率為94.81%～96.82%，其中唯一比對率為62.11%～74.21%，所有樣品的GC含量均在49%以上，測序質量較好，可進行下一步試驗分析(表2)。

表2 測序數據統計和參考基因組上的比對率Tab.2 Sequencing data statistics and mapped ratios to the reference genome

2.3 差異表達基因數量分析

各比較組中差異表達基因數目如圖2所示。對照條件下，CH891與02428(CK3_CH891 vs.CK3_02428)之間存在4 911個差異表達基因，其中上調表達基因1 647個，下調表達基因3 264個；鎘處理條件下，CH891與02428(Cd3_ CH891 vs.Cd3_02428)之間存在2 542個差異表達基因，其中有1 499個上調表達基因，1 043個下調表達基因；CH891鎘處理后與對照相比(Cd3_CH891 vs.CK3_CH891)共存在1 952個差異表達基因，其中上調表達基因1 203個，下調表達749個基因；02428鎘處理后與對照相比(Cd3_02428 vs.CK3_02428)共有1 792個差異表達基因，其中有516個上調表達基因，1 276個下調表達基因。

A、B、C和D.火山分布圖；E.上下調差異表達基因個數。A，B，C and D.Volcano plot；E.The number of up and down-regulated DEGs.

2.4 qRT-PCR驗證

為驗證RNA-Seq數據的準確性，以OsActin作為內參基因，隨機挑選了12個差異表達基因(LOC_Os03g08470、LOC_Os04g39880、LOC_Os09g36700、LOC_Os10g40720、LOC_Os02g02400、LOC_Os04g3830、LOC_Os06g21590、LOC_Os09g31430、LOC_Os01g18744、LOC_Os06g13280、LOC_Os03g03600、LOC_Os08g06110)，使用特異性引物進行qRT-PCR分析，將qRT-PCR結果與差異表達基因的Log2(Fold Change)值進行比較，結果表明，qRT-PCR結果與RNA-Seq趨勢一致(圖3)，表明轉錄組測序結果準確可靠。

圖3 差異表達基因轉錄組與qRT-PCR比較Fig.3 Comparison of differentially expressed gene transcriptome and qRT-PCR results

2.5 差異表達基因分類

通過對同一水稻品種在不同環境和不同水稻品種在相同環境的差異表達基因進行比較，共獲得4個比較組，分別為Cd3_CH891 vs.Cd3_02428、CK3_CH891 vs.CK3_02428、Cd3_CH891 vs.CK3_CH891和Cd3_02428 vs.CK3_02428。在4個比較組中共存在7 204個差異表達基因，這些差異表達基因可以分為15個亞組(圖4)。去掉CK3_CH891 vs.CK3_02428與鎘無關的基因，剩余的7個亞組可以分為3個類別：CH891特異表達的基因(SCR3)、02428特異表達的基因(RCR3)、2個品種均有表達但表達量存在差異的基因(CCR3)。其中，SCR3有849個，RCR3有676個，CCR3有770個(表3)。

圖4 4個比較組差異表達基因的韋恩圖Fig.4 Venn diagrams for DEGs in the four comparison groups

以|log2(Fold change)| ≥2并且至少在一個比較組的一組數據中FPKM≥2為篩選條件，分別對3個類別的基因進行篩選，共獲得554個主效的鎘相關基因，其中194個來自SCR3，154個來自RCR3，206個來自CCR3。

表3 差異表達基因的分類Tab.3 Classification of three categories of DEGs

2.6 基因功能注釋

2.6.1 SCR3基因功能注釋 GO富集分析顯示，轉錄調節和DNA模板、細胞核、膜的整體結構、質膜相關的基因、RNA結合、轉錄活性和序列特異性DNA結合、金屬離子結合相關的基因在SCR3中均有富集(圖5)。KEGG富集分析中，異黃酮生物合成、磷酸肌醇代謝、黃酮類生物合成和花青素生物合成等途徑在SCR3中富集，其中異黃酮生物合成和花青素生物合成途徑顯著富集(圖6)。

1.生物進程；2.轉錄調節、DNA模板；3.轉錄、DNA模板；4.蛋白磷酸化；5.代謝過程；6.傷害反應；7.氧化還原反應；8.脫落酸反應；9.油菜素類固醇反應；10.脂質轉運；11.缺氧反應；12.蛋白質泛素化；13.蛋白酶體介導的泛素化依賴性蛋白質分解代謝過程；14.冷反應；15.多細胞生物發育；16.種子休眠中組織胚胎發育；17.葉發育；18.光刺激反應；19.蛋白質折疊；20.類黃酮葡萄糖醛酸化；21.胞外分泌；22.次級代謝生物合成；23.防御反應；24.赤霉素反應；25.甲基化作用；26.核；27.膜整體結構；28.質膜；29.細胞質；30.細胞外區域；31.胞液；32.葉綠體；33.膜；34.高爾基體；35.線粒體；36.細胞壁；37.植物細胞壁；38.內質網；39.質膜固定組成成分；40.胞間連絲；41.分子功能；42.RNA結合；43.轉錄活性、序列特異性DNA結合；44.金屬離子結合；45.蛋白結合；46.ATP結合；47.羧基酯水解酶；48.蛋白絲氨酸/蘇氨酸激酶活性；49.鋅離子結合；50.水解酶活性。

1.α亞麻酸代謝；2.萜類骨架生物合成；3.類固醇生物合成；4.淀粉和蔗糖代謝；5.倍半萜和三萜類生物合成；6.自噬調節；7.RNA聚合酶；8.RNA降解；9.嘧啶代謝；10.卟啉和葉綠素代謝；11.植物激素信號轉導；12.光合作用；13.磷脂酰肌醇信號系統；14.類苯丙酸生物合成；15.苯丙氨酸、酪氨酸和色氨酸生物合成；16.磷酸戊糖途徑；17.N-聚糖生物合成；18.異黃酮生物合成；19.磷酸肌醇代謝；20.糖基磷脂酰肌醇生物合成；21.硫代葡萄糖苷生物合成；22.類黃酮生物合成；23.醚類脂代謝；24.二帖生物合成；25.角質、木栓質和蠟質的生物合成；26.類胡蘿卜素；27.光合生物的碳固定；28.生物合成；29.堿基切除；30.花青素生物合成。

2.6.2 RCR3基因功能注釋 GO富集分析顯示，轉錄調控和DNA模板、細胞核、質膜、膜整體結構、和胞外區域相關的基因、轉錄因子活性、序列特異性DNA結合、蛋白質結合相關的基因在RCR3都有富集(圖7)。KEGG富集分析顯示，核糖體、自噬調節和胰島素抗性等途徑在RCR3中富集，其中核糖體和自噬調節富集程度最為顯著(圖8)。

2.6.3 CCR3基因功能注釋 GO富集分析顯示，轉錄調節和DNA模板、細胞核、細胞質、膜整體結構、質膜相關的基因、轉錄因子活性、序列特異性的DNA結合、蛋白質結合、DNA結合、鋅離子結合相關的基因在CCR3中富集(圖9)。KEGG富集分析顯示，α-亞麻酸代謝、吞噬體、過氧化物酶體、檸檬烯和蒎烯的降解以及脂肪酸降解途徑在CCR3中富集，其中α-亞麻酸代謝和脂肪酸降解途徑極顯著富集(圖10)。

1.生物過程；2.轉錄調控、DNA模板；3.轉錄、DNA模板；4.氧化還原反應；5.防御反應；6.蛋白磷酸化；7.細胞壁組織；8.脫落酸反應；9.傷害反應；10.鹽脅迫反應；11.蛋白質水解；12.多細胞生物發育；13.一維細胞生長；14.氧化應激反應；15.冷反應；16.獨腳金內酯生物合成；17.茉莉酸反應；18.蛋白酶體介導的泛素依賴性蛋白質分解代謝過程；19.信號轉導；20.脂質轉運；21.類黃酮生物合成過程；22.卡里金反應；23.涉及多維細胞生長的植物型細胞壁修飾；24.蛋白轉運體；25.胞吐作用中的囊泡對接；26.核；27.質膜；28.膜整體結構；29.胞外區域；30.細胞質；31.線粒體；32.膜；33.葉綠體；34.細胞溶質；35.胞間連絲；36.細胞壁；37.細胞組分；38.液泡；39.植物細胞壁；40.細胞內有膜的細胞器；41.分子功能；42.轉錄活性、序列特異性DNA結合；43.蛋白結合；44.DNA結合；45.蛋白絲氨酸/蘇氨酸激酶活性；46.ATP結合；47.金屬離子結合；48.激酶活性；49.核糖體的結構成分；50.鋅離子結合。

1.mRNA監控途徑；2.β-丙氨酸代謝；3.玉米素生物合成；4.纈氨酸、亮氨酸和異亮氨酸降解；5.囊泡運輸中的SNARE相互作用；6.核糖體；7.自噬調節；8.丙酸代謝；9.植物-病原體相互作用；10.植物激素信號轉導；11.苯丙酸類生物合成；12.苯丙氨酸、酪氨酸和色氨酸的生物合成；13.磷酸戊糖；14.氧化磷酸化；15.其他類型O-聚糖生物合成；16.其他聚糖降解；17.N聚糖生物合成；18.胰島素抗性；19.糖苷磷酸肌醇生物合成；20.醣酵解/糖異生；21.果糖和甘露糖；22.二帖生物合成；23.角質、木栓質和蠟質的生物合成；24.TCA循環；25.類胡蘿卜素生物合成；26.碳代謝；27.堿基切除修復；28.抗壞血酸鹽和醛酸鹽；29.精氨酸和脯氨酸代謝；30.花青素生物合成。

1.轉錄調節、DNA模板；2.生物進程；3.轉錄、DNA模板；4.傷害反應；5.蛋白質泛素化；6.氧化還原反應；7.冷反應；8.脫落酸反應；9.水楊酸反應；10.寡肽轉運；11.防御反應；12.氧化應激反應；13.脂肪酸生物合成；14.脫水反應；15.參與細胞內蛋白質代謝過程的蛋白質水解；16.光刺激反應；17.一維細胞生長；18.乙炔激活信號轉導；19.赤霉酸介導的信號轉導；20.次級代謝生物合成；21.鹽脅迫反應；22.幾丁質反應；23.調控防御反應；24.其他有機物反應；25.木質素生物合成；26.核；27.細胞質；28.膜整體結構；29.質膜；30.膜；31.葉綠體；32.細胞內；33.胞外區域；34.線粒體；35.質膜固定成分；36.細胞溶質；37.高爾基體；38.質體；39.葉綠體基質；40.細胞壁；41.轉錄因子活性、序列特異性DNA結合；42.分子功能；43.蛋白質結合；44.DNA結合；45.鋅離子結合；46.序列特異性DNA結合；47.轉運蛋白活性；48.轉移酶活性，轉移糖酰基；49.蛋白泛素化轉移酶活性；50.轉錄調控區DNA結合。

1.α亞麻酸代謝；2.纈氨酸、亮氨酸和異亮氨酸降解；3.銻烯類、二芳基庚烷類和姜酚生物合成；4.植物病原體互作；5.吞噬體；6.過氧化物酶體；7.煙酸鹽和煙酰胺代謝；8.N聚糖生物合成；9.賴氨酸降解；10.亞油酸代謝；11.檸檬烯和蒎烯降解；12.異黃酮類生物合成；13.糖磷脂生物合成；14.糖磷脂生物合成-神經節系列；15.糖酵解/糖異生；16.甘氨酸、絲氨酸和蘇氨酸代謝；17.甘油磷脂代謝；18.半乳糖代謝；19.脂肪酸延伸；20.脂肪酸降解；21.脂肪酸生物合成；22.二帖生物合成；23.角質、木栓質和蠟質的生物合成；24.TCA循環；25.植物晝夜節律；26.碳代謝；27.油菜素類固醇生物合成；28.精氨酸生物合成；29.氨酰基tRNA生物合成；30.丙氨酸、天冬氨酸和谷氨酸代謝。

3 結論與討論

鎘作為一種重金屬污染物，嚴重影響水稻的正常生長發育。水稻鎘抗性機制十分復雜，因品種、遺傳背景不同而存在差異。本研究通過對耐鎘粳稻品種02428和鎘敏感秈稻品種CH891進行鎘處理，并對第3天的幼苗進行差異表達基因分析。結果表明，CH891和02428的代謝途徑在鎘處理3 a后存在差異，表明代謝途徑的差異可能是秈稻和粳稻對鎘脅迫反應不同的原因。因此，深入分析這些代謝途徑中差異表達的基因對于解釋秈稻和粳稻鎘脅迫差異是必要的。

3.1 SCR3中的差異表達基因

差異表達基因的代謝途徑富集分析表明，異黃酮生物合成和花青素生物合成途徑是SCR3的主要富集通路。有研究表明，異黃酮在植物防御反應中起到重要作用[27]。在鎘脅迫下，植物中的細胞壁受到破壞，使其易受到病原體的侵害，從而誘導了異黃酮的表達以抵抗病害[28]。本研究結果可能為水稻在異黃酮水平上對鎘脅迫的響應提供新的見解。花青素是一種抗氧化劑，能清除植物細胞內的氧自由基，減輕活性氧的毒性[29]。Roychoudhury等[30]發現，花青素生物合成能被CdCl2誘導，與本研究結果相一致，表明花青素作為一種抗氧化劑，與水稻鎘脅迫反應相關，可能減輕鎘對水稻的毒害作用。

3.2 RCR3中的差異表達基因

就RCR3而言，不同水稻品種的KEGG代謝通路在鎘脅迫后也存在不同。Bai等[31]研究發現，與核糖體相關的DEGs在冷應激信號轉導中起著關鍵作用。核糖體途徑在本研究中富集于RCR3中，表明核糖體途徑很可能在水稻鎘脅迫下被誘導并發揮作用。自噬可以將應激引起的氧化蛋白和損傷的胞內物質運輸到空泡中降解，從而減少毒性物質在細胞質中的積累[32]，它在植物的發育和對脅迫的響應中起著重要的作用[33]。本研究發現，RCR3中誘導了自噬途徑的調節，進一步驗證了自噬調節在水稻鎘脅迫中發揮了一定的作用。這2條通路是本研究在水稻鎘脅迫下發現的未被大量報道的途徑。

3.3 CCR3中的差異表達基因

代謝途徑差異同樣存在于CCR3中。α-亞麻酸代謝和脂肪酸降解可概括為脂肪酸代謝，它與CCR3中的鎘反應有關。已有研究表明，重金屬能特異性地改變植物體內脂肪酸代謝，鎘能提高脂質過氧化物水平，導致水稻體內脂肪酸含量發生強烈變化[34-35]。本研究中，α-亞麻酸代謝和脂肪酸降解途徑在CCR3中顯著富集，進一步證明脂肪酸代謝可能是水稻對鎘脅迫的應激反應之一。

綜上所述，本研究表明，水稻鎘響應機制存在品種的特異性。SCR3中富集的信號通路大多與次級代謝過程相關，而RCR3中則多為蛋白質代謝。這表明CH891可能通過誘導次生代謝過程對鎘產生應答，從而增強其免疫力和逆境抗性。而02428中蛋白質的產生和降解過程受到影響，從而保護自身不受鎘的侵害。