PEG和NaCI脅迫下毛竹萌發(fā)種子的MicroRNAS表達(dá)譜分析

王曉靜 王濤 楊凱 李潞濱

關(guān)鍵詞：毛竹；種子萌發(fā)；PEG；NaCI；miRNA

毛竹（PhyHostachys edulis （Garr.）H.deLehaie）屬于禾本科（Gramineae）竹亞科（Bambusoideae）剛竹屬（PhyHostachys Sieb.Et Zucc.），是一種多年生木本散生竹，廣泛分布在460 N以南地區(qū)，是我國(guó)種植面積最廣的竹種之一，具有重要的經(jīng)濟(jì)、生態(tài)和文化價(jià)值。毛竹開花后產(chǎn)生典型穎果類果實(shí)，薄而干的膜質(zhì)果皮與種皮緊密貼合，在種子生物學(xué)和苗木培育上毛竹穎果通常被稱為種子。毛竹種子萌發(fā)是實(shí)生苗造林的基礎(chǔ)，通過毛竹種子萌發(fā)進(jìn)行實(shí)生苗造林能降低運(yùn)輸成本，便于移栽，有利于提高造林整齊度、成活率和遺傳多樣性，延長(zhǎng)竹林高產(chǎn)時(shí)間，且實(shí)生苗竹林萌筍多、竹鞭生長(zhǎng)較快、竹材利用率更高，因此，開展毛竹種子萌發(fā)研究對(duì)于竹類資源的培育和應(yīng)用等具有重要意義。

隨著氣候變化的加劇，干旱和由土壤初級(jí)或次級(jí)鹽堿化造成的鹽脅迫已經(jīng)成為全球性問題，并對(duì)種子萌發(fā)、植物生長(zhǎng)等造成嚴(yán)重影響。干旱和鹽脅迫抑制了毛竹種子的萌發(fā)，進(jìn)而影響了毛竹實(shí)生苗的生長(zhǎng)，限制了通過毛竹種子萌發(fā)進(jìn)行實(shí)生苗培育和應(yīng)用。目前，研究者逐步開展了干旱或鹽脅迫對(duì)毛竹種子萌發(fā)率、種子活力、種子壽命和萌發(fā)后幼苗生長(zhǎng)等生理方面的研究，但相關(guān)的分子調(diào)控機(jī)制鮮有報(bào)道。

MicroRNA（miRNA）是一類長(zhǎng)度約為21～24nt的小分子非編碼RNA，能夠通過互補(bǔ)配對(duì)抑制或降解靶基因進(jìn)行轉(zhuǎn)錄后調(diào)控。miRNA廣泛參與了植物的生長(zhǎng)發(fā)育，在種子萌發(fā)和非生物脅迫響應(yīng)中具有重要的調(diào)控作用。在竹類植物中，研究者已經(jīng)在組學(xué)水平開展了關(guān)于毛竹葉片、花發(fā)育、雄蕊發(fā)育、莖稈快速生長(zhǎng)和竹筍發(fā)育等的miRNA研究，但毛竹種子中miRNA的表達(dá)及其在干旱和鹽脅迫下萌發(fā)時(shí)的調(diào)控仍是未知的，因此，本研究通過small RNA測(cè)序?qū)γ穹N子露白階段的miRNA進(jìn)行系統(tǒng)分析，并探究其在不同PEG或NaCI脅迫下萌發(fā)時(shí)的表達(dá)模式，挖掘毛竹種子萌發(fā)階段調(diào)控干旱或鹽脅迫抗性潛在的關(guān)鍵miRNA，以期為后續(xù)的研究提供基礎(chǔ)和參考。

1材料與方法

1.1試驗(yàn)材料

