硒處理下冰菜的轉錄組響應及相關基因功能分析

姚運法，練冬梅，林碧珍，賴正鋒，洪建基

(福建省農業科學院亞熱帶農業研究所，福建 漳州 363005)

【研究意義】冰葉日中花(MesembryanthemumcrystallinumLinn.)俗稱冰菜，屬于番杏科日中花屬(MesembryanthemumL.)的一年生或二年生草本植物，其莖匍匐、葉肉質，具有較高食用價值[1]。原產非洲南部，世界各地均有引種栽培，近年來引種中國[2]。冰菜在環境適應性方面，具有耐干旱、鹽堿、喜陽光、怕寒及水澇等特點，其適合生長的溫度是5～25 ℃[3]。冰菜作為一種新型特菜資源，研究冰菜外源硒對冰菜生長和生理變化機制具有重要意義。【前人研究進展】硒是一種重要的微量元素，適量硒能清除植物體過量自由基[4]，促進根系生長發育[5]，有利于増強植物體新陳代謝，提高生命力，抵御逆境傷害，延緩組織成熟衰老等[6-8]；高濃度硒會促進植物過氧化作用，易造成植物體損傷[9]。研究表明，硒能改善多種作物的可食用部分品質[4]。蔬菜食用部位硒的含量濃度在0.0008～5.37 mg/kg范圍內，均值為0.067 mg/kg[10]，其中葉菜中硒的平均濃度低于推薦值0.025 mg/kg[11]。蔬菜中的硒含量隨著土壤中硒含量的增加而増加，對外源硒的響應能力因植物種類而異，薛瑞玲等[12]發現，低濃度硒鹽能提高小白菜的抗氧化作用，促進葉綠素合成和葉片生長，高濃度時則相反。因此，適宜濃度的硒對植物的生長發育有益。【本研究切入點】近年來，在中國的冰菜設施栽培越來越廣泛，開展冰菜富硒栽培模式具有重要經濟價值。【擬解決的關鍵問題】本試驗對冰菜以不同濃度(0、10和100 mmol/L)外源硒處理，利用Illumina Hi-seq 2500高通量測序技術，分析不同濃度硒處理下冰菜響應機制及其植物激素信號轉導途徑相關差異基因表達情況，篩選冰菜硒響應的相關基因及其表達機制，在轉錄組水平闡述冰菜對硒處理的響應機制，挖掘發現富硒的冰菜資源，為富硒冰菜育種積累理論基礎。

1 材料與方法

1.1 供試材料

試驗材料為冰菜，種植于福建省農業科學院亞熱帶農業研究所試驗農場。2019年12月11日育苗，2020年1月10日定植。2020年2月15日，以亞硒酸鈉為外源硒，不同濃度(M0: 0 mmol/L; M1:10 mmol/L；M2: 100 mmol/L)亞硒酸鈉處理冰菜。2020年2月22日取嫩莖葉，每處理設置3個重復，共9個樣品，置于-80 ℃超低溫冰箱備用。

1.2 試驗方法

1.2.1 轉錄組測序 利用Trizol(Invitroge)試劑分別提取硒處理(M1和M2)和對照(M0)冰菜嫩莖葉總RNA，采用Nanodrop(ThermoFisher)、Qubit 2.0、Aglient 2100(Bioanalyzer)技術檢測RNA樣品的純度、濃度和完整性等，構建cDNA文庫，再分別使用QubiM2.0、Agilent 2100和Q-PCR對文庫的質量進行檢測。合格后，基于SBS技術，使用Illumina Hiseq高通量測序平臺對cDNA文庫進行測序，產出高質量Raw Data，采用截除Reads中的測序接頭以及引物序列方式，過濾低質量值數據，獲得高質量Clean Data。利用Trinity軟件對高質量的Clean Data進行序列組裝。

