泡桐根系應(yīng)答鉛鋅礦脅迫的轉(zhuǎn)錄組分析

韓良澤，陳明利，陳永華，陳 昊，賈惠雁，歐陽(yáng)林男

（1. 中南林業(yè)科技大學(xué) a. 環(huán)境科學(xué)與工程學(xué)院；b. 林學(xué)院，湖南 長(zhǎng)沙 410004；2. 國(guó)家林業(yè)和草原局桉樹(shù)研究開(kāi)發(fā)中心，廣東 湛江 524000）

植物對(duì)重金屬的耐性往往是多基因控制的復(fù)雜性狀，響應(yīng)重金屬脅迫的基因通常涉及到金屬離子轉(zhuǎn)運(yùn)、螯合作用、離子區(qū)隔化、抗氧化系統(tǒng)等諸多方面[1]。植物應(yīng)答重金屬脅迫研究有助于揭示植物對(duì)重金屬的應(yīng)答機(jī)理，從而可以采取相應(yīng)的措施去降低重金屬對(duì)植物的毒害。轉(zhuǎn)錄組學(xué)作為基因組與蛋白質(zhì)功能組的重要橋梁，其研究已經(jīng)成為揭示生物生長(zhǎng)發(fā)育規(guī)律以及尋找逆境脅迫機(jī)制等方面的最佳手段。轉(zhuǎn)錄組所測(cè)定的對(duì)象是某些細(xì)胞在特定狀態(tài)下的mRNA，而泡桐尚未有公布的參考基因組，因此，對(duì)于無(wú)參基因的轉(zhuǎn)錄組，可以通過(guò)高通量測(cè)序技術(shù)對(duì)某種物種某一器官在特定狀態(tài)下進(jìn)行轉(zhuǎn)錄組測(cè)序，組裝得到基因，再進(jìn)行功能注釋。這種方法相比較常規(guī)的分子手段，具有高效、價(jià)格低、覆蓋面廣等優(yōu)點(diǎn)[2-3]。因此可以用此技術(shù)去發(fā)掘特定狀態(tài)下的關(guān)鍵基因。

白花泡桐隸屬于玄參科泡桐屬，是一種喜光、耐寒性不強(qiáng)、對(duì)瘠薄土壤有較強(qiáng)適應(yīng)性的速生樹(shù)種。課題組前期在野外調(diào)查中發(fā)現(xiàn)泡桐在鉛鋅礦中可以存活[4]，說(shuō)明泡桐有較強(qiáng)的適應(yīng)能力。在礦上生長(zhǎng)的泡桐，其體內(nèi)的重金屬含量明顯高于背景區(qū)[5]，可當(dāng)作金屬富集植物，用于尾礦的生態(tài)修復(fù)。泡桐根系作為與污染源直接接觸的部位，其對(duì)重金屬的吸收、積累及耐性機(jī)制等方面已有研究[6]，但是，在重金屬脅迫下，泡桐根轉(zhuǎn)錄組調(diào)節(jié)機(jī)制及相關(guān)基因表達(dá)尚未有明確研究，因此，以泡桐為植物材料，以鉛鋅礦、紅壤為基質(zhì)，進(jìn)行盆栽實(shí)驗(yàn)，并對(duì)所得的泡桐根進(jìn)行轉(zhuǎn)錄組測(cè)序，分析泡桐根部響應(yīng)重金屬脅迫的相關(guān)差異表達(dá)基因，以期進(jìn)一步研究受重金屬脅迫差異表達(dá)的基因與自然礦床或開(kāi)采礦山污染土壤的關(guān)系，并為泡桐應(yīng)用于礦山生態(tài)修復(fù)提供基礎(chǔ)理論依據(jù)。

1 材料與方法

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

供試材料白花泡桐購(gòu)置于河南某苗圃，選取粗細(xì)一致的泡桐根置于紅壤中進(jìn)行土培生芽，待幼苗長(zhǎng)至株高10 cm 左右時(shí)選取高度一致的幼苗移栽。鉛鋅礦取自于郴州市資興鉛鋅礦尾庫(kù)，紅壤取自于中南林業(yè)科技大學(xué)生態(tài)園0 ～20 cm 表層土壤，鈣鎂磷肥購(gòu)置于湖南長(zhǎng)沙花卉市場(chǎng)。