試驗(yàn)所用毛竹種子于2020年9月收集自廣西壯族自治區(qū)桂林市靈川縣，剝除外稃后挑選大小均一、顆粒飽滿、色澤明亮的完整種子備用（圖1a）。種子的千粒質(zhì)量為21.85g±0.04g，含水量為12.05%±0.35%。種子表面消毒后在培養(yǎng)皿中使用紙床法進(jìn)行萌發(fā)。聚乙二醇（PEG6000）是一種親水性大分子物質(zhì)，是模擬植物干旱脅迫的常用試劑之一，本研究使用PEG6000模擬干旱脅迫環(huán)境，同時(shí)使用氯化鈉（NaCI）模擬鹽脅迫。使用完全隨機(jī)實(shí)驗(yàn)設(shè)計(jì)，將種子分為5組（A～E）：A組為對(duì)照組，B組、C組模擬干旱脅迫，D組、E組模擬鹽脅迫。參考前期預(yù)實(shí)驗(yàn)的結(jié)果，分別在培養(yǎng)皿中加入H20（A）、10% PEG（B）、15%PEG（C）、50mmoI.L-1 NaCI （D）和100mmoI.L-1NaCI（E）．每組萌發(fā)200粒種子，并重復(fù)3次。每天觀察并記錄種子萌發(fā)狀態(tài)，在第4天選擇生長(zhǎng)狀態(tài)一致且均到達(dá)種皮破裂階段的種子（圖1b）。對(duì)于每個(gè)處理，在3個(gè)重復(fù)中各取5粒種子并將15粒種子混合取樣后用于建庫(kù)。取樣后使用液氮速凍并置于-80℃保存。

1.2建庫(kù)、測(cè)序與數(shù)據(jù)過濾

使用植物多糖多酚RNA提取試劑盒（DP441，TIANGEN）提取樣品總RNA。在建庫(kù)時(shí)首先根據(jù)small RNA的結(jié)構(gòu)特點(diǎn)將總RNA中small RNA加接頭并進(jìn)行純化，隨后進(jìn)行反轉(zhuǎn)錄、PCR擴(kuò)增和聚丙烯酰胺凝膠電泳，切膠回收目的條帶即得到small RNA的cDNA文庫(kù)。文庫(kù)有效濃度合格后進(jìn)行測(cè)序，使用lllumina高通量測(cè)序平臺(tái)NovaSeq6000進(jìn)行單端測(cè)序，測(cè)序長(zhǎng)度為50bp。對(duì)獲得的原始數(shù)據(jù)首先進(jìn)行過濾，去除接頭序列、低質(zhì)量序列、含N比例>10%的序列、5接頭污染的序列、沒有3接頭序列和插入片段的序列及含連續(xù)的/T/G/C的序列；以及長(zhǎng)度異常的序列。

1.3miRNA鑒定

首先將過濾后的數(shù)據(jù)比對(duì)到miRbase數(shù)據(jù)庫(kù)（V22）鑒定已知miRNA；將未比對(duì)至miRBase的數(shù)據(jù)比對(duì)至Rfam數(shù)據(jù)庫(kù)中的ncRNA序列（rRNA、tRNA、snRNA、snoRNA）和參考基因組重復(fù)序列；隨后對(duì)于上述均未比對(duì)中的cleanreads，使用mirEvo和miRdeep2進(jìn)行新miRNA預(yù)測(cè)。對(duì)于已知和新的miRNA序列，將其前體序列比對(duì)至miRBase數(shù)據(jù)庫(kù)進(jìn)行microRNA家族的鑒定。

1.4表達(dá)水平分析

為了探究毛竹種子萌發(fā)露白階段miRNA的表達(dá)，首先統(tǒng)計(jì)各樣本中已知和新miRNA的readscount．進(jìn)行TPM（Transcripts per million）歸一化轉(zhuǎn)化（TPM=readscount×107／library size）。隨后使用edgeR軟件進(jìn)行差異表達(dá)分析，以鑒定露白階段毛竹種子響應(yīng)干旱或鹽脅迫的miRNA，差異表達(dá)miRNA篩選條件為llog2foldchangel≥1且p<0.05。

1.5靶基因預(yù)測(cè)及功能富集分析

