鄒承武， 宋 瑋， 宋慈安， 雷良奇

(1. 廣西大學(xué)農(nóng)學(xué)院， 廣西 南寧 530004； 2. 廣東工業(yè)大學(xué)計算機學(xué)院， 廣東 廣州 510006；3. 桂林理工大學(xué)： a. 廣西隱伏金屬礦產(chǎn)勘查重點實驗室， b. 地球科學(xué)學(xué)院， 廣西 桂林 541004)

基于RNA-Seq技術(shù)分析在金屬礦上部生長的芒萁的差異表達基因

鄒承武1， 宋 瑋2， 宋慈安3a，3b，①， 雷良奇3a，3b

(1. 廣西大學(xué)農(nóng)學(xué)院， 廣西 南寧 530004； 2. 廣東工業(yè)大學(xué)計算機學(xué)院， 廣東 廣州 510006；3. 桂林理工大學(xué)： a. 廣西隱伏金屬礦產(chǎn)勘查重點實驗室， b. 地球科學(xué)學(xué)院， 廣西 桂林 541004)

以生長于廣西大廠錫多金屬礦上部(重金屬脅迫區(qū))和未受礦化或污染影響的礦區(qū)外圍(對照區(qū))的芒萁〔Dicranopterispedata(Houtt.) Nakaike〕為實驗材料，對芒萁葉片進行轉(zhuǎn)錄組高通量測序，并對組裝得到的unigenes經(jīng)NCBI官方非冗余蛋白質(zhì)序列數(shù)據(jù)庫(Nr)、NCBI官方非冗余核苷酸序列數(shù)據(jù)庫(Nt)、KEGG直系同源數(shù)據(jù)庫(KO)、Swiss-Prot數(shù)據(jù)庫(Swiss-Prot)、蛋白質(zhì)家族數(shù)據(jù)庫(Pfam)、基因功能分類體系數(shù)據(jù)庫(GO)和真核生物直系同源序列數(shù)據(jù)庫(KOG)進行注釋，同時分析重金屬脅迫區(qū)和對照區(qū)芒萁葉片間的差異表達unigenes。結(jié)果顯示：測序獲得19.56 Gb clean data，其中，重金屬脅迫區(qū)和對照區(qū)芒萁葉片分別含10.14和9.42 Gb clean data。組裝得到的250 582個unigenes中有120 097個unigenes得到注釋，占unigenes總數(shù)的47.93%。與對照區(qū)相比較，重金屬脅迫區(qū)芒萁葉片中上調(diào)和下調(diào)差異表達unigenes分別有208和620個，其中120個上調(diào)差異表達unigenes注釋為代謝過程，占所有上調(diào)差異表達unigenes的57.69%；285個下調(diào)差異表達unigenes注釋為催化活性，占所有下調(diào)差異表達unigenes的45.97%。重金屬脅迫區(qū)芒萁葉片中15個unigenes與重金屬轉(zhuǎn)運和耐受相關(guān)，其中c44988_g1和c84121_g1的相對表達量分別極顯著和顯著高于對照區(qū)。研究結(jié)果顯示：芒萁響應(yīng)自然金屬礦化或礦山重金屬污染的基因可以用于生物地球化學(xué)找礦和土壤重金屬污染檢測。

芒萁； RNA-Seq； 重金屬脅迫； 差異表達基因

芒萁〔Dicranopterispedata(Houtt.) Nakaike〕隸屬于里白科(Gleicheniaceae)芒萁屬(DicranopterisBernh.)，是一種普遍生長于熱帶、亞熱帶紅壤丘陵、荒坡林緣的蕨類植物，耐酸、耐旱、耐貧瘠，屬于典型的酸性土壤指示植物[1]。在金、錫、鉬和鉈等金屬礦床上部生長的芒萁，其體內(nèi)與成礦有關(guān)的元素明顯高于背景區(qū)，產(chǎn)生了生物地球化學(xué)異常，因此可以作為生物地球化學(xué)找礦的有效指示植物之一[2-6]。在受礦區(qū)開采污染土壤中生長的芒萁，其地上部分能吸收大量重金屬元素，屬于金屬富集型植物，可用于污染土壤的植物修復(fù)[7-8]。說明芒萁受到高含量金屬環(huán)境的脅迫后，體內(nèi)會聚集大量的金屬元素。

大量研究結(jié)果表明:植物體內(nèi)聚集的金屬元素特別是重金屬元素，可顯著影響植物細胞內(nèi)的基因表達水平[9-13]。為了揭示芒萁受金屬脅迫后的轉(zhuǎn)錄調(diào)控機制，明確芒萁響應(yīng)金屬脅迫的相關(guān)基因和重要通路，作者分別采集了生長在廣西大廠錫多金屬礦上部和未受礦化或污染影響的礦區(qū)外圍的芒萁植株作為實驗材料，進行轉(zhuǎn)錄組測序并從頭組裝(denovoassembly)，分析芒萁體內(nèi)響應(yīng)重金屬脅迫的相關(guān)差異表達基因，以期進一步研究受重金屬脅迫差異表達的基因與自然礦床或開采礦山污染土壤的關(guān)系，并為芒萁應(yīng)用于生物地球化學(xué)找礦和植物修復(fù)環(huán)境治理提供基礎(chǔ)理論依據(jù)。