1.2 實(shí)驗(yàn)設(shè)計(jì)

實(shí)驗(yàn)在中南林業(yè)科技大學(xué)苗圃內(nèi)進(jìn)行，將鉛鋅礦、紅壤按照質(zhì)量比分別進(jìn)行攪拌成2 個(gè)實(shí)驗(yàn)組R（紅壤土+0.3 kg 鈣鎂磷肥）、T（鉛鋅礦渣+0.3 kg 鈣鎂磷肥），每盆儲(chǔ)量10 kg，每盆1 株，每個(gè)處理3 次重復(fù)，定期澆水，除蟲(chóng)，實(shí)驗(yàn)周期為1 a。土壤重金屬含量見(jiàn)表1。

表1 盆栽土壤重金屬質(zhì)量分?jǐn)?shù)?Table 1 Heavy metal content in potted soil

1.3 實(shí)驗(yàn)方法

1.3.1 RNA 提取與鑒定

將種植一年后的泡桐根部洗凈，經(jīng)液氮速凍后，用于總RNA 的提取。使用Qubit?2.0 Flrometer（加利福尼亞州Life Technologies）中的Qubit?RNA 分析試劑盒檢測(cè)RNA 濃度。使用納米光度計(jì)?分光光度計(jì)（Implen，CA，USA）檢測(cè)RNA 純度，使用Agilent Bioanalyzer 2100 系統(tǒng)（美國(guó)加利福尼亞州Agilent Technologies）的RNA nano 6000 分析試劑盒評(píng)估RNA 完整性。

1.3.2 文庫(kù)構(gòu)建及轉(zhuǎn)錄組測(cè)序

每個(gè)樣品取總RNA1.5 μg，用帶有Oligo(dT)的磁珠富集mRNA。隨后加入Fragmentation buffer將mRNA 打斷成短片段，以此mRNA 為模板，用六堿基隨機(jī)引物（Random hexamers）合成cDNA 第一鏈，隨后加入緩沖液、dNTPs、DNA polymerase I 和RNase H 合成cDNA 第二鏈，再用AMPure XP beads 純化雙鏈cDNA。將純化的雙鏈cDNA 進(jìn)行末端修復(fù)、加A 尾并連接測(cè)序接頭后，再用AMPure XP beads 選擇用于測(cè)序的片段，并進(jìn)行PCR 擴(kuò)增，用AMPure XP beads 純化PCR 產(chǎn)物，即獲得最終的測(cè)序文庫(kù)。測(cè)序由北京諾禾致源生物信息科技有限公司完成。

1.3.3 測(cè)序數(shù)據(jù)質(zhì)量評(píng)估及轉(zhuǎn)錄本拼接

將基于測(cè)序平臺(tái)所得到的原始數(shù)據(jù)進(jìn)行內(nèi)部Perl 腳本處理，其中包括去除有接頭的reads、去除低質(zhì)量的reads、去除含有超過(guò)10%的不確定堿基的reads。同時(shí)，通過(guò)Q20、Q30、GC 含量對(duì)clean reads 進(jìn)行評(píng)估。使用Trinity[7]軟件對(duì)質(zhì)量合格的reads 進(jìn)行組裝，默認(rèn)情況min_kmer_cov 設(shè)置為2，所有其他參數(shù)設(shè)置為默認(rèn)值。

1.3.4 基因功能注釋