為了探究毛竹萌發(fā)種子中miRNA的潛在功能，對(duì)本研究中已知和新miRNA的靶基因進(jìn)行預(yù)測(cè)，并進(jìn)行靶基因的KEGG和GO功能富集；同時(shí)通過對(duì)差異表達(dá)miRNA靶基因的功能富集，進(jìn)一步分析毛竹萌發(fā)種子中miRNA響應(yīng)干旱或鹽脅迫的潛在通路或途徑。使用軟件TargetFinder和psRNAtarget進(jìn)行miRNA靶基因的預(yù)測(cè)，并對(duì)預(yù)測(cè)結(jié)果取交集。使用OmicShare在線平臺(tái)（（https：／／www．omicshare．com／tooIs））進(jìn)行靶基因KEGG和GO功能富集，富集顯著性閾值為p<0.05。

1.6差異表達(dá)miRNA驗(yàn)證

總RNA提取方法同1.2，反轉(zhuǎn)錄使用miRcute增強(qiáng)型miRNA cDNA第一鏈合成試劑盒（KR211，TlANGEN），熒光定量檢測(cè)使用增強(qiáng)型熒光定量試劑盒（FP411，TIANGEN）。熒光定量正向引物見表1，反向引物使用試劑盒通用引物，內(nèi)參引物為U6。熒光定量檢測(cè)儀器為qTOWER3.0（AnaIytik jene），每個(gè)基因進(jìn)行3次生物學(xué)重復(fù)和3次技術(shù)重復(fù)，采用2方法計(jì)算相對(duì)表達(dá)水平。

2結(jié)果與分析

2.1small RNA測(cè)序結(jié)果

small RNA測(cè)序結(jié)果表明：每個(gè)文庫(kù)產(chǎn)出數(shù)據(jù)約0.9～1.8GB，過濾后每個(gè)樣本中分別有12984751條（A， H20）、14982423條（B，10%PEG）、13401719條（C，15% PEG）、13995450條（D，50mmoI.L-1NaCI）和11452796條（E，100mmoI.L-1NaCI）Clean Reads（表2）.

CleanReads的Q20值約為98%，Q30值均≥93%，與毛竹參考基因組的比對(duì)率≥88.42%，表明測(cè)序數(shù)據(jù)質(zhì)量較好。具體數(shù)據(jù)產(chǎn)出及比對(duì)情況見表2。

2.2miRNA鑒定與miRNA家族分析

對(duì)small RNA數(shù)據(jù)庫(kù)中miRNA進(jìn)行比對(duì)和預(yù)測(cè)，結(jié)果顯示：每個(gè)樣品中有0.08%～0.11%的序列被注釋為已知miRNA，0.10%～0.11%的序列被預(yù)測(cè)為novel miRNA。對(duì)miRNA長(zhǎng)度進(jìn)行統(tǒng)計(jì)發(fā)現(xiàn)本研究成熟miRNA長(zhǎng)度主要分布在21nt和24nt（圖2）。本研究共鑒定miRNA成熟體序列508條，前體序列845條，其中，已知miRNA成熟體序列數(shù)量為246，前體序列數(shù)量為574，鑒定出novel miRNA成熟體數(shù)量共有262條，前體序列數(shù)量為271。

對(duì)miRNA前體序列進(jìn)行家族分析，共鑒定到包含514條前體序列的45個(gè)家族，包括MIR159、MIR166、MIR156、MIR408、MIR399、MIR530等，其中，有10條新miRNA前體序列也被鑒定屬于已知家族，如novel 57屬于MIR408，novel 14和novel 66屬于MIR169 2。統(tǒng)計(jì)每個(gè)家族包含的前體成員數(shù)目表明，本研究毛竹萌發(fā)種子中最大的miRNA家族為MIR159，包含54個(gè)成員，其次為MIR166 （48）、MIR156 （46）、MIR167 1（37）、MIR396（32）等（詳細(xì)數(shù)據(jù)未列出）。

2.3miRNA表達(dá)水平分析

毛竹萌發(fā)種子中表達(dá)量前5的已知miRNA分別為phe-miR166a-3p、phe-miR159a.1、phe-miR319a-3p.2-3p、phe-miR6478、phe-miR156a-5p，表達(dá)量前5的新miRNA分別為novel 1、novel 199、novel 311、novel 241、novel 2，其中，novel 1在5個(gè)樣本所有miRNA中表達(dá)量均為最高。表達(dá)量前10的新miRNA折疊最小自由能值為-187.0～-13.8