1.2.2 測序數據組裝與分析 將冰菜轉錄組測序Clean Data與組裝得到的Transcript或Unigene庫進行序列比對，比對到Transcript或Unigene的Reads稱為Mapped Reads，Mapped Reads將用于后續的分析。

1.2.3 Unigene表達量估計與差異表達分析 采用Bowtie將測序得到的Reads與Unigene庫進行比對，根據比對結果，結合RSEM進行表達量水平估計。利用FPKM值(Fragments Per Kilobase of transcript per Million mapped reads)表示對應Unigene的表達豐度。使用EBSeq進行差異表達分析。采用Benjamini-Hochberg方法對原假設試驗得到的顯著性P值進行校正，校正后P值，即偽發現率(False discovery rate)小于0.01且差異倍數(Fold Change,FC)大于等于2作為篩選標準。其中，FC表示兩樣品(組)間FPKM的比值。

1.2.4 DEGs篩選和韋恩圖統計 根據對照與處理樣品之間表達水平的相對高低，DEGs可分為下調基因(Down gene)和上調基因(Up gene)。然后對兩組 (M0 與 M1; M0 與 M2) 數據DEGs進行維恩圖分析。

1.2.5 Unigene功能注釋與DEGs富集分析 使用Blast軟件將Unigene序列與COG、GO、KEGG、KOG、Pfam、Swiss-Prot、eggNOG、NR等8大數據庫比對，使用KOBAS 2.0得到Unigene在KEGG中的KEGG Orthology結果，預測完Unigene的氨基酸序列之后使用HMMER軟件與Pfam數據庫比對，獲得Unigene的注釋信息。系統分析基因產物代謝途徑及功能，并將DEGs比對到KEGG數據庫，得到DEGs的代謝途徑。

1.2.6 硒處理對冰菜植物激素信號轉導途徑分析 對硒處理條件下冰菜植物激素信號轉導途徑DEGs(關鍵酶基因)進行轉錄組數據庫檢索，分析植物激素信號轉導途徑關鍵基因表達情況，研究不同硒處理濃度對植物激素信號轉導途徑中間產物影響。

1.2.7 植物激素信號轉導途徑關鍵差異基因驗證 表1顯示，取1 μg冰菜葉片總RNA，利用反轉錄試劑盒反轉錄成cDNA，利用Primer Premier 6軟件，設計并合成引物，使用qRT-PCR分別檢測M0、M1、M2處理與植物信號轉導代謝途徑關鍵差異表達基因，設置3個重復；統計9個基因(c38220為內參)在待測樣品中的Ct值，計算相對表達量。

2 結果與分析

2.1 數據組裝及分析

表2顯示，經高通量測序和質量控制分析，共獲得23.63 Gb Clean data，其中M0 Clean data為8.47 Gb，M1 Clean data為7.61 Gb，M2 Clean data為7.55 Gb，堿基百分比均達到92.0%以上。本次測序數據質量好，利于后續分析。由表3顯示，轉錄本序列(Transcipt)組裝出223 477條轉錄本序列，平均長度為1843.9 bp，N50長2916 bp；由單基因序列(Unigene)組裝出76 806條，平長度為759.0 bp，N50長1513 bp。

表1 實時熒光定量PCR引物

表2 有效數據評估統計

表3 組裝結果統計分析

表4 單基因序列注釋統計表

表5 差異基因表達情況

2.2 Unigene功能注釋

如表4所示，在76 806條單基因序列中，共獲得28 172條Unigene的注釋結果，占單基因序列總數36.68%。其中與COG數據庫比對，獲得9096條基因注釋，占總注釋基因的32.29%；與KOG數據庫比對，獲得16 137條基因注釋，占總注釋基因的57.28%；與GO數據庫比對，獲得10 493條基因注釋，占總注釋基因的37.25%；與KEGG數據庫比對，獲得15 636條基因注釋，占總注釋基因的55.50%；與Pfam數據庫比對，獲得18 324條基因注釋，占總注釋基因的65.04%；與Swiss-Prot數據庫比對，獲得17371條基因注釋，占總注釋基因的61.66%；與eggNOG數據庫比對，獲得25 240條基因注釋，占總注釋基因的89.59%；與NR數據庫比對，獲得26 405條基因注釋，占總注釋基因的93.73%。