將組裝獲得的基因序列與NCBI 官方非冗余蛋白質(zhì)序列數(shù)據(jù)庫(kù)（NCBI non-redundant protein sequence database，Nr）、NCBI 官方非冗余核苷酸序列數(shù)據(jù)庫(kù)（NCBI non-redundant nucleotide sequence database，Nt）、KEGG 直系同源數(shù)據(jù)庫(kù)（KEGG Orthology database，KO）、Swiss-Prot 數(shù) 據(jù) 庫(kù)（Swiss-Prot）、蛋白質(zhì)家族數(shù)據(jù)庫(kù)（Protein family database，Pfam)、基因功能分類(lèi)體系數(shù)據(jù)庫(kù)（Gene Ontology database，GO）和真核生物直系同源序列 數(shù) 據(jù) 庫(kù)（euKaryoticOrthology Groups database，KOG)進(jìn)行 Blast 比對(duì)（其中 Nr、Nt、Swiss-Prot 數(shù)據(jù)庫(kù)的比對(duì)e值≤1e-5，KOG 比對(duì)e值≤1e-3，Go比對(duì)e值 =1e-6，KEGG 比對(duì)e值 =1e-10）。

1.3.5 SSR 監(jiān)測(cè)及引物設(shè)計(jì)

使 用 MISA(http://pgrc.ipk-gatersleben.de/misa/misa.html) 鑒定轉(zhuǎn)錄組中SSR 序列，并用Primer3(https://sourceforge.net/projects/primer3/獲取軟件信息)設(shè)計(jì)擴(kuò)增SSR 序列的引物。

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

使用DESeq對(duì)2個(gè)樣品組進(jìn)行差異表達(dá)分析，使用q值調(diào)整P值[8]。設(shè)置了 qvalue ＜ 0.005 ＆|log2（foldchange）|＞1 作為顯著差異表達(dá)的閾值。分別用GOseq[9]和KOBAS[10]對(duì)獲得的差異表達(dá)基因進(jìn)行通路富集性分析。

1.3.7 差異表達(dá)基因的qRT-PCR 驗(yàn)證

挑選轉(zhuǎn)錄組測(cè)序得到的差異基因，以Actin 基因作為內(nèi)標(biāo)基因，引物對(duì)為：Actin-F（5′-TTGTGCTTAGTGGTGGGT-3′） 和 Actin-R（5′-CTGGAAGGTGCTGAGGGAT-3′）。 通 過(guò) 使用Primer6.0 軟件設(shè)計(jì)引物并進(jìn)行熒光定量分析。以 2 μgRNA-seq 測(cè)序的根總 RNA 為模板，合成cDNA 第一鏈。qRT-PCR 按照以下步驟進(jìn)行擴(kuò)增：95 ℃預(yù)變性 10 min，95 ℃變性 5 s，72 ℃退火 40 s，40 個(gè)循環(huán)。按照 2?ΔΔCT[11]計(jì)算基因相對(duì)表達(dá)量。

2 結(jié)果與分析

2.1 測(cè)序數(shù)據(jù)質(zhì)量評(píng)估

提取泡桐根部RNA 后由Nova-seq 高通量測(cè)序，共獲得366 701 960 條reads 序列，總堿基數(shù)為53.90 Gb，平均GC 含量達(dá)44.58%。序列質(zhì)量評(píng)估值Q20、Q30 均達(dá)90%以上（表2），說(shuō)明通過(guò)測(cè)序所得到原始序列數(shù)據(jù)較好，可用于后續(xù)的序列組裝。用Trinity 軟件對(duì)得到的clean reads進(jìn)行組裝，共獲得311 822 條轉(zhuǎn)錄本以及126 061條基因序列，最短轉(zhuǎn)錄本長(zhǎng)度為301 bp，最長(zhǎng)轉(zhuǎn)錄本長(zhǎng)度為96 533 bp，平均單條轉(zhuǎn)錄本長(zhǎng)度為1 420 bp，N50 為2 195 bp，最短基因和最長(zhǎng)基因均和轉(zhuǎn)錄本長(zhǎng)度一樣，平均單條基因?yàn)? 068 bp，基因序列在長(zhǎng)度為200 ～500 bp 區(qū)間數(shù)量最多，約為總數(shù)的36%（表3）。

表2 樣品測(cè)序數(shù)據(jù)統(tǒng)計(jì)?Table 2 Summary of sample sequencing data

表3 泡桐根部基因和轉(zhuǎn)錄本長(zhǎng)度分布Table 3 Distribution of all-genes and transcripts length of Paulownia fortunei root

2.2 基因功能注釋與分類(lèi)