1 材料和方法

1.1 材料

供試芒萁植株分別采自廣西大廠錫多金屬礦上部(以下簡稱為重金屬脅迫區(qū))土壤和未受礦化或污染影響的礦區(qū)外圍(以下簡稱為對照區(qū))土壤。其中，重金屬脅迫區(qū)生長的芒萁植株較對照區(qū)矮小，葉面泛黃伴褐色銹斑。

隨機選取重金屬脅迫區(qū)和對照區(qū)同一生長時期的芒萁植物各3株,剪取同一部位、大小一致的新鮮葉片，同一區(qū)域的芒萁葉片各稱取0.1 g，剪碎并混合均勻，置于5 mL凍存管中，立即投入液氮中保存，用于總RNA的提取。

在采集芒萁葉片的同時，采集芒萁植株根際土壤樣品并測定土壤中相關(guān)礦化元素含量。測定結(jié)果表明：重金屬脅迫區(qū)土壤中砷、銻、鎢、鉍、錫、鉛、銀和鎳的含量分別為對照區(qū)土壤的39、38、33、27、26、16、11和11倍。

1.2 方法

1.2.1 總RNA提取及檢測方法 根據(jù)TRIzol?Reagent試劑(美國Thermo Fisher Scientific公司)產(chǎn)品說明書進行實驗操作，提取芒萁葉片總RNA。每100 mg 芒萁葉片加入1 mL TRIzol?Reagent，迅速混勻，加入0.2 mL三氯甲烷，劇烈震蕩30 s，室溫靜置5 min，4 ℃條件下 12 000 r·min-1離心15 min。為獲得高質(zhì)量的總RNA，避免觸碰中間雜蛋白層，小心地吸取上層水相約0.6 mL于1.5 mL Eppendorf管中，加入10 μL DNase Ⅰ(RNase free)消化DNA，輕柔混勻30 s，室溫孵育10 min，加入等體積三氯甲烷，劇烈震蕩30 s，室溫靜置5 min，4 ℃ 條件下12 000 r·min-1離心15 min。小心吸取上層水相0.4 mL，加入等體積異丙醇，室溫放置10 min，4 ℃條件下12 000 r·min-1離心10 min，然后用體積分數(shù)75%乙醇清洗沉淀2次，4 ℃條件下8 000 r·min-1離心1 min，棄乙醇，冰上自然干燥沉淀，最后用RNase-free水溶解沉淀。利用Agilent 2100生物分析儀(美國Agilent公司)對總RNA的質(zhì)量和濃度進行檢測和分析。

1.2.2 cDNA文庫構(gòu)建及RNA-Seq測序方法 樣品檢測合格后，分別取等體積的經(jīng)DNase Ⅰ消化的芒萁葉片總RNA，用帶有Oligo(dT)的磁珠(美國Thermo Fisher Scientific公司)富集poly(A) mRNA，隨后加入5×Fragmentation buffer(美國Illumina公司)在ThermoMixer?F1.5恒溫混勻儀(德國Eppendorf公司)中于室溫條件下將mRNA打斷成短片段，以打斷的mRNA為模板，用六堿基隨機引物(random hexamers)合成第1鏈cDNA，然后加入緩沖液、dNTPs、DNA polymerase Ⅰ和RNase H合成第2鏈cDNA，再用AMPure XP Beads試劑盒(美國Illumina公司)純化雙鏈cDNA。純化后的雙鏈cDNA先進行末端修復(fù)、加A尾并連接測序接頭，再用AMPure XP Beads試劑盒進行片段大小選擇，最后使用ABI 9700型PCR擴增儀(美國Applied Biosystems公司)進行PCR擴增。PCR反應(yīng)體系總體積為50.0 μL，包含10.0 μL 5×Phusion buffer、1.0 μL PCR Primer PE 1.0、1.0 μL PCR Primer PE 2.0、0.5 μL Phusion DNA Polymerase、0.5 μL 25 mmol·L-1dNTPs Mix、7.0 μL ddH2O和30.0 μL已經(jīng)連接測序接頭的連接產(chǎn)物。擴增程序為：98 ℃預(yù)變性30 s；98 ℃變性10 s，65 ℃退火30 s，72 ℃延伸30 s，15個循環(huán)；最后在72 ℃延伸5 min。擴增反應(yīng)結(jié)束后降溫至4 ℃，用AMPure XP Beads試劑盒純化PCR產(chǎn)物，得到最終的cDNA文庫。