2.3 DEGs的篩選與功能注釋

2.3.1 DEGs篩選 如表5所示，通過DEGs的篩選，獲得冰菜硒處理DEGs 570個(M0 與 M1)和580個(M0 與 M2)，其中M0 與 M1組表達量上調的基因493個，下調基因77個；M0 與 M2組表達量上調的基因293個，下調基因287個。通過對兩組(M0 與 M1，M0 與 M2)DEGs維恩圖分析，硒處理(M0 與 M1)非共表達基因468個，其中上調表達基因447個，下調表達基因21個；硒處理(M0 與 M2)非共表達基因478個，其中上調表達基因247個，下調表達基因231個；硒處理共表達基因102個，其中上調表達基因為46個，下調表達基因56個(圖1)。

表6顯示，將M0與M1和M0與M2兩組DEGs單基因序列分別注釋到COG、GO、KEGG、KOG、Pfam、Swiss-Prot、eggNOG、NR8個數據庫，分別獲得426和428個基因，其中COG數據庫分別208和144個，GO數據庫分別241和239個，KEGG數據庫分別237和138個，KOG數據庫分別309和168個，Pfam數據庫分別366和330個，Swiss-Prot數據庫分別289和325個，eggNOG數據庫分別404和381個；NR數據庫分別367個和422個，其中M0與M1組在eggNOG數據庫注釋比最高，超過94.8%；M0與M2組在NR數據庫注釋比最高，超過98.6%。

表6 差異表達基因注釋數量統計

2.3.2 GO功能注釋 如表7所示，對不同硒濃度處理M1、M2與M0(CK)DEGs進行GO分類統計，兩組DEGs均歸到43個功能小類。生物學過程中DEGs主要集中在代謝過程、細胞過程、單一生物過程、生物調節、刺激應答和細胞成分或生物合成等6個功能小類占比高，兩組DEGs差值來看，單一生物過程和刺激應答生物學過程，隨著硒處理濃度增加，DEGs分別增加29和22個；細胞過程和細胞成分或生物合成生物過程，隨著硒處理濃度增加，DEGs分別減少45和24個。細胞組分過程中DEGs主要在細胞器、復雜大分子、細胞成分、細胞、細胞器成分、膜成分、膜結構7個功能小類占比高；兩組DEGs差值來看：隨著硒處理程度增加，細胞器、復雜大分子、細胞成分、細胞、細胞器成分、細胞組分過程DEGs減少，且DEGs均在50個以上；膜成分、膜結構和細胞器成分等過程比值較高，分別差42、35和18個DEGs。分子功能過程中DEGs在催化活性、結合活性、轉運活性和分子結合活性4個功能小類占比最高。兩組DEGs比值差異來看：其中催化活性過程與硒處理濃度呈正相關，DEGs個數達69個；結構分子活性過程與硒處理濃度呈負相關，DEGs數目達59個。

表7 硒處理差異表達基因GO功能注釋比較

續表7Continued table 7

2.3.3 KOG功能注釋 將M0與M1和M0與M2兩組DEGs注釋到KOG數據庫，對注釋結果進行直系同源分類，均獲得22個功能分類，其中未知功能(S)、無機離子轉運與代謝(P)、胞質運輸、分泌和囊泡運輸(U)、碳水化合物轉運與代謝(G)、能源的產生和轉化(C)、次生代謝產物的合成、轉運和代謝(Q)、脂質轉運與代謝(I)、翻譯、核糖體結構和生物合成(J)、氨基酸轉運與代謝(E)、轉錄(K)、無機離子轉運與代謝(P)共11個功能分類DEGs差異較多。