2.2.1 功能注釋序列比較

將獲得的基因序列在常用數(shù)據(jù)庫(kù)（Nr，Nt，Pfam，KOG/COG，Swiss-Prot，KEGG，GO） 中進(jìn)行序列對(duì)比，最終獲得96 133 條有注釋信息的基因序列，注釋率為76.25%（表4）。經(jīng)Nr 數(shù)據(jù)庫(kù)比對(duì)后發(fā)現(xiàn)，泡桐根部基因序列與芝麻Sesamum indicum基因序列同源性最高，達(dá)到20.9%，其后依次為栓皮櫟Quercus variabilis（12.5%），紫 花 風(fēng) 鈴Handroanthus impetiginosus（10.6%）和蓖麻Ricinus communis（5.3%），與其他物種同源的有32 416 條（50.7%）。注釋后相似度達(dá)60%～80%、80%～95%、95%～100%的基因分別有23 337 條（36.5%）、18 222 條（28.5%）、4 795 條（7.5%）。

表4 泡桐根部Unigene 序列在各數(shù)據(jù)庫(kù)的功能注釋情況Table 4 Functional annotations of Unigene sequences in databases

2.2.2 GO 數(shù)據(jù)庫(kù)注釋與分類(lèi)

GO 作為轉(zhuǎn)錄組數(shù)據(jù)研究中一套國(guó)際標(biāo)準(zhǔn)化的基因功能描述的分類(lèi)系統(tǒng)，可用來(lái)描述基因及其蛋白產(chǎn)物的特性。將126 061 條基因與GO 數(shù)據(jù)庫(kù)對(duì)比，共有68 840 條基因序列獲得了注釋，涉及3 個(gè)功能大類(lèi)、55 個(gè)功能小類(lèi)，從生物過(guò)程（Biological process）功能分類(lèi)來(lái)看，主要集中在細(xì)胞過(guò)程（Cellular process）、代謝過(guò)程（Metabolic process）、 單 一 有 機(jī) 體 過(guò) 程（Single-organism process）。從細(xì)胞組分（Cellular component）功能分類(lèi)來(lái)看，細(xì)胞（cell）和細(xì)胞部分（Cell part）所占比例較高。從分子功能（Molecular function）分類(lèi)來(lái)看，泡桐根主要涉及結(jié)合（Binding）、催化活性（Catalytic activity）、轉(zhuǎn)運(yùn)活性（Transporter activity）等，其中，僅結(jié)合和催化兩大類(lèi)約占總數(shù)的80%以上（圖1）。

泡桐作為一種耐性植物，其根部涉及到的響應(yīng)重金屬脅迫的基因眾多，如金屬離子響應(yīng)（Stress response to metal ion）543 條、金屬酶活性調(diào)節(jié)（Metalloenzyme regulator activity）263 條、金屬群結(jié)合（Metal cluster binding）248 條、氧含量響應(yīng)（Responseto oxygen levels）854 條、谷胱甘肽酶轉(zhuǎn)移活性（Glutathione transferase activity）762條、細(xì)胞壁增厚防御響應(yīng)（Defense response by cell wall thickening）942 條等。多樣的通路為泡桐的抗重金屬育種提供理論基礎(chǔ)。

2.2.3 KOG 數(shù)據(jù)庫(kù)注釋與分類(lèi)

KOG 作為NCBI 針對(duì)真核生物的蛋白同源數(shù)據(jù)庫(kù)，目前有4 852 個(gè)分類(lèi)。將某個(gè)蛋白質(zhì)與所有KOGs 條目中的蛋白質(zhì)進(jìn)行比對(duì)，可將其歸入適當(dāng)?shù)腒OG 簇，用于直系同源基因功能分類(lèi)。泡桐根部基因涉及的功能類(lèi)別較為全面，可將獲得到的29 886 條基因序列分為25 類(lèi)（圖2）。其中蛋白質(zhì)翻譯后修飾與轉(zhuǎn)運(yùn)、分子伴侶基因最多（Posttranslational modification, protein turnover,chaperones，3 854 條）；其次為一般功能基因（General function prediction only，3 619 條 ）；翻譯，核糖體結(jié)構(gòu)和生物合成基因（Translation,ribosomal structure and biogenesis，3 281 條），而細(xì)胞運(yùn)動(dòng)基因（Cell motility，29 條）和胞外結(jié)構(gòu)類(lèi)基因（Extracellular structures，60 條）較少。