cDNA文庫構(gòu)建完成后，先使用Qubit?2.0熒光定量儀(美國Thermo Fisher Scientific公司)進行初步定量，稀釋至1.5 ng·μL-1，隨后使用Agilent 2100型生物分析儀(美國Agilent 公司)采用DNA 1000分析試劑盒(美國Agilent公司)對cDNA文庫的insert size進行檢測，insert size符合預(yù)期后，使用Illumina qPCR定量操作指南中的qPCR方法對構(gòu)建的cDNA文庫的有效濃度進行準確定量(cDNA文庫有效濃度大于2 nmol·L-1)，以保證cDNA文庫的質(zhì)量。使用Illumina HiSeqTM4000平臺(美國Illumina公司)對cDNA文庫進行測序，測序在北京諾禾致源生物信息科技有限公司完成。

1.2.3 Unigene序列組裝 使用組裝軟件Trinity[14]將Illumina高通量測序得到的clean reads按順序進行denovo組裝獲得unigene。

1.2.4 Unigene序列注釋、功能分類及生物學(xué)通路分析 使用BLAST軟件利用NCBI官方非冗余蛋白質(zhì)序列數(shù)據(jù)庫(NCBI non-redundant protein sequence database，Nr)、NCBI官方非冗余核苷酸序列數(shù)據(jù)庫(NCBI non-redundant nucleotide sequence database，Nt)、KEGG直系同源數(shù)據(jù)庫(KEGG Orthology database，KO)、Swiss-Prot數(shù)據(jù)庫(Swiss-Prot)、蛋白質(zhì)家族數(shù)據(jù)庫(Protein family database，Pfam)、基因功能分類體系數(shù)據(jù)庫(Gene Ontology database，GO)和真核生物直系同源序列數(shù)據(jù)庫(euKaryotic Orthology Groups database，KOG)等對unigene進行注釋，同時使用Blast2GO[15]以及Nr注釋結(jié)果進行GO注釋和GO功能分類統(tǒng)計，再聯(lián)合使用InterProScan 5[16]進行InterPro注釋分析，從而得到unigene的功能注釋信息。

1.2.5 差異表達unigene分析 使用RSEM(RNA-Seq by expectation-maximization)[17]將clean data中的reads映射(mapping)到通過Trinity軟件組裝獲得的轉(zhuǎn)錄物上，Bowtie 2[18]不允許錯配，使用DEGseq[19]進行差異表達unigene檢測。所得到的差異表達unigene分別使用goseq軟件[20]和GO富集軟件[21]進行分析。

1.2.6 差異表達unigene的驗證 選取轉(zhuǎn)錄組測序獲得的差異表達unigene，以芒萁的actin基因作為內(nèi)參基因，根據(jù)基因序列使用Primer 6.0軟件設(shè)計引物并由生工生物工程(上海)股份有限公司合成。使用轉(zhuǎn)錄組測序的RNA樣品進行反轉(zhuǎn)錄獲得第1鏈cDNA，然后通過熒光定量PCR法對選取unigene的表達量進行定量分析。熒光定量PCR實驗使用羅氏熒光定量試劑盒(FastStart Universal SYBR Green Master，瑞士Roche公司)并參考說明書進行操作，并在羅氏LightCycler? 480Ⅱ?qū)崟r熒光定量PCR系統(tǒng)(瑞士Roche公司)進行反應(yīng)。擴增程序為：95 ℃預(yù)變性10 min；95 ℃變性5 s，72 ℃退火40 s，45個循環(huán)。擴增反應(yīng)結(jié)束后，以芒萁的actin基因進行校正，并將對照區(qū)芒萁葉片中待測unigene的轉(zhuǎn)錄水平設(shè)為“1”，重金屬脅迫區(qū)芒萁葉片中待測unigene的相對表達量采用2-ΔΔCT方法[22]進行計算。

2 結(jié)果和分析

2.1 轉(zhuǎn)錄組測序概況及序列組裝結(jié)果分析

用Illumina HiSeqTM4000平臺對構(gòu)建的芒萁葉片cDNA文庫進行高通量測序，共獲得135 646 166條raw reads，對其進行數(shù)據(jù)質(zhì)控，丟棄5 259 234條不合格的reads，獲得測序數(shù)據(jù)19.59 Gb clean data，其中重金屬脅迫區(qū)芒萁葉片含10.14 Gb clean data，對照區(qū)芒萁葉片含9.42 Gb clean data，平均Q30分別為92.94%和93.29%。

利用組裝軟件Trinity[14]對clean reads進行混合拼裝，獲得的轉(zhuǎn)錄物平均長度為567 bp，最大長度為17 363 bp，選取拼接得到的最長轉(zhuǎn)錄物作為unigenes，獲得250 582個unigenes，平均長度為455 bp，N50長度為515 bp，轉(zhuǎn)錄物與unigenes長度分布見圖1。由圖1可見：長度為200～300 bp的轉(zhuǎn)錄物和unigenes較多，分別有158 839和149 360個；長度為301～500 bp的轉(zhuǎn)錄物和unigenes明顯減少，分別有63 027和54 594個；長度為501～1 000 bp的轉(zhuǎn)錄物和unigenes分別有35 444和26 179個；長度為1 001～2 000 bp的轉(zhuǎn)錄物和unigenes分別有24 874和13 266個；長度大于2 000 bp的轉(zhuǎn)錄物和unigenes分別有15 636和7 183個。

