基于轉錄組測序對紫花苜蓿鉛脅迫相關基因的富集分析

郭 強，王英哲，王 昆，徐 博

（1. 吉林農業(yè)大學動物科學技術學院，吉林 長春 130118；2. 吉林省農業(yè)科學院，吉林 長春 130118；3. 黑龍江省二九一農場畜牧獸醫(yī)科，黑龍江 集賢 155923）

金屬污染物如鉛(Pb)、鎘(Cd)和汞(Hg)是非必需元素，可以通過人類活動釋放到環(huán)境中，對動物和植物有毒害作用并且可以通過食物鏈威脅人類健康[1-2]。鉛是一種劇毒的重金屬，能夠被植物吸收，抑制根和芽的生長，導致葉片萎黃，影響其他生理過程，如有絲分裂、蒸騰和DNA合成[3-4]等。因此，恢復土壤環(huán)境狀態(tài)是當務之急，土壤清洗、土壤淋溶和植物修復等技術已被用于修復重金屬污染的土壤，其中植物修復被認為是一種高效、經濟、環(huán)保的清潔方法[5]。

了解基因的表達和調控在植物應對重金屬毒性反應以及植物修復涉及的遺傳和分子機制中顯得尤為重要。植物重金屬脅迫相關研究已經確定了一組基因家族和相關蛋白質參與金屬運輸，包括鐵響應轉運蛋白(IRT)、陽離子擴散促進劑(CDFs)、P型ATP酶、ATP結合盒(ABC)轉運蛋白、天然抗性相關巨噬細胞蛋白(NRAMP)和Fbox家族[6]。參與重金屬反應的基因在不同植物中的研究較廣泛[7-10]；然而，基因組網絡的基因表達對重金屬脅迫的反應仍然未知。最近的研究使用新一代測序技術鑒定了許多microRNA，確認為重金屬脅迫敏感的相關基因[11-13]。這些發(fā)現(xiàn)表明，植物中重金屬誘導的反應是復雜的，并且涉及多種生理過程和代謝途徑。

紫花苜蓿(Medicago sativa)是豆科植物，具有高度分枝的根系、高生物量，擁有對多種環(huán)境的廣泛適應性。這些特點使紫花苜蓿成為一種理想的土壤改良植物。了解紫花苜蓿的分子基礎對重金屬的反應將有助于其在植物修復中的應用。

1 材料和方法

1.1 試驗材料和處理條件

試驗采取先發(fā)芽后移栽的方式進行苗期培養(yǎng)。選擇20～30粒種子(隴東苜蓿，由吉林省農業(yè)科學院提供)放到鋪有3層濾紙的培養(yǎng)皿中，適量澆水后放入HPG-400HX人工氣候箱中培養(yǎng)，培養(yǎng)條件為(光照12 000 lx，14 h，溫度25 ℃，濕度75%；無光照，10 h，溫度20 ℃，濕度75%)。待種子發(fā)芽高度為1 cm時進行移栽，移栽到水培箱中(水培箱由箱體加浮板組成，浮板規(guī)格為28 cm × 18 cm，箱體規(guī)格為32 cm × 22 cm × 17 cm。營養(yǎng)液為Hoagland營養(yǎng)液[14]，每箱用量為3 L稀釋液，每兩箱栽種一份材料，每份材料栽種12棵苗，即每箱為6棵。將栽種好的材料分為兩組，一組為對照組，一組為處理組，培養(yǎng)30 d。30 d后進行脅迫處理，處理組以200 mg·L-1的硝酸鉛進行脅迫，用量為每3 L Hoagland營養(yǎng)稀釋液澆灌200 mL硝酸鉛溶液，處理96 h (記為Pb96)。對照組用蒸餾水處理，用量為每3 L Hoagland營養(yǎng)稀釋液澆灌200 mL蒸餾水，處理96 h (記為C96)。處理后取植株根部(預試驗得知根部較地上部效果明顯) 0.15 g進行試驗指標測定。取下的根經過液氮處理后保存于-80 ℃?zhèn)溆茫魈幚磉M行3次重復。本研究選取水培方法參考丁振中等[15]、程貝等[16]的研究確定。

1.2 RNA提取以及轉錄組測序