2.2.4 KEGG 數(shù)據(jù)庫(kù)注釋分析與分類(lèi)

KEGG 是系統(tǒng)分析基因產(chǎn)物功能及其所在細(xì)胞代謝途徑的數(shù)據(jù)庫(kù)。據(jù)KEGG 數(shù)據(jù)庫(kù)的注釋信息進(jìn)一步將泡桐根部基因序列進(jìn)行pathway 注釋，共有28 747條基因序列富集于不同的代謝途徑中，其中包括細(xì)胞過(guò)程、環(huán)境信息處理、遺傳信息處理、代謝、有機(jī)系統(tǒng)5 個(gè)代謝通路大類(lèi)、環(huán)境適應(yīng)、翻譯、碳水化合物代謝等19 個(gè)代謝途徑及核糖體、氨基酸的生物合成、碳代謝等131 個(gè)代謝分支。在二級(jí)分類(lèi)中，富集于碳水化合物（Carbohydrate metabolism）（4 793 條）、翻譯（Transcription）（4 456 條）途徑的基因占大多數(shù)。在各分支中，富集于核糖體（Ribosome）、碳代謝（Carbon metabolism）、氨基酸生物合成（Biosynthesis of amino acids）的基因相對(duì)較多（表5）。以上這些代謝通路可能在泡桐根部生理代謝過(guò)程中具有重要作用。

圖1 GO 注釋分類(lèi)Fig. 1 GO annotation classification

2.2.5 基因表達(dá)水平分析

本實(shí)驗(yàn)采用RSEM 軟件[12]將 Trinity 拼接得到的轉(zhuǎn)錄組作為參考序列（ref），將每個(gè)樣品的clean reads 比對(duì)回ref（mapping），發(fā)現(xiàn)各樣品均有較高的注釋率（表6）。隨后將RSEM 對(duì)bowtie的比對(duì)結(jié)果進(jìn)行統(tǒng)計(jì)，進(jìn)一步得到了每個(gè)樣品比對(duì)到每個(gè)基因上的readcount 數(shù)目，并對(duì)其進(jìn)行FPKM 轉(zhuǎn)換，結(jié)果顯示，大部分樣品的FPKM 值位于0 ～0.1 區(qū)間和0.3 ～3.57 區(qū)間（表7）。

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

2.3.1 差異基因的篩選

本實(shí)驗(yàn)采用DESeq2[13]對(duì)在紅壤、鉛鋅礦中生長(zhǎng)的泡桐根部細(xì)胞轉(zhuǎn)錄組進(jìn)行基因差異表達(dá)分析，將在不同樣本間進(jìn)行兩兩比較時(shí)表達(dá)量變化閾值為 padj ＜ 0.05 且 |log2FoldChange|＞ 1 的 基 因 定義為差異表達(dá)基因（圖3）。在紅壤與鉛鋅礦之間，共篩選到4 297 個(gè)差異表達(dá)基因，其中4 143 個(gè)基因下調(diào)，154 個(gè)基因上調(diào)。

2.3.2 差異表達(dá)基因的GO 富集性分析

圖2 KOC 注釋分類(lèi)Fig. 2 KOC annotation classification diagram