KJ·mol-1（表3、4）。

對(duì)每個(gè)樣本中表達(dá)水平前10的已知miRNA和新miRNA的進(jìn)行TPM值統(tǒng)計(jì)，結(jié)果表明：排名前10的已知miRNA在每個(gè)樣本總miRNA中的讀序豐度（total miRNA reads）占45%～51.8%，排名前10的新miRNA在每個(gè)樣本總miRNA中的讀序豐度中占23.1%～29.5%。PEG和NaCI脅迫下這些miRNA在毛竹露白階段種子表達(dá)水平均較高，推測(cè)其在毛竹種子萌發(fā)的調(diào)控中可能具有保守的重要作用。

2.4差異表達(dá)分析

通過差異表達(dá)分析，本研究共鑒定181個(gè)差異miRNA，其中，有84個(gè)上調(diào)表達(dá)，97個(gè)下調(diào)表達(dá)（圖3a）。與對(duì)照組（A）相比，10%PEG（B）、15%PEG（C）、50mmoI.L-1NaCI（D）、100mmoI.L-1NaGI（E）

4種脅迫下分別有20、26、41和24個(gè)差異表達(dá)miRNA（DEmiRNA），此外在B-C和D-E比較組分別有16和54個(gè)DEmiRNA（圖3b）。

依據(jù)miRNA的差異表達(dá)情況和表達(dá)水平，對(duì)本研究具有高豐度且顯著差異表達(dá)的miRNA進(jìn)行聚焦。與對(duì)照組相比，10% PEG、15% PEG、50mmoI.L-1 NaCI、100mmoI.L-1 NaCI中表達(dá)水平最高的DEmiRNA分別為novel 14、novel 311、novel 14、phe-miR159a.1，表達(dá)水平最高的已知DEmiRNA分別為phe-miR171e-5p、phe-miR3630-3p、phe-miR171e-5p和phe-miR159a.1，這些miRNA在毛竹種子萌發(fā)露白階段中大量積累，同時(shí)在PEG或NaCI脅迫下差異表達(dá)，可能參與了種子萌發(fā)露白階段miRNA對(duì)干旱或鹽脅迫的調(diào)控。

2.5miRNA靶基因預(yù)測(cè)

通過Targetfinder預(yù)測(cè)，共獲得505個(gè)miRNA的31729個(gè)靶基因，psRNAtarget預(yù)測(cè)結(jié)果中501個(gè)miRNA可以靶向22215個(gè)基因。對(duì)2個(gè)軟件的預(yù)測(cè)結(jié)果取交集后共有17666對(duì)miRNA-靶基因關(guān)系，包含489個(gè)miRNA和8812個(gè)靶基因。預(yù)測(cè)的miRNA-靶基因關(guān)系中，novel 153的靶基因數(shù)目最多，有188個(gè)靶基因，其次是phe-miR396b（167個(gè)）、phe-miR396e-5p（161個(gè)）、phe-miR396h （141個(gè)）、phe-miR164b（138個(gè)）.對(duì)同- miRNA家族成員的靶基因數(shù)目進(jìn)行統(tǒng)計(jì)，靶基因數(shù)目前十的家族中MIR396家族靶基因數(shù)目最多。此外，根據(jù)預(yù)測(cè)結(jié)果8812個(gè)靶基因中每個(gè)基因能夠受到1～28個(gè)miRNA靶向，表明在毛竹種子萌發(fā)過程存在復(fù)雜的miRNA調(diào)控網(wǎng)絡(luò)。具有miRNA-gene關(guān)系最多的基因中，排名前20的基因有7個(gè)屬于GAMYB（GibberellinM、YB）家族，飛2個(gè)屬于SPL（Squamosapromoter-binding-like）家族，推測(cè)GAM、YB和SPL基因家族在毛竹萌發(fā)種子miRNA調(diào)控網(wǎng)絡(luò)中可能具有重要意義。