□： 轉(zhuǎn)錄物Transcript； ■： Unigene.A： 200-300； B： 301-500； C： 501-1 000； D： 1 001-2 000； E： 大于2 000 Above 2 000.圖1 組裝獲得芒萁葉片中轉(zhuǎn)錄物與unigenes的長度分布圖Fig. 1 Distribution diagram of length of assembled transcripts and unigenes in leaves of Dicranopteris pedata (Houtt.) Nakaike

2.2 Unigenes功能注釋分析

為了獲得芒萁葉片中unigenes的功能信息，將組裝得到的250 582個unigenes進行NCBI官方非冗余蛋白質(zhì)序列數(shù)據(jù)庫(Nr)、NCBI官方非冗余核苷酸序列數(shù)據(jù)庫(Nt)、KEGG直系同源數(shù)據(jù)庫(KO)、Swiss-Prot數(shù)據(jù)庫(Swiss-Prot)、蛋白質(zhì)家族數(shù)據(jù)庫(Pfam)、基因功能注釋數(shù)據(jù)庫(GO)和真核生物直系同源序列數(shù)據(jù)庫(KOG)的BLAST注釋，各數(shù)據(jù)庫分別注釋了80 536、30 096、37 632、85 035、81 450、83 611和52 557個unigenes，分別占unigenes總數(shù)的32.14%、12.01%、15.02%、33.93%、32.50%、33.37%和20.97%。統(tǒng)計結(jié)果顯示：共有120 097個unigenes得到注釋，占unigenes總數(shù)的47.93%；有130 485個unigenes未被注釋，占unigenes總數(shù)的52.07%。

GO數(shù)據(jù)庫對基因及其產(chǎn)物的描述分為生物過程(biological process)、細胞組分(cellular component)和分子功能(molecular function)。對獲得GO注釋的83 611個unigenes進行分類統(tǒng)計(圖2)。在生物過程中，43 673個unigenes注釋為細胞過程(cellular process)，44 569個unigenes注釋為代謝過程(metabolic process)。在分子功能中，38 014個unigenes注釋為催化活性(catalytic activity)，38 333個unigenes注釋為結(jié)合(binding)。

KEGG是一系列分子作用與互作通路網(wǎng)絡(luò)，包括新陳代謝(metabolism)、遺傳信息處理(genetic information processing)、環(huán)境信息處理(environmental information processing)、細胞過程(cellular processes)、生物系統(tǒng)(organismal systems)和人類疾病(human diseases)。為了進一步研究芒萁中活躍的生物學(xué)通路，本研究使用KAAS(KEGG Automatic Annotation Server)[23]將37 632個unigenes定位到KO數(shù)據(jù)庫中除人類疾病以外的5種分類(圖3)。在二級分類中較為豐富的通路分別為碳水化合物代謝(carbohydrate metabolism)與翻譯(translation)，分別有6 007和5 282個unigenes參與通路，分別占KO數(shù)據(jù)庫注釋unigenes總數(shù)的15.96%和14.04%；氨基酸代謝(amino acid metabolism)通路也較豐富，有4 095個unigenes參與通路，占KO數(shù)據(jù)庫注釋unigenes總數(shù)的10.88%。這些通路可能在芒萁生化代謝中具有重要作用。

2.3 差異表達unigenes檢測

圖2 組裝獲得的芒萁葉片中unigenes的GO功能分類結(jié)果Fig. 2 Result of GO function classification of assembled unigenes in leaves of Dicranopteris pedata (Houtt.) Nakaike

圖中數(shù)據(jù)表示unigenes數(shù)量Datums in the diagram indicate number of unigenes. A： 新陳代謝Metabolism； B： 生物系統(tǒng)Organismal systems； C： 遺傳信息處理Genetic information processing； D： 環(huán)境信息處理Environmental information processing； E： 細胞過程Cellular processes.圖3 組裝獲得的芒萁葉片中unigenes的KO代謝通路分類結(jié)果Fig. 3 Result of KO metabolic pathway classification of assembled unigenes in leaves of Dicranopteris pedata (Houtt.) Nakaike

以對照區(qū)芒萁葉片中unigenes表達水平為參照，檢測重金屬脅迫區(qū)芒萁葉片中差異表達unigenes的表達情況，選取q<0.005并且log2(fold change)>1的差異表達unigenes評估重金屬脅迫區(qū)和對照區(qū)芒萁葉片中的差異表達unigenes，結(jié)果見圖4。在重金屬脅迫區(qū)芒萁葉片中上調(diào)差異表達unigenes有208個，下調(diào)差異表達unigenes有620個，下調(diào)差異表達unigenes是上調(diào)差異表達unigenes的2.98倍。

2.4 差異表達unigenes的GO功能分類分析