為尋找響應(yīng)鉛鋅礦脅迫基因，本研究利用GOseq 比對(duì)紅壤和鉛鋅礦的差異基因發(fā)現(xiàn)，在鉛鋅礦中，大部分差異基因呈下調(diào)趨勢(shì)，極少一部分基因呈上調(diào)趨勢(shì)（圖4），說(shuō)明泡桐根部組織主要通過(guò)降低一些基因的表達(dá)量來(lái)應(yīng)對(duì)鉛鋅脅迫。結(jié)合（Binding）作為分子功能的主要類(lèi)別，富集其中的差異表達(dá)基因高達(dá)1 655 條，主要包括蛋白質(zhì)結(jié)合698 條，金屬離子結(jié)合381 條，過(guò)渡金屬離子結(jié)合250 條，鋅離子結(jié)合181 條等。在生物過(guò)程類(lèi)別中，細(xì)胞蛋白代謝（Cellular protein metabolic process）所含差異表達(dá)基因最多（482 條），另外，本研究發(fā)現(xiàn)細(xì)胞成分組織或生物發(fā)生（Cellular component organization or biogenesis）、細(xì)胞蛋白修飾（Cellular protein modification process）在生物過(guò)程類(lèi)別中含有相對(duì)較多的應(yīng)答脅迫差異基因。同時(shí)在細(xì)胞組分類(lèi)別中，細(xì)胞器（Organelle）、細(xì)胞內(nèi)細(xì)胞器（Intracellular organelle）、膜結(jié)合細(xì)胞器（Membrane-bounded organelle）所含差異基因數(shù)量位居前3 位。

表5 注釋Unigene 數(shù)量最多的10 個(gè)代謝通路Table 5 Ten metabolic pathways annotated by Unigene

表6 Reads 與參考序列比對(duì)Table 6 Reads and reference sequence alignment

2.3.3 差異表達(dá)基因KEGG 富集性分析

為揭示泡桐根部受到鉛鋅礦脅迫時(shí)差異表達(dá)基因代謝途徑和信號(hào)轉(zhuǎn)導(dǎo)途徑的富集性，本研究將差異表達(dá)基因序列在KEGG 數(shù)據(jù)庫(kù)中進(jìn)行富集性分析發(fā)現(xiàn)，所有差異表達(dá)基因共富集于108 個(gè)KEGG 代謝通路，其中富集差異表達(dá)基因數(shù)最多的7 個(gè)代謝通路依次為：核糖體（Ribosome）、內(nèi)質(zhì)網(wǎng)蛋白加工（Protein processing in endoplasmic reticulum）、剪接體（Spliceosome）、胞吞作用（Endocytosis）、RNA 轉(zhuǎn) 運(yùn)（RNA transport）、吞噬代謝（Phagosome）、mRNA 監(jiān)測(cè)（mRNA surveillance pathway）（表8）。這些通路涉及到泡桐根在鉛鋅脅迫下蛋白質(zhì)加工、遺傳信息加工、信號(hào)防御等方面，表明泡桐根通過(guò)調(diào)控生命活動(dòng)等各方面基因來(lái)應(yīng)答鉛鋅脅迫。

表7 FPKM 區(qū)間分布Table 7 FPKM interval distribution

圖3 樣品間基因差異表達(dá)分析火山Fig. 3 Volcanic maps of gene differential expression analysis between samples

核糖體在蛋白質(zhì)的翻譯過(guò)程中扮演著重要角色，代謝通路富集性分析發(fā)現(xiàn)此途徑中共包含111個(gè)差異表達(dá)基因，在鉛鋅礦脅迫下，108 條基因表達(dá)下調(diào)，3 條基因表達(dá)上調(diào)，所有下調(diào)基因均與核糖體蛋白合成有關(guān)，涉及到26 個(gè)小亞基核糖體蛋白，42 個(gè)大亞基核糖體蛋白，且有實(shí)驗(yàn)表明，單個(gè)核糖體蛋白減少會(huì)導(dǎo)致生物發(fā)育遲緩，活性降低[14]。另外在Ni 脅迫下，纖維細(xì)胞中與核糖體蛋白相關(guān)的基因表達(dá)呈下降趨勢(shì)[15]，說(shuō)明植物可能在重金屬脅迫下減少蛋白質(zhì)合成，減緩細(xì)胞生長(zhǎng)。

2.3.4 鉛鋅礦脅迫下的轉(zhuǎn)錄因子分析