2.6差異表達(dá)miRNA靶基因功能富集分析

根據(jù)靶基因預(yù)測(cè)結(jié)果，A-B、A-C、A-D、A-E、B-C、D-E 6個(gè)比較組差異表達(dá)miRNA分別有609、548、1011、812、460和1866個(gè)靶基因。GO富集結(jié)果顯示：6個(gè)組合差異表達(dá)miRNA的靶基因分別顯著富集在520（A-B）、516（A-C）、620（A-D）、446（A-E）、483（B-G）、610（D-E）個(gè)條目中，總計(jì)1989個(gè)涉及細(xì)胞成分、生物學(xué)過程和分子功能的條目，推測(cè)miRNA可以通過這些靶基因廣泛地參與不同的生物學(xué)過程或分子功能等調(diào)控。表5為不同比較組GO富集前10的生物學(xué)過程條目，與對(duì)照組相比，10% PEG（B）、15% PEG（C）．50mmoI.L-1 NaCI（D）、100mmoI.L-1 NaCI（E）脅迫下DEmiRNA靶基因富集最顯著的生物學(xué)過程分別為跨膜轉(zhuǎn)運(yùn)（GO：0055085）、DNA復(fù)制正向調(diào)控（G0：0045740）、發(fā)育過程（G0：0032502）、海藻糖生物合成（GO：0005992），推測(cè)這些生物學(xué)過程對(duì)于miRNA調(diào)控PEG和NaCI脅迫下的種子萌發(fā)具有重要作用。

KEGG富集結(jié)果（表6）表明：苯丙烷生物合成途徑在A-B、A-C、A-E 3個(gè)組合中被顯著富集，推測(cè)苯丙烷生物合成途徑對(duì)于毛竹種子萌發(fā)露白階段miRNA響應(yīng)PEG和NaCI脅迫可能具有重要意義，但毛竹種子萌發(fā)過程miRNA與苯丙烷合成途徑關(guān)鍵基因的靶向關(guān)系及其是否共同調(diào)控干旱或鹽脅迫下的種子萌發(fā)需要進(jìn)一步的研究和驗(yàn)證。此外，果糖和甘露糖代謝、玉米素生物合成途徑僅在鹽脅迫下（A-D和A-E）顯著富集，油菜素內(nèi)酯生物合成和脂肪酸降解途徑僅在A-B中顯著富集，表明在毛竹種子露白階段miRNA響應(yīng)10%PEG、15% PEG、50mmoI.L-1NaCI、100mmoI.L-1 NaCI脅迫的調(diào)控通路也存在差異。

2.7差異表達(dá)miRNA qPCR驗(yàn)證

隨機(jī)選取10個(gè)在不同樣本間存在差異表達(dá)的miRNA進(jìn)行qPCR驗(yàn)證，miRNA的表達(dá)量在不同處理組之間被上調(diào)或下調(diào)，如phe-miR6478在10% PEG脅迫下比對(duì)照組上調(diào)1.45倍，phe-miR3630-3p在100mmoI.L-1 NaCI脅迫下上調(diào)1.33倍，phe-miR171e-5p在10% PEG、15% PEG、50 mmoI.L-1 NaCI、100mmoI.L-1 NaCI脅迫下分別下調(diào)0.23、0.50、0.56和0.61倍（圖4）。miRNA熒光定量結(jié)果和small RNA測(cè)序數(shù)據(jù)在表達(dá)趨勢(shì)上整體一致，表明測(cè)序數(shù)據(jù)的可靠性。

3討論