對重金屬脅迫區(qū)芒萁葉片中上調(diào)和下調(diào)差異表達unigenes進行GO功能分類富集分析，結(jié)果分別見圖5和圖6。上調(diào)差異表達unigenes中有28個富集的GO條目，其中有120個差異表達unigenes注釋為代謝過程(metabolic process)，占所有上調(diào)差異表達unigenes的57.69%。下調(diào)差異表達unigenes中有18個富集的GO條目，下調(diào)差異表達unigenes中有285個注釋為催化活性(catalytic activity)，占所有下調(diào)差異表達unigenes的45.97%。

2.5 與重金屬轉(zhuǎn)運和耐受相關(guān)的unigenes及其差異表達情況

對重金屬脅迫區(qū)芒萁葉片中組裝獲得的unigenes進行BLASTp注釋，結(jié)果見表1。結(jié)果顯示：15個unigenes與重金屬轉(zhuǎn)運和耐受相關(guān)。

padj： 調(diào)整后的p值A(chǔ)djusted p-value. ： 下調(diào)差異表達unigenes Down-regulated differentially expressed unigenes； ： 除上調(diào)和下調(diào)以外的差異表達unigenes Differentially expressed unigenes without up- and down-regulation； ： 上調(diào)差異表達unigenes Up-regulated differentially expressed unigenes. 橫向虛線以下及縱向虛線間為無顯著性的差異表達unigenes Differentially expressed unigenes below horizontal dashed line and between vertical dashed lines without significance.圖4 重金屬脅迫區(qū)芒萁葉片中差異表達unigenes的火山圖Fig. 4 Volcano plot of differentially expressed unigenes in leaves of Dicranopteris pedata (Houtt.) Nakaike in heavy metal stress area

圖中不包括被歸類為“未知”功能的unigenes The unigenes classified as function “unknown” are not included in the figure.圖5 重金屬脅迫區(qū)芒萁葉片中上調(diào)差異表達unigenes的GO功能分類結(jié)果Fig. 5 Result of GO function classification on up-regulated differentially expressed unigenes in leaves of Dicranopteris pedata (Houtt.) Nakaike in heavy metal stress area

圖中不包括被歸類為“未知”功能的unigenes The unigenes classified as function “unknown” were not included in the figure.圖6 重金屬脅迫區(qū)芒萁葉片中下調(diào)差異表達unigenes的GO功能分類結(jié)果Fig. 6 Result of GO function classification on down-regulated differentially expressed unigenes in leaves of Dicranopteris pedata (Houtt.) Nakaike in heavy metal stress area

表1 重金屬脅迫區(qū)芒萁葉片中BLASTp注釋的與重金屬轉(zhuǎn)運和耐受相關(guān)的unigenes

Table 1 Unigenes descripted by BLASTp and associated with heavy metal transport and tolerance in leaves ofDicranopterispedata(Houtt.) Nakaike in heavy metal stress area

Unigene序列長度/bpLengthofsequence序列號Accessionnumber期望值Expectedvalue注釋Description表達水平1)Regulation1)c13489＿g1519AIG13049 12 90E-12假定重金屬轉(zhuǎn)運P1B-ATPase2Putativeheavymetaltrans?portingP1B?ATPase2URc93031＿g1842KEH18420 14 68E-10重金屬相關(guān)結(jié)構(gòu)域蛋白Heavymetal?associateddomainproteinURc91259＿g11516KEH36986 14 17E-24重金屬轉(zhuǎn)運/解毒超家族蛋白Heavymetaltransport/detox?ificationsuperfamilyproteinURc54437＿g1509KGN46770 16 93E-57假定蛋白Csa＿6g133790HypotheticalproteinCsa＿6G133790DRc107663＿g1237NP＿001052304 12 08E-29重金屬相關(guān)異戊二烯植物蛋白26Heavymetal?associatedisoprenylatedplantprotein26S c11915＿g1722XP＿001697267 14 25E-70一半大小的ABC轉(zhuǎn)運蛋白，膜蛋白Half?sizeABCtrans?porter，membraneproteinURc132091＿g1447XP＿001703119 17 62E-32線粒體上的一半大小的ABC轉(zhuǎn)運蛋白，膜蛋白Mitochon?drialhalf?sizeABCtransporter，membraneproteinDRc13999＿g1297XP＿005649995 11 19E-11重金屬遷移HeavymetaltranslocationURc207040＿g1295XP＿005829762 14 15E-11重金屬輸出蛋白hMT1，ABC超家族HeavymetalexporterHMT1，ABCsuperfamilyDRc44988＿g11093XP＿006448292 11 44E-16重金屬輸出蛋白hMT1，ABC超家族HeavymetalexporterHMT1，ABCsuperfamilyDRc92284＿g11026XP＿006846435 14 41E-35重金屬相關(guān)異戊二烯植物蛋白26Heavymetal?associatedisoprenylatedplantprotein26DRc84121＿g11318XP＿008783549 12 91E-34重金屬相關(guān)異戊二烯植物蛋白26Heavymetal?associatedisoprenylatedplantprotein26URc177177＿g1241XP＿009348170 15 31E-12重金屬耐受樣蛋白Heavymetaltoleranceprotein?likeURc115986＿g1252XP＿010525447 11 82E-26重金屬相關(guān)異戊二烯植物蛋白26Heavymetal?associatedisoprenylatedplantprotein26DRc237853＿g1246XP＿011399202 15 66E-18線粒體ATP結(jié)合盒亞家族B成員6ATP?bindingcassettesub?familyBmember6，mitochondrialUR