在泡桐根轉(zhuǎn)錄組注釋信息中，一共鑒定得到37 個(gè)家族的3 506 個(gè)的轉(zhuǎn)錄因子（表9），通過(guò)對(duì)這些轉(zhuǎn)錄因子進(jìn)行深入分析，發(fā)現(xiàn)HSF，C2H2，MYB 家族是注釋信息較多的種類(lèi)，此外，還有一些諸如bZIP[16]、ARF、WRKY[17]等與環(huán)境相關(guān)的轉(zhuǎn)錄因子家族也得到鑒定。眾多的轉(zhuǎn)錄因子共同調(diào)控泡桐根部功能基因，形成完整網(wǎng)絡(luò)來(lái)應(yīng)答鉛鋅礦脅迫。

2.3.5 差異表達(dá)基因qRT-PCR 驗(yàn)證

為了驗(yàn)證轉(zhuǎn)錄組數(shù)據(jù)的準(zhǔn)確性，本研究隨機(jī)選擇6 個(gè)差異基因進(jìn)行qRT-PCR 驗(yàn)證，其對(duì)應(yīng)的引物如表10 所示。結(jié)果顯示，熒光定量PCR 與RNA-seq 檢測(cè)得到的基因相對(duì)表達(dá)量存在一定差異，但是兩者反應(yīng)的上下調(diào)趨勢(shì)一致（圖5），因此表明轉(zhuǎn)錄組測(cè)序得到的結(jié)果真實(shí)可靠。

3 結(jié)論與討論

新一代測(cè)序技術(shù)已經(jīng)應(yīng)用于許多無(wú)參考基因的植物包括紅花[18]、甘薯[19]等，目前，關(guān)于泡桐的遺傳背景研究較少，本研究應(yīng)用RNA-seq 高通量測(cè)序技術(shù)對(duì)鉛鋅礦生長(zhǎng)區(qū)和紅壤生長(zhǎng)區(qū)的泡桐根進(jìn)行轉(zhuǎn)錄組測(cè)序，共得到366 701 960 條reads，總堿基數(shù)為53.90 Gb，經(jīng)過(guò)濾得到259 921 130 條clean reads，且Q20、Q30 均達(dá)90%以上，說(shuō)明本次得到的測(cè)序數(shù)據(jù)質(zhì)量較好，再經(jīng)Trinity 拼接組裝得到126 061 個(gè)基因序列。這些數(shù)據(jù)揭示了泡桐根部轉(zhuǎn)錄組遺傳信息，為尋找響應(yīng)鉛鋅礦脅迫基因提供基礎(chǔ)。

將數(shù)字基因表達(dá)譜與高通量測(cè)序技術(shù)相結(jié)合是目前針對(duì)無(wú)參物種尋找差異基因最高效、最便捷的方法[20]，為尋找響應(yīng)鉛鋅礦脅迫基因，本研究以紅壤中泡桐根基因表達(dá)水平為參照來(lái)檢測(cè)鉛鋅礦脅迫下泡桐根中差異表達(dá)基因表達(dá)情況，發(fā)現(xiàn)在兩個(gè)樣品中共有4 297 個(gè)差異表達(dá)基因，其中4 143 個(gè)基因下調(diào)，154 個(gè)基因上調(diào)。說(shuō)明鉛鋅礦脅迫已經(jīng)影響了基因的正常表達(dá)水平，并且主要以下調(diào)基因表達(dá)模式來(lái)應(yīng)對(duì)重金屬脅迫。通過(guò)GO富集分析，在鉛鋅脅迫下，顯著富集下調(diào)基因主要涉及蛋白質(zhì)結(jié)合、金屬離子結(jié)合、過(guò)度金屬離子結(jié)合、鋅離子結(jié)合、細(xì)胞蛋白代謝、細(xì)胞蛋白修飾、細(xì)胞器等，其上調(diào)基因主要和結(jié)合有關(guān)。推測(cè)鉛鋅脅迫引起了泡桐根中與離子轉(zhuǎn)運(yùn)以及蛋白合成相關(guān)的基因發(fā)生變化，導(dǎo)致泡桐根部金屬離子穩(wěn)態(tài)以及蛋白合成發(fā)生變化。