兩組KOG功能注釋差值分析，一般功能預測(R)、次生代謝產物的合成、轉運和代謝(Q)、碳水化合物轉運與代謝(G)、氨基酸轉運與代謝(F)功能分類DEGs數分別為34、14、13、12個；翻譯、核糖體結構和生物合成(J)、能源的產生和轉化(C)、無機離子轉運與代謝(P)、細胞骨架(Z)、翻譯后修飾、蛋白質轉換、伴侶(O)功能分類DEGs數分別為70、27、13、10、10個(表8)。

表8 硒處理差異表達基因KOG功能注釋比較

2.3.4 KEGG功能注釋 表9顯示，將兩組DEGs與KEGG數據庫比對，M0與M1和M0與M2分別有431和351條基因得到注釋，分別富集在104和101條代謝通路，其中M0與M1組中，硒處理對核糖體(68個)、碳代謝(35個)、氧化磷酸化(30個)、氨基酸的生物合成(27個)、植物激素信號轉導(26個)、苯丙烷類生物合成(19個)、植物—病原互作(18個)、光合作用的固碳作用(17個)、內質網蛋白加工(16個)、嘌呤代謝(15個)、丙酮酸代謝(15個)等11個代謝途徑具有主要影響；M0與M2組中，硒處理對植物激素信號轉導(34個)不飽和脂肪酸的生物合成(30個)、苯丙烷類生物合成(30個)、脂肪酸代謝(18個)4個代謝途徑具有主要影響。M0與M1和M0與M2兩組DEGs差值分析表明，不飽和脂肪酸的生物合成、脂肪酸代謝、苯丙烷類生物合成、植物激素信號轉導代謝途徑其DEGs增加分別為27、12、11和8個；核糖體、碳代謝、氨基酸的生物合成、氧化磷酸化、植物—病原互作、內質網蛋白加工代謝途徑其DEGs分別減少66、32、25、21、9和9個。

表9 硒處理差異表達基因KEGG功能注釋比較

2.4 植物激素信號轉導途徑分析

生長素(Auxin)信號轉導途徑中，生長素吲哚乙酸蛋白(Auxin/indole-3-Acetic Acid，AUX/IAA)有2個差異拷貝，M1基因表達無顯著差異，M2處理均表現為下調表達；生長素上調小RNA(Small auxin up RNA,SAUR)共16個差異拷貝，M1和M2處理均具有上調和下調基因；生長素響應家族基因(Gretchen Hagen3，GH3)共有2個拷貝，M1和M2均表現下調(c32701)或上調(c32788)。細胞分裂素(Cytokinine，CK)信號轉導過程中，磷酸轉運蛋白(Arabidopsis histidine phosphotransfer proteins，AHP)(c31251、c36579)在M1和M2處理下均表現為上調表達，AHP(c39901)對M1處理表現下調，M2處理表現不敏感。A型擬南芥響應調節因子(Type-A Arabidopsis response regulators，A-ARR)對M1和M2處理顯著上調作用。赤霉素(Gibberellin，GA)信號轉導過程中，赤霉素不敏感矮稈蛋白1(Gibberellin insensitive Dwarf 1，GID1)(c39849)在M1和M2處理均表現為上調表達，轉錄因子(Transcription Factor，TF)(c37269)在M1和M2處理均表現為下調表達。脫落酸(Abscisic acid，ABA)信號轉導途徑中：蔗糖非發酵1相關蛋白激酶2(Sucrose nonfermenting 1-related protein kinase2，SnRK2)(c37312)在M1處理表現為下調表達，M2處理表現不顯著表達，轉導受體蛋白PYR/PYL/RCAR家族基因影響比較復雜，PYR/PYL/RCAR復合體(c26211)在M1處理表現為上調表達，M2處理表現為下調表達；PYR/PYL/RCAR復合體(c34313)在M1處理表達不顯著，M2處理顯著下調表達；PYR/PYL/RCAR(c28738)在M1和M2均表現上調表達；2C類蛋白磷酸酯酶(Protein phosphatase 2C,PP2C)有7個拷貝具有表達量差異表現，無規律；ABRE結合因子(ABRE binding factors，ABF)有2個差異表達拷貝，其中c34215在M1處理表現下調表達，M2處理表達差異不顯著，c32284在M1處理表達差異不顯著，M2處理表達顯著下調。乙烯(Ethylene，ET)信號轉導途徑中，ETR和ERF1/2各有1個拷貝，在M1和M2處理下均表現為顯著下調表達。茉莉酸(Jasmonic acid，JA)信號轉導途徑中，骨髓細胞瘤病蛋白2(Myeolcytomatosis proteins 2，MYC2)(c36305)和茉莉酮酸酯ZIM結構域蛋白(Jasmonate ZIM-domain，JAZ)(c28853)在M1處理下均表現顯著上調表達，M2處理下表達無顯著差異。水楊酸(Salicylic acid，SA)信號轉導途徑中，TGACG基序結合因子(TGACG motif-binding factor，TGA)和病程相關蛋白-1(Pathogenesis-related protein-1，PR-1)在M1處理條件下，表達量無顯著差異表達，M2處理下，TGA表達量下調，PR-1表達量上調，上調表達量增加12倍。另外，油菜素類固醇(Brassinosteroid，BRs)信號轉導途徑在M1和M2處理下，均無顯著差異表達基因(表10)。