取Pb處理和對照植物根部的總RNA，使用RNAprep Pure Plant(Bio-Teke，China)試劑盒提取Pb96和C96。使用Qubit 3.0 (Thermo Scientific，USA)和Agilent 2100生物分析儀(Agilent Technologies，USA)評估RNA的濃度和質量。

使用NEBNext1 Ultra TM Directional RNA Library Prep Kit (NEB，UK)試劑盒根據說明制備cDNA文庫。分別使用Qubit2.0和Agilent 2100對文庫進行檢測，使用q-PCR方法對文庫的有效濃度進行準確定量，以保證文庫質量，準備好的cDNA文庫在Hiseq 4000(Illumina)上進行測序。

1.3 測序數(shù)據的過濾、組裝

在Q20條件下使用FastQC工具(http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc/)去除其中的接頭序列及低質量序列，獲得高質量的去冗余序列 (Clean Reads)。采用 Trinity(版本2.2.0，默認參數(shù))軟件進行序列組裝，獲得紫花苜蓿的功能基因(Unigenes)。

1.4 差異表達基因(differentially expressed genes，DEGs)的功能注釋與富集分析

使用FPKM方法計算基因表達水平并標準化。使用Student's t-test檢驗不同文庫之間的差異，P值根據Benjamini和Yekutieli[17]的方法進行調整。最后，以 FDR ＜ 0.005 且|log2(fold-change)| ≥ 1 為條件，篩選出Pb96和C96文庫中的DEGs。