圖4 鉛鋅礦和紅壤差異表達(dá)基因GO 分類(lèi)結(jié)果Fig. 4 Result of GO function classification on differentially expressed genes in lead-zinc mine and red soil

表8 差異基因數(shù)最多的KEGG termTable 8 KEGG term with the largest number of differentially expressed genes

根系作為重金屬深入植物的主要途徑，其對(duì)重金屬脅迫涉及到的代謝途徑是多種多樣的。有研究發(fā)現(xiàn)，當(dāng)蘿卜暴露于含鉛的環(huán)境中，鉛主要積累于根部，且可以通過(guò)細(xì)胞壁增厚來(lái)抵制鉛進(jìn)入原生質(zhì)體[21]。另外，根細(xì)胞壁被認(rèn)為是活躍的代謝位點(diǎn)，面對(duì)重金屬脅迫時(shí)，可以產(chǎn)生絲裂原活化蛋白激酶MAPKs（Mitogen activated protein kinase）和鈣結(jié)合蛋白CDPK（Calcium-binding related protein kinase）[22-23]。近年來(lái)，有研究發(fā)現(xiàn)，谷胱甘肽（GSH）可以通過(guò)抑制鉛誘導(dǎo)產(chǎn)生的ROS，來(lái)保護(hù)植物[24]?；比~蘋(píng)中的鉛積累可以增加谷胱甘肽數(shù)量，提高谷胱甘肽酶活性，且能提高根部SmGS 的基因表達(dá)量[25]。谷胱甘肽也可以對(duì)東南景天鉛脅迫進(jìn)行解毒[26]。金屬轉(zhuǎn)運(yùn)蛋白可以將金屬離子轉(zhuǎn)運(yùn)到細(xì)胞外或者螯合至液泡內(nèi)，這些轉(zhuǎn)運(yùn)過(guò)程往往涉及到精確基因的表達(dá)調(diào)控。目前發(fā)現(xiàn)的這些基因主要來(lái)自不同家族轉(zhuǎn)錄因子，例如一些ABC 轉(zhuǎn)運(yùn)蛋白可以直接參與運(yùn)輸離子或其他溶質(zhì)，甚至還參與調(diào)節(jié)通道。NARMP 可以對(duì)金屬二價(jià)陽(yáng)離子進(jìn)行轉(zhuǎn)運(yùn)[27]。本研究中，我們發(fā)現(xiàn)了大量與重金屬相關(guān)的轉(zhuǎn)錄因子如C2H2、MYB 等，這些轉(zhuǎn)錄因子在響應(yīng)并抵擋重金屬脅迫中起到重要作用。泡桐根部組織內(nèi)與重金屬脅迫的相關(guān)的代謝通路多種多樣，為了確定哪些代謝通路更有研究?jī)r(jià)值，我們對(duì)這些差異基因進(jìn)行KEGG 富集分析發(fā)現(xiàn)，這些代謝途徑主要涉及到核糖體、內(nèi)質(zhì)網(wǎng)蛋白加工、剪接體、胞吞作用、RNA 轉(zhuǎn)運(yùn)、mRNA 監(jiān)測(cè)、吞噬代謝等，可能這些代謝通路與重金屬積累、解毒有關(guān)，是以后研究的重點(diǎn)。

表9 轉(zhuǎn)錄因子分析Table 9 Analysis of transcription factors

表10 差異表達(dá)基因qRT-PCR 引物Table 10 qRT-PCR primers used in the validation of differentially expressed genes

圖5 差異表達(dá)基因的qRT-PCR 驗(yàn)證Fig. 5 Validation of differentially expressed genes using qRT-PCR

綜上所述，本研究通過(guò)高通量測(cè)序技術(shù)對(duì)生長(zhǎng)在鉛鋅礦和生長(zhǎng)在紅壤的的泡桐根部進(jìn)行轉(zhuǎn)錄組測(cè)序，獲得了質(zhì)量較好的轉(zhuǎn)錄本和基因，并對(duì)顯著差異表達(dá)基因進(jìn)行功能注釋，獲得了大量與重金屬積累及解毒的相關(guān)信息，為泡桐的抗重金屬生長(zhǎng)提供參考與幫助。