1)UR： 上調(diào)Up-regulation； S： 相似Similar； DR： 下調(diào)Down-regulation.

從15個與重金屬轉(zhuǎn)運和耐受相關(guān)的unigenes中選擇5個組裝獲得的序列長度較長的unigenes，分別為c93031_g1、c91259_g1、c44988_g1、c92284_g1和c84121_g1。RNA-Seq檢測結(jié)果(圖7-A)顯示：c44988_g1在重金屬脅迫區(qū)芒萁葉片中的相對表達量只有對照區(qū)的37%，與對照區(qū)差異極顯著(q<0.001)；而c84121_g1在重金屬脅迫區(qū)芒萁葉片中的相對表達量為對照區(qū)的3.2倍，與對照區(qū)差異顯著(q<0.005)；其余3個與重金屬轉(zhuǎn)運和耐受相關(guān)的unigenes的相對表達量在重金屬脅迫區(qū)和對照區(qū)間的差異均不顯著。

以芒萁的actin基因為內(nèi)參基因，對5個unigenes的相對表達量進行熒光定量PCR檢測(圖7-B)。結(jié)果顯示：熒光定量PCR和RNA-Seq檢測得到的5個unigenes的相對表達量存在一定差異，但是二者反映出的上調(diào)和下調(diào)的趨勢一致。

□： 對照區(qū)The control area； ■： 重金屬脅迫區(qū)Heavy metal stress area. U1: c44988-g1; U2: c84121-g1; U3: c91259-g1; U4: c92284-g1; U5: c93031-g1.A： RNA-Seq檢測RNA-Seq detection； B： 熒光定量PCR檢測Fluorescence quantitative PCR detection.圖7 重金屬脅迫區(qū)芒萁葉片中與重金屬轉(zhuǎn)運和耐受相關(guān)的unigenes的差異表達分析Fig. 7 Differentially expressed analysis on unigenes associated with heavy metal transport and tolerance in leaves of Dicranopteris pedata (Houtt.) Nakaike in heavy metal stress area

3 討 論

不同種類金屬脅迫可導(dǎo)致植物某些基因表達水平出現(xiàn)不同的變化[9-13，24]，這些基因可能潛在指示某種金屬的脅迫作用，或者同種金屬不同濃度的脅迫作用，找到響應(yīng)金屬脅迫的差異表達基因有望為利用植物找礦開辟新途徑。用常規(guī)的分子生物學(xué)手段尋找響應(yīng)金屬脅迫的差異表達基因不僅要求研究對象有基因組數(shù)據(jù)支持，而且工作量大，耗資巨大。RNA-Seq高通量測序技術(shù)，不但可以快速全面地了解樣本間的差異表達基因，而且測序價格越來越低，為將來利用該技術(shù)尋找指示礦藏的植物基因提供了可能。

芒萁的遺傳背景研究較少，在NCBI上尚未找到已經(jīng)公開報道的芒萁參考基因組或轉(zhuǎn)錄組信息。本研究應(yīng)用RNA-Seq高通量測序技術(shù)對重金屬脅迫區(qū)和對照區(qū)的芒萁進行轉(zhuǎn)錄組測序，獲得19.56 Gb clean data，其中重金屬脅迫區(qū)和對照區(qū)芒萁葉片中分別含10.14和9.42 Gb clean data，獲得的芒萁轉(zhuǎn)錄組信息不但能用于尋找響應(yīng)金屬脅迫相關(guān)基因，還為芒萁遺傳信息的相關(guān)研究提供了數(shù)據(jù)參考。

將轉(zhuǎn)錄組測序分析結(jié)果與數(shù)字基因表達譜相結(jié)合是篩選、分析無參考基因組物種差異表達基因最有效的方法之一。本研究設(shè)定對照區(qū)芒萁葉片中的基因表達水平為參照來檢測重金屬脅迫區(qū)芒萁葉片中差異表達unigenes的表達情況，重金屬脅迫區(qū)芒萁樣品中上調(diào)差異表達unigenes有208個，下調(diào)差異表達unigenes有620個，說明重金屬脅迫影響了基因的正常表達水平。推測芒萁能富集和耐受高濃度的重金屬元素可能與某些重金屬轉(zhuǎn)運和耐受相關(guān)的基因調(diào)控有關(guān)。本研究對重金屬脅迫區(qū)芒萁葉片中的unigenes進行BLASTp注釋，明確了這些unigenes在細胞代謝過程的具體分類和生物通路，共有15個unigenes被注釋為與重金屬轉(zhuǎn)運和耐受相關(guān)，其中c44988_g1和c84121_g1的相對表達量分別極顯著和顯著高于對照區(qū)，這2個unigenes可以作為潛在的用于植物地球化學(xué)找礦的候選基因。