表10 硒處理下冰菜植物激素信號轉導途徑關鍵基因差異表達

續表10Continued table 10

2.5 植物信號轉導途徑關鍵差異基因驗證分析

對表10中部分DEGs進行熒光實時定量驗證，共驗證AUX/IAA、SAUR、GH3、AHP、MYC2、JAZ、PR-1等7個。以冰菜Actin (c38220) 基因為內參基因，進行qRT-PCR驗證。由圖2可知，7個DEGs基因qRT-PCR分析得到的相對表達量與轉錄組表達譜分析趨勢完全一致，但表達的變化大小存在一定差異，說明基因表達譜的分析結果基本可靠，其中PR-1(c41466)基因分別M0和M2表達量差異極大，為冰菜水楊酸信號轉導途徑硒處理響應關鍵基因克隆提供研究基礎。

3 討 論

硒處理對植物激素信號轉導途徑中除油菜素類固醇外，其它7個代謝途徑均產生重要影響。生長素IAA通過細胞質膜上生長素輸入載體AUX1轉運至細胞內，其受體TIR1與AUX/IAA互作，促進AUX/IAA泛素化而降解，解除生長素對ARF活性的抑制作用[13-14]。隨著硒處理程度加深，冰菜葉片中AUX/IAA被泛素化降解，解除AUX對下游ARF活性的抑制作用，促進生長素信號轉導；另外，SAUR和GH3對硒處理的響應較為復雜，且不具有規律性。胞外細胞分裂素與質膜上受體蛋白(Cytokinin response1,CRE1)結合，使得CRE1發生自磷酸化，然后其磷酸基團被傳遞至細胞質中的AHP上，磷酸化的AHP進入到細胞核中，再將磷酸基團轉移至B型擬南芥應答調節子(B-Arabidopsis Response Regulators,B-ARR)，作為正調控轉錄因子，磷酸化的B-ARR激活下游基因A-ARR等表達，促進細胞分裂，A-ARR作為負調控因子，其累積反過來抑制B-ARR活性，從而抑制細胞分裂[15-16]。冰菜硒處理過程中，AHP對硒處理敏感，但對其濃度變化不敏感；A-ARR對硒處理顯著上調作用，但對濃度變化不敏感；硒處理M1與M2處理均能顯著性的提高A-ARR的表達量，顯著的抑制細胞分裂素的細胞分裂作用，抑制冰菜生長和發育作用。赤霉素信號轉導過程中，隨著硒處理濃度增加而顯著提高，GID1與DELLA蛋白結合，異構化DELLA蛋白空間構象，異構化后的DELLA蛋白在GID2作用下，進行泛素化修飾后被26S蛋白酶降解，其下游TF得到激活[17-18]，促進冰菜莖葉生長，但硒處理也對赤霉素信號轉導途徑下游TF(c37269)表達量顯著抑制作用。脫落酸信號轉導途徑中，低濃度硒處理(10 mmol/L)導致SnRK2表達量顯著降低，從而負調控ABA信號轉導下游轉錄因子ABF，抑制氣孔關閉，導致冰菜易脫水[19-20]；高濃度硒處理(100 mmol/L)SnRK2蛋白表達量顯著恢復。另外，硒處理對ABA信號轉導受體蛋白PYR/PYL/RCAR家族基因影響比較復雜，無顯著規律。乙烯信號轉導途徑中，乙烯信號分子傳到內質網乙烯受體蛋白(Ethylene receptor,ETR)上，ETR與下游乙烯抑制蛋白(Constitutive triple response 1,CTR1)互作，通過激酶磷酸化下游乙烯不敏感蛋白(Ethylene insensive 2,EIN2)，從而抑制EIN2下游的乙烯響應途徑；此外，EIN2還受到EIN3結合F盒蛋白(EIN3-binding F-box protein1/2,EBF1/2)介導的泛素化蛋白酶降解調控，進而抑制ERF等下游乙烯響應基因[21-22]。隨著硒處理濃度的提高，ETR和ERF1/2表達量降低，抑制EIN2和ERF下游的乙烯響應途徑，抑制冰菜成熟和衰老[23]。茉莉酸信號轉導途徑中，JA在JAR1酶催化下將異亮氨酸(Isoleucine,Ile)加到JA上，獲得活化性JA-Ile，活性JA-Ile大量積累，結合JA受體蛋白—冠毒素不敏感蛋白(Coronatine-insensitive protein 1,COI-1)，COI-1結合JAZ形成復合體，促進JAZ泛素化而降解，進而下游MYC2等轉錄因子釋放，激活JA下游響應基因表達，調控植物非生物逆境響應[24]。低濃度硒處理(10 mmol/L)促進JAZ和MYC2表達量上調，調節硒對冰菜的逆境脅迫；高濃度硒處理(100 mmol/L)對JAZ和MYC2表達量影響不顯著，表明高濃度硒處理對茉莉酸信號轉導途徑無顯著影響。水楊酸信號轉導途徑中，SA與細胞質中受體蛋白—非發病相關蛋白表達子(Nonexpresser of Pathogenesis-Related genes,NPR1)受體結合，以寡聚體形式的NPR1解離為單體形式，單體NPR1進入細胞核，與轉錄因子TGA互作，共同激活下游基因表達，激活植物抗非生物逆境能力[25-26]。本實驗中，TGA和PR-1在低硒處理(10 mmol/L)條件下，表達量無上下調表達，高硒處理(100 mmol/L)下，TGA表達量下調，PR-1表達量上調，其中PR-1上調表達量增加12倍，激活冰菜硒響應能力。經實時定量PCR分析，其中7個DEGs得到的相對表達量與轉錄組表達譜分析趨勢完全一致。

4 小 結

利用高通量轉錄組測序技術，前期項目組初步對冰菜鹽脅迫相關代謝分析和關鍵基因挖掘[27]，在此基礎上，繼續開展冰菜硒處理試驗，研究植物信號轉導途徑響應硒處理的分子機制，為下一步研究硒調控過程的轉錄因子及其關鍵功能基因、信號轉導途徑提供參考。以上鹽脅迫、硒處理逆境冰菜分子機制的研究，有助于冰菜抗逆品種的改良及選育，并指導富硒冰菜生產實踐上具有重要意義。