在本研究中，miR166、miR159、miR319、miR156、miR396、miR167、miR168、miR894、miR160等在對(duì)照組和不同濃度的PEG和NaCI脅迫條件下均具有較高表達(dá)水平，推測(cè)這些miRNA可能在毛竹種子萌發(fā)中有重要作用，其中，miR159、miR156、miR396、miR160已經(jīng)被證明參與了種子萌發(fā)調(diào)控。在大麥（Hordeum vulgare L）萌發(fā)種子中，miR156、miR166、miR167、miR168等也具有高表達(dá)水平，同時(shí)miR5071被大量積累，而在本研究中未檢測(cè)到miR5071的表達(dá)，表明不同物種間miRNA對(duì)種子萌發(fā)的調(diào)控具有特異性。在已有研究中，miR402、miR163和miR417被證明能夠參與種子萌發(fā)期的干旱或鹽脅迫響應(yīng)，而在本研究對(duì)照和處理組中均未檢測(cè)到miR402、miR163和miR417的表達(dá)，推測(cè)這些miRNA在本研究中可能不是響應(yīng)干旱或鹽脅迫的關(guān)鍵miRNA，其分布和表達(dá)可能與物種有關(guān)。

與對(duì)照組相比，本研究中miRNA主要在單一脅迫下顯著差異表達(dá)，僅有2個(gè)miRNA能夠同時(shí)響應(yīng)2種脅迫，表明miRNA在10% PEG、15%PEG、50mmoI.L-1 NaCI、100 mmoI.L-1 NaCI4種脅迫下的調(diào)控存在差異。本研究豐度較高且顯著差異表達(dá)的miRNA中，phe-miR171e-5p能夠響應(yīng)10% PEG和50mmoI.L-1 NaCI脅迫，phe-miR3630-3p在15%PEG脅迫下顯著上調(diào)表達(dá)，phe-miR159a.1在100mmoI.L-1 NaCI脅迫下差異表達(dá)，推測(cè)這3個(gè)miRNA對(duì)于毛竹種子萌發(fā)響應(yīng)PEG或NaGI脅迫可能具有重要意義。

前人的研究中，PEG脅迫下甘藍(lán)型油菜（Brassica napus L）萌發(fā)種子中的miR171比對(duì)照組顯著下調(diào)表達(dá)；此外，在擬南芥（Arabidopsisthaliana （L.） Heynh.）中過表達(dá)桑樹（Morus albaL_） mno-miR171可提高轉(zhuǎn)基因植株在NaCI和甘露醇脅迫下的種子發(fā)芽率，這與本研究phe-miR171e-5p能夠響應(yīng)PEG和NaCI脅迫是一致的。miR171主要通過靶向GRAS家族基因影響GA和ABA信號(hào)通路，參與胚胎發(fā)生潛能的維持、花藥發(fā)育、芽分枝和復(fù)葉形態(tài)的調(diào)控、頂端優(yōu)勢(shì)的調(diào)控等，其是否通過GRAS家族參與毛竹種子萌發(fā)過程中對(duì)干旱和鹽脅迫的調(diào)控需要進(jìn)一步研究。在苜蓿（Medicago sativaL．）中，miR3630被鑒定為響應(yīng)干旱脅迫的miRNA，本研究中phe-miR3630-3p在15% PEG脅迫下顯著上調(diào)表達(dá)，但目前關(guān)于miR3630的研究較少，本研究為miR3630響應(yīng)干旱脅迫提供了新的證據(jù)。此外，phe-miR159a.1在A-E比較組顯著下調(diào)。擬南芥中miR159通過介導(dǎo)GAMYB家族MYB101和MYB33的轉(zhuǎn)錄本切割，參與種子萌發(fā)過程糊粉層細(xì)胞的細(xì)胞程序性死亡，推測(cè)miR159對(duì)于100mmoI.L-1 NaCI脅迫下毛竹種子的萌發(fā)也可能具有重要的調(diào)控作用。

4結(jié)論

本研究首次系統(tǒng)鑒定了毛竹萌發(fā)種子中的miRNA，并對(duì)其在PEG和NaCI脅迫下毛竹萌發(fā)露白階段種子中的表達(dá)模式進(jìn)行研究，探究了響應(yīng)不同PEG和NaCI脅迫的差異表達(dá)miRNA。下一步將對(duì)本研究中涉及的重要miRNA進(jìn)行靶基因的驗(yàn)證，并對(duì)其調(diào)控機(jī)制和調(diào)控功能進(jìn)行深入研究和探討。