本研究中，芒萁的轉(zhuǎn)錄物中僅47.93%的unigenes被注釋，可能由于芒萁的基因組和轉(zhuǎn)錄組還沒有相關(guān)研究，在數(shù)據(jù)庫中沒有相應(yīng)的參考序列，而芒萁作為比較古老的植物，與很多已進行基因組測序的物種親緣關(guān)系較遠，所以根據(jù)其他物種的基因信息進行注釋獲得的注釋基因數(shù)量較少。此外，由于芒萁轉(zhuǎn)錄物可能比較復(fù)雜，常規(guī)組裝手段獲得的unigenes平均長度太短，數(shù)量太大，影響后續(xù)的數(shù)據(jù)分析，今后需要進一步優(yōu)化組裝的參數(shù)，或采用第三代測序技術(shù)，提高測序讀長，從而使更多的芒萁基因得到注釋。對未能獲得注釋的響應(yīng)重金屬脅迫的unigenes也可以開展基因功能研究，以期獲得能夠穩(wěn)定響應(yīng)某種重金屬脅迫的基因用于植物地球化學(xué)找礦和土壤重金屬污染檢測。

綜上所述，生長于重金屬脅迫區(qū)的芒萁與對照區(qū)的unigenes的表達差異較大，這些差異表達unigenes可以作為植物地球化學(xué)找礦和土壤重金屬污染檢測的候選基因，但還需要后續(xù)的研究進一步驗證。

[1] 陸樹剛. 蕨類植物學(xué)[M]. 北京： 高等教育出版社， 2007： 71-74.

[2] 季峻峰， 崔衛(wèi)東， 孫承轅. 湖南黃金洞金礦床植物地球化學(xué)勘查的初步研究[J]. 物探與化探， 1992， 16(6)： 470-473.

[3] 馬躍良. 廣東省河臺金礦生物地球化學(xué)特征及遙感找礦意義[J]. 礦物學(xué)報， 2000， 20(1)： 80-86.

[4] 陳代演， 鄒振西， 任大銀. 植物找礦法在尋找鉈礦床中的初步應(yīng)用[J]. 礦物巖石地球化學(xué)通報， 2000， 19(4)： 397-400.

[5] 徐金鴻， 徐瑞松， 夏 斌. 廣東鼎湖山斑巖鉬礦區(qū)生物地球化學(xué)特征[J]. 地球與環(huán)境， 2006， 34(1)： 23-28.

[6] 陳 楊， 蒙夢平， 宋慈安. 植物對土壤中微量元素的吸收與轉(zhuǎn)移及對生物地球化學(xué)異常形成的影響： 以廣西鹽田嶺錫石硫化物礦床為例[J]. 地球與環(huán)境， 2012， 40(2)： 208-218.

[7] 劉足根， 楊國華， 楊 帆， 等. 贛南鎢礦區(qū)土壤重金屬含量與植物富集特征[J]. 生態(tài)學(xué)雜志， 2008， 27(8)： 1345-1350.

[8] 楊勝香， 田啟建， 梁士楚， 等. 湘西花垣礦區(qū)主要植物種類及優(yōu)勢植物重金屬蓄積特征[J]. 環(huán)境科學(xué)， 2012， 33(6)： 2038-2045.

[9] LOUIE M， KONDOR N， DEWITT J G. Gene expression in cadmium-tolerantDaturainnoxia： detection and characterization of cDNAs induced in response to Cd2+[J]. Plant Molecular Biology， 2003， 52： 81-89.

[10] ROUT J R， SAHOO S L. Antioxidant enzyme gene expression in response to copper stress inWithaniasomniferaL.[J]. Plant Growth Regulation， 2013， 71： 95-99.

[11] HUANG B， XIN J， YANG Z， et al. Suppression subtractive hybridization (SSH)-based method for estimating Cd-induced differences in gene expression at cultivar level and identification of genes induced by Cd in two water spinach cultivars[J]. Journal of Agricultural and Food Chemistry， 2009， 57： 8950-8962.

[12] TADDEI S， BERNARDI R， SALVINI M， et al. Effect of copper on callus growth and gene expression ofinvitro-cultured pith explants ofNicotianaglauca[J]. Plant Biosystems， 2007， 141： 194-203.

[13] NEVRTALOVA E， BALOUN J， HUDZIECZEK V， et al. Expression response of duplicatedmetallothionein3 gene to copper stress inSilenevulgarisecotypes[J]. Protoplasma， 2014， 251： 1427-1439.

[14] GRABHERR M G， HAAS B J， YASSOUR M， et al. Trinity： reconstructing a full-length transcriptome without a genome from RNA-Seq data[J]. Nature Biotechnology， 2013， 29： 644-652.

[15] CONESA A， G?TZ S. Blast2GO： a comprehensive suite for functional analysis in plant genomics[J]. International Journal of Plant Genomics， 2008(2008)： 619832.

