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

王曉靜 王濤 楊凱 李潞濱

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

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

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

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

1材料與方法

1.1試驗材料

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

1.2建庫、測序與數據過濾

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

1.3miRNA鑒定

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

1.4表達水平分析

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

1.5靶基因預測及功能富集分析

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

1.6差異表達miRNA驗證

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

2結果與分析

2.1small RNA測序結果

small RNA測序結果表明：每個文庫產出數據約0.9～1.8GB，過濾后每個樣本中分別有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%，與毛竹參考基因組的比對率≥88.42%，表明測序數據質量較好。具體數據產出及比對情況見表2。

2.2miRNA鑒定與miRNA家族分析

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

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

2.3miRNA表達水平分析

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

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

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

2.4差異表達分析

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

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

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

2.5miRNA靶基因預測

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

2.6差異表達miRNA靶基因功能富集分析

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

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

2.7差異表達miRNA qPCR驗證

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

3討論

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

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

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

4結論

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