利用BLAST程序(E-value = 10-5)將DEGs序列分別與COG、GO、KEGG(數(shù)據庫進行比對，獲得DEGs的注釋信息。

1.5 q-PCR驗證

PCR反應體系包括：SYBR premix Ex TaqTM(2X)10 μL，上游引物和下游引物各1 μL，cDNA 模板2 μL，RNase Free dH2O 6 μL。PCR 程序為：95 ℃ 預變性10 min；95 ℃ 變性 15 s，60 ℃ 退火 20 s，72 ℃ 延伸40 s；40個循環(huán)；80 ℃ 1 s收集熒光，讀板；之后72 ℃延伸10 min，最后從65 ℃開始，以每步0.5 ℃的速度升高至95 ℃，每個溫度保持1 s，繪制溶解曲線。每個樣品3個生物學重復。選用和ac-tin(CpActin)作為內參基因。

2 結果與分析

2.1 轉錄組測序與組裝

提取的RNA總濃度在324.8 ng·μL-1以上，表明RNA質量合格。Pb脅迫處理96 h與對照組的根部進行轉錄組測序后共獲得100 007條Unigene序列，所有序列N50值為1 097，在所有的Unigene中，長度在200～300 nt的最多，占Unigene序列總數(shù)的37.21%，長度超過2 000 nt的序列達到6 297條，占Unigene序列總數(shù)的6.30%(表1)。研究結果表明隨著Unigene序列長度的增加，Unigene序列呈現(xiàn)減少的趨勢，表明測序結果良好。

2.2 差異表達基因的鑒定

以FDR ≤ 0.01和|log2Ratio| ≥ 1作為DEG篩選條件，Pb96處理組與C96對照組相比，96 h Pb處理后共得到2 242個DEGs，占總Unigene序列的2.24%。其中，1 321個DEGs上調表達，921個DEGs下調表達。

表1 功能基因長度以及數(shù)量Table 1 Functional gene length and quantity

2.3 差異表達基因的功能注釋

通過BLAST(Evalue = 10-5)比對分析，2 242個DEGs序列GO數(shù)據庫中注釋到的DEGs數(shù)目最多，達到2 162個，占DEGs總數(shù)的96.4%；最少的為KEGG數(shù)據庫，注釋到的DEGs數(shù)目為448個，占DEGs總數(shù)的19.9%(表2)。

表2 差異表達基因的功能注釋結果Table 2 Functional annotation results of differentially expressed genes

2.4 差異表達基因的富集分析

COG數(shù)據庫的構建源自于細菌、藻類和真核生物的系統(tǒng)進化關系，利用該數(shù)據庫可以對基因產物進行直系同源分類。被注釋到COG數(shù)據庫的696條DEGs覆蓋所有 21個類別，其中“碳水化合物的運輸和新陳代謝”(Carbohydrate transport and metabolism)數(shù)量最多，為70個(10.06%)；具有“翻譯、核糖體結構和生物發(fā)生”(Translation，ribosomal structure，and biogenesis)、“次生代謝產物的生物合成、轉運和分解代謝”(Secondary metabolites biosynthesis，transport，and catabolism)和“僅通用功能預測”(General function prediction only)功能的DEGs數(shù)量較高，分別為 63 個(9.05%)、64個(9.20%)和68個(9.77%)(圖1)。可以看出，測序獲得的轉錄本參與了多數(shù)代謝過程，這些測序數(shù)據將是研究紫花苜蓿功能基因的可利用資源。

KEGG數(shù)據庫可以用來預測DEGs在細胞中參與的代謝途徑和功能，了解基因參與的生物學過程。將DEGs進行通路分析，發(fā)現(xiàn)共有448個DEGs參與到97個通路當中(部分展示，圖2)。其中DEGs參與最多的是“核糖體”(Ribosome)通路ko03010，共有48個DEGs參與；其次是“光合作用”(Photosynthesis)通路ko00195、“碳代謝”(Carbon metabolism)通路ko01200、“淀粉和蔗糖代謝”(Starch and sucrose metabolism)通路ko00500等。由此可知，植物在鉛脅迫狀態(tài)下，光合作用發(fā)生改變，同時體內的淀粉蔗糖代謝以及碳代謝均有明顯變化。

GO數(shù)據庫是一個結構化的標準生物學注釋系統(tǒng)，其信息適用于各物種。GO分類體系分為3 個大類，包括生物學過程(Biological process)、細胞元件(Cellular component)和分子功能(Molecular function)，以及57 個亞類。本研究中共有 2 242個DEGs被注釋到GO數(shù)據庫，分布于3個大類的24個亞類中，其中生物學過程中有9個亞類，細胞元件中有6個亞類，分子功能中有9個亞類(圖3)。被富集到物質代謝 (Metabolic process)245個 (GO:0008152)、蛋白結合(binding)214個(GO:0005488)、催化活性(Catalytic activity)202個(GO:0003824)、細胞過程 (Cellular process)162個(GO:0009987)、細胞(Cell)122個(GO:0005623)5個亞類的DEGs數(shù)量明顯較其他亞類多。說明植物在受到鉛脅迫時，植物體內蛋白含量和酶活性有明顯的變化。

通過注釋得到主要調節(jié)耐鉛機制的差異表達基因如表3所列，其中下調表達有12個，上調表達有10個。選擇其中差異最為顯著的c53338.graph_c0、c40437.graph_c1、c60665.graph_c0、c44469.graph_c0、c56258.graph_c0、c54935.graph_c1和c46893.graph_c0為候選基因，有待進一步分析。

2.5 qRT-PCR驗證結果

選取參與調節(jié)耐鉛性的10個差異表達基因，進行qRT-PCR驗證。結果表明，C96表達量均小于Pb96，這與高通量測序結果表達趨勢一致(圖4)。

圖1 注釋到COG數(shù)據庫的基因數(shù)目Figure 1 Notes the number of genes in the COG database

圖2 注釋到KEGG數(shù)據庫的基因數(shù)目Figure 2 Notes the number of genes in the KEGG database

3 討論與結論

鉛是世界范圍內廣泛存在的土壤污染物，對食物安全與人類健康構成巨大威脅[18-19]。不同的物理、化學和生物方法都有被用來解決環(huán)境中的重金屬污染問題，其中，植物修復被認為是一種高效、經濟、環(huán)保的方法[20-21]。但是，這些基因參與了重金屬毒性和分子的保護，這種保護的潛在機制在植物中仍未定義。

本研究通過BLAST軟件與GO、COG、KEGG數(shù)據庫對比最終獲得100 007條Unigene序列，鑒定出差異表達基因2 242個。在GO數(shù)據庫中注釋到2 162條，在COG數(shù)據庫中注釋到696條，在KEGG數(shù)據庫中注釋到448條。其余未被注釋的序列可能是紫花苜蓿所特有的基因序列，仍有待于進一步的分析。在GO注釋中的物質代謝、催化活性2個通路在黃海燕[22]的研究中同樣被注釋到，說明紫花苜蓿在鉛脅迫條件下，酶活性起到相對較高的作用。KEGG注釋中的“碳代謝”、“淀粉和蔗糖代謝”在宋麗莉[23]的紫花苜蓿抗鉛研究中也同樣被注釋到，說明淀粉蔗糖的代謝與合成能夠調節(jié)紫花苜蓿的耐鉛性。而在COG數(shù)據庫中注釋到的“翻譯、核糖體結構和生物發(fā)生”、“次生代謝產物的生物合成、轉運和分解代謝”功能的DEGs數(shù)量較高，這與齊曉[24]的研究結果相似，表明紫花苜蓿在抗鉛反應中主要調節(jié)作用在于物質的合成與分解。

圖3 注釋到GO數(shù)據庫的基因數(shù)目Figure 3 Notes the number of genes in the GO database

圖4 實時熒光定量PCRFigure 4 Real-time PCR

王少峽等[25]研究發(fā)現(xiàn)，DREB轉錄因子由逆環(huán)境脅迫誘導產生后，可激活其他多達 12個依賴DRE順式作用元件的抗逆功能基因，引起脯氨酸及蔗糖含量提高，從而增強植株對多種逆境的抵抗性。而在本研究中同樣注釋到“淀粉和蔗糖代謝”通路，表明差異表達基因促進了淀粉蔗糖代謝過程，也充分證明了淀粉蔗糖代謝對紫花苜蓿的抗鉛性至關重要。

次生代謝產物是植物長期進化中適應環(huán)境的產物，其公認的生態(tài)學功能主要是抗病、抗蟲、抗環(huán)境脅迫[26]等。而在COG數(shù)據庫中同樣注釋到“次生代謝產物的生物合成、轉運和分解代謝”，且注釋到的差異基因數(shù)量為64個，位居第二，表明次生代謝產物對于紫花苜蓿適應逆境環(huán)境具有較為重要的作用。

n e n t n i o p o a l r t e r，a t e r c h a n n e l a c t i v i n t e g r a l c o m 性活i t y，酶接合泡液連結o A合糖央-C結多a l s u b u n i t s中酸丁離子，分運運，轉o t p e p t i d e)屬成i c l a r g e r i b o s o m程e t a b o l i s m，金轉體膜s t a s i s，w分過跨、整謝膜p h o s p h o r y l a t i o n，成r t e r a c t i v i t y，的跨)鍵i c l a r g e r i b o s o m子e o 體膜代g l y o x y l a t e c y c l e，鹽，整肽，e s，c y t o p l a s m，的G e n e f u n c t i o n 酸活白s t r a t e-s p e c i f i c t r a n s m e m b r a n e t r a n s p o d i n g 膜離分r a n e t r a n s p o化是r a n e t r a n s p o r t，a m i n o a c i d t r a n s m e m b 亞動a t e r h o m不中a c t s o n c a r b o n-n i t r o g e n (b u t n亞e s，c y t o p l a s m 丁活，基p o l y s a c c h a r i d e b i n (但蛋，酸體磷成體道e t a b o l i s m，程環(huán)-氮質運整通過性糖基核循轉白謝的膜n e n t，s u b碳能膜酸代p o u n d s，c a r b o h y d r a t e m d i n g水體i t y的蛋功大f r i b o s o m ，質跨分c t u r a l c o m 醛p o大p o n e n t s o f r i b o s o m，態(tài)物因性i t y，i n乙胞結A T 糖性穩(wěn)性合l c o m合脒，P b i n d i n g，p r o t e i n活核基，異水，a c e t i c a c i d m化g e n e s 性T P 特胞程水成線活，y d r o l a s e a c t i v r a l c o m b r a n e，l y c o s y物構s t r u 蛋質細過l i g a s e a c t i v碳于，分p o n e n t s o謝e t a l i o n b，A白物f O 酶表r e s s e d ，胞o s p h a t e a l d o l a s e a c t i v i t y N u t r i e n t s t o r a g e a c t i v i t i e s，運b l y 縮成b r a n e i n t e g用的b l y，o A 底a c t i v i t y，t r a n s m e m b r a n e t r a n s p o r t 結膜p l a n t-t y p e t o n o p l a s t，c e l l w，r t，c e n t r a l v a c u o l e性代醋l i g a s e a c t i v i t y，作-g e m 活醛分泡膜酸合體e m酸液構化性磷解r o c e s s，h a n d m糖型p o n e n t s o f t h e m n d s i n l i n e a r r u t h e n i u m成跨列酶酸，基激基因- 二結l l y e x p活r a l c o m核體e m整p o n e n t s，物性e A b u t y r a t e-C d r o l y s i s o b r a n e基異i f f e r e n t i a 酸轉酶氨的，u n i t a s s e m的植活O-糖h y o f t h e m達糖氨，體裝膜，酶解表果水水H y d r o l a s e a c t i v i t y, f r u c t o s e-d i p h/蘇運糖動組b r a n e ，分接o e n z y m水，，核r e o n i n e k i n a s e a c t i v i t y，酸N u c l e o s o m e a s s e m活基膜成b r a n e c o m連o l i s m，，p o s i t i o n性3 差程氨表L i s t o f d 膜，c o l y s i s p存裝亞泡A e t a b性跨活過e m 絲合儲成e m組大液p e t o n o p l a s t，m組b r a n e t r a n s p o酶活酸酶解結酵e g l y 質物體體型t-t y 的-輔酶r a t e m解糖T h轉b o 白基蛋P r o t e i n s e r i n e/t h i n t e g r a l c o m氨A m i n o a c i d t r a n s m e m b A A b i n d i n g，s t r u c t u養(yǎng)小糖物膜酸t(yī) y 解e r a l l c o m b l e 3 o f t h e m R N R N 營核核R i b o s o m e l a r g e s u b s u b u n i t s植P l a n質P l a s m a m e m t r a n s m乙A c e t i c a c i d - C b u 水H y d r o l a s e a c t i v i t y，o v T a比g 2 r a t i o 1 8值4 4 2 0 7 4 6 0 1 4 5 3 4 6 9 4 1 8 6 5 8 8 4 5 6 5 8 8 4 5 3 1 8 0 2 8 3 1 8 0 2 8 0 9 5 6 2 4 9 9 2 4 5 9 6 2 6 4 6 5 6 9 5 6 6 1 l o-1.0 6 3-1.6 3 3-3.6 3 1-3.6 2 5-3.6 2 5-3.5 5 8 5 8 2 8 3-3.5-3.5 5 8-3.4 5 1-2.9 5 2-2.6 9 3-1.9 6 1 M P K P b 9 6_F 6.1 5 1 0 6.6 5 0.1 3 0.2 5 0.0 8 0.1 4 0.0 4 0.0 5 0.2 3 0.0 7 0.7 1 0.7 9 K M F P C 9 6_1 9 8.6 9 1 8.9 7 2.2 2 6 1.9 1 2.9 9 0.8 7 1.2 6 3.3 2 0.7 4 4.5 3.0 7 e n e名U n i g i g e n e I D r a p h_c 0.g r a p h_c 1.g r a p h_c 0.g r a p h_c 0.g r a p h_c 0.g r a p h_c 0.g r a p h_c 0.g r a p h_c 0.g r a p h_c 0.g r a p h_c 0.g r a p h_c 0.g r a p h_c 0.g U n 4 7 8 2 c 5 0 4 3 7 c 4 3 9 9 2 c 4 0 6 6 5 c 6 0 9 0 7 c 4 3 8 1 1 c 3 1 6 9 4 c 3 0 9 8 4 c 3 5 5 0 0 c 3 4 6 3 8 c 4 4 4 6 9 c 4 0 2 9 5 c 4

?

在上下調基因中(表3)基因的主要功能是調節(jié)質膜、酶活以及跨膜轉運和離子結合，其中在酶活調節(jié)過程中，主要調節(jié)輔酶、水解酶以及氧化還原酶。在調節(jié)跨膜轉運過程中主要調節(jié)水的轉運和離子轉運。同時在基因調解過程中還有其他部分細胞器的參與。綜上所述，通過數(shù)據庫的比對發(fā)現(xiàn)，物質的合成與分解、淀粉蔗糖代謝以及次生代謝產物3大過程均對于植物耐鉛性起到調節(jié)作用。

本研究中將差異表達基因與3大數(shù)據庫對比發(fā)現(xiàn)許多與耐鉛相關的通路，在紫花苜蓿中獲得了大量響應鉛脅迫的DEGs，并在轉錄組水平鑒定了參與調控紫花苜蓿鉛脅迫反應的基因，從而有利于從整體水平加深了解紫花苜蓿耐鉛反應表達特性，為進一步研究這些響應鉛脅迫的基因的功能提供參考，也可為紫花苜蓿耐鉛功能研究提供新的基因資源。