[16] JONES P， BINNS D， CHANG H Y， et al. InterProScan 5： genome-scale protein function classification[J]. Bioinformatics， 2014， 30： 1236-1240.

[17] LI B， DEWEY C N. RSEM： accurate transcript quantification from RNA-Seq data with or without a reference genome[J]. BMC Bioinformatics， 2011， 12： 323.

[18] LANGMEAD B， SALZBERG S L. Fast gapped-read alignment with Bowtie 2[J]. Nature Methods， 2012， 9： 357-359.

[19] WANG L， FENG Z， WANG X， et al. DEGseq： an R package for identifying differentially expressed genes from RNA-seq data[J]. Bioinformatics， 2010， 26： 136-138.

[20] YOUNG M D， WAKEFIELD M J， SMYTH G K， et al. goseq： Gene Ontology testing for RNA-seq datasets[J]. Nature Methods， 2012， 9： 357-359.

[21] BOYLE E I， WENG S， GOLLUB J， et al. GO∶∶TermFinder： open source software for accessing Gene Ontology information and finding significantly enriched Gene Ontology terms associated with a list of genes[J]. Bioinformatics， 2004， 20： 3710-3715.

[22] LIVAK K J， SCHMITTGEN T D. Analysis of relative gene expression data using real-time quantitative PCR and the 2-ΔΔCTmethod[J]. Methods， 2001， 25： 402-408.

[23] MORIYA Y， ITOH M， OKUDA S， et al. KAAS： an automatic genome annotation and pathway reconstruction server[J]. Nucleic Acids Research， 2007， 35： W182-W185.

[24] GHELICH S， ZARINKAMAR F， SOLTANI B M， et al. Effect of lead treatment on medicarpin accumulation and on the gene expression of key enzymes involved in medicarpin biosynthesis inMedicagosativaL[J]. Environmental Science and Pollution Research， 2014， 21： 14091-14098.

(責任編輯： 張明霞)

Analysis on differentially expressed genes inDicranopterispedatagrown on metal deposit based on RNA-Seq technology

ZOU Chengwu1， SONG Wei2， SONG Cian3a，3b，①， LEI Liangqi3a，3b

(1. College of Agriculture， Guangxi University， Nanning 530004， China； 2. School of Computers， Guangdong University of Technology， Guangzhou 510006， China； 3. Guilin University of Technology： a. Guangxi Key Laboratory of Hidden Metallic Ore Deposits Exploration， b. College of Earth Sciences， Guilin 541004， China)，J.PlantResour. &Environ.， 2017， 26(2)： 1-9

TakingDicranopterispedata(Houtt.) Nakaike grown on tin polymetallic deposits (heavy metal stress area) in Dachang of Guangxi and on the periphery of mining area without mineralization or pollution (the control area) as experimental materials， the transcriptome high-throughput sequencing was conducted on the leaves ofD.pedata， and assembled unigenes were descripted by NCBI non-redundant protein sequence database (Nr)， NCBI non-redundant nucleotide sequence database (Nt)， KEGG Orthology database (KO)， Swiss-Prot database (Swiss-Prot)， Protein family database (Pfam)， Gene Ontology database (GO)， and euKaryotic Orthology Groups database (KOG). Meanwhile， differentially expressed unigenes between leaves ofD.pedatain heavy metal stress area and the control area were analyzed. The results show that 19.56 Gb clean data are acquired by sequencing， in which, leaves ofD.pedatain heavy metal stress area and the control area include 10.14 and 9.42 Gb clean data， respectively. 120 097 unigenes from 250 582 assembled unigenes are descripted, with 47.93% of total number of unigenes. Compared with the control area， there are 208 and 620 up- and down-regulated differentially expressed unigenes in leaves ofD.pedatain heavy metal stress area， respectively. In which， 120 up-regulated differentially expressed unigenes are descripted in metabolism process， accounting for 57.69% of all up-regulated differentially expressed unigenes， and 285 down-regulated differentially expressed unigenes are descripted in catalytic activity， accounting for 45.97% of all down-regulated differentially expressed unigenes. Fifteen unigenes in leaves ofD.pedatain heavy metal stress area are associated with heavy metal transport and tolerance， in which relative expression of c44988_g1 and c84121_g1 are extremely significantly and significantly higher than those in the control area, respectively. It is suggested that genes inD.pedatain response to natural metal mineralization or heavy metal pollution of mine could be used for biogeochemical prospecting and soil heavy metal pollution detection.

Dicranopterispedata(Houtt.) Nakaike； RNA-Seq； heavy metal stress； differentially expressed genes

2016-12-06

國家自然科學(xué)基金資助項目(41363003； 40972220)

鄒承武(1983—)，男，廣西賀州人，博士，主要從事植物分子生物學(xué)方面的研究。

①通信作者E-mail： gldysca@126.com

Q946-33； Q948.119

A

1674-7895(2017)02-0001-09

10.3969/j.issn.1674-7895.2017.02.01