王曉靜,王 濤,楊 凱,李潞濱

(1. 中國(guó)林業(yè)科學(xué)研究院林業(yè)研究所,林木遺傳育種國(guó)家重點(diǎn)實(shí)驗(yàn)室,國(guó)家林業(yè)和草原局林木培育重點(diǎn)實(shí)驗(yàn)室, 北京 100091;2.運(yùn)城學(xué)院,山西 運(yùn)城 044011;3.北京植物園,北京市花卉園藝工程技術(shù)研究中心/ 城鄉(xiāng)生態(tài)環(huán)境北京實(shí)驗(yàn)室,北京 100093;4.北京農(nóng)學(xué)院植物科學(xué)與技術(shù)學(xué)院,北京 100096)

毛竹(Phyllostachysedulis)是禾本科(Gramineae)竹亞科(Bambusoideae)剛竹屬(Phyllostachys)的多年生木本散生竹,是我國(guó)種植面積最廣、應(yīng)用及研究最為廣泛和深入的竹種,具有重要的經(jīng)濟(jì)、生態(tài)和文化價(jià)值[1-2]。使用毛竹種子進(jìn)行實(shí)生苗造林是毛竹造林的重要方式,具有竹鞭生長(zhǎng)較快、萌筍多、竹材利用率高、運(yùn)輸方便、成本低等優(yōu)點(diǎn),有利于提高造林整齊度、成活率和遺傳多樣性[3-5]。但毛竹種子壽命短,其種子萌發(fā)容易受到環(huán)境的影響,限制了實(shí)生苗造林及其研究[6-12]。因此,開展毛竹種子萌發(fā)相關(guān)研究對(duì)毛竹種質(zhì)資源的保存和應(yīng)用具有重要意義。

隨著全球氣候的極端變化,干旱和鹽脅迫等非生物脅迫變得越來(lái)越頻繁,增加了植物生長(zhǎng)和發(fā)育面臨的挑戰(zhàn)。前期研究結(jié)果表明,干旱和鹽脅迫是抑制毛竹種子萌發(fā)的重要環(huán)境因素[13],研究者已經(jīng)對(duì)干旱和鹽脅迫下毛竹種子的萌發(fā)率、活力指數(shù)、氧化酶活性以及萌發(fā)后幼苗的生長(zhǎng)開展了研究[6, 10, 12],但對(duì)分子水平上的相關(guān)調(diào)控機(jī)制仍然缺乏了解。因此,亟須挖掘調(diào)控毛竹種子萌發(fā)和萌發(fā)期抗旱抗鹽的因子,為竹類資源的遺傳改良、培育和應(yīng)用提供依據(jù)。

circRNAs是一類由轉(zhuǎn)錄本反向可變剪接產(chǎn)生的內(nèi)源非編碼RNA,在真核生物中非常豐富[14]。circRNAs的環(huán)狀結(jié)構(gòu)相對(duì)穩(wěn)定,能夠參與調(diào)節(jié)線性轉(zhuǎn)錄本表達(dá)、選擇性剪接以及與miRNA或蛋白質(zhì)的相互作用等[15]。近年來(lái)對(duì)植物中與circRNAs的研究逐漸興起,目前PlantcircBase數(shù)據(jù)庫(kù)(Version 7.0,更新于2022-03-15)已經(jīng)收錄了擬南芥(Arabidopsisthaliana)、水稻(Oryzasativa)、玉米(Zeamays)、小麥(Triticumaestivum)等21種植物的171 118個(gè)circRNAs[16]。據(jù)報(bào)道circRNAs參與了植物生長(zhǎng)、發(fā)育和脅迫響應(yīng)[17-18]。Zhou等[19]通過CRISPR-Cas9技術(shù)發(fā)現(xiàn)水稻circRNAOs05circ02465的無(wú)效突變體在種子萌發(fā)期間表現(xiàn)出較高的耐鹽性,表明其能夠參與脅迫下種子萌發(fā)的調(diào)控。竹類植物中,Wang等[20]對(duì)毛竹莖稈快速生長(zhǎng)階段的circRNAs進(jìn)行了鑒定和分析,Li等[21]報(bào)道了circRNAs在3月齡的毛竹幼苗中對(duì)干旱脅迫的響應(yīng)。然而,竹類植物中與circRNAs相關(guān)的研究還存在大量空白,特別是在干旱或鹽脅迫下circRNAs如何參與毛竹種子萌發(fā)調(diào)控研究鮮見報(bào)道。

本研究使用聚乙二醇(PEG)和NaCl模擬干旱和鹽脅迫,通過構(gòu)建鏈特異性文庫(kù)首次對(duì)毛竹萌發(fā)種子中的circRNAs進(jìn)行鑒定,對(duì)PEG和NaCl脅迫下毛竹種子萌發(fā)階段circRNAs的表達(dá)進(jìn)行分析,并對(duì)其源基因進(jìn)行功能富集,以期為深入開展竹類植物種子在非生物脅迫條件下萌發(fā)的分子調(diào)控機(jī)制研究奠定基礎(chǔ)。

1 材料與方法

1.1 試驗(yàn)材料

試驗(yàn)所用毛竹種子于2020年9月由廣西壯族自治區(qū)林業(yè)科學(xué)研究院收集,剝除外稃后挑選健康飽滿且大小均勻的種子備用(圖1a)。種子千粒質(zhì)量為(21.85±0.04)g,含水率為(12.05±0.35)%。將種子表面消毒后轉(zhuǎn)移到底部有兩層無(wú)菌濾紙的無(wú)菌培養(yǎng)皿中,每個(gè)培養(yǎng)皿放置20粒種子。試驗(yàn)分為5組(A1—E1):分別在培養(yǎng)皿中加入H2O(A1)、10%(質(zhì)量分?jǐn)?shù),下同) PEG(B1)、15% PEG(C1)、50 mmol/L NaCl(D1)和100 mmol/L NaCl(E1),每組萌發(fā)200粒種子,重復(fù)3次。在(26 ± 2) ℃的持續(xù)黑暗環(huán)境下進(jìn)行萌發(fā),每2 d更換濾紙和培養(yǎng)液。觀察和記錄種子萌發(fā)狀態(tài),在第4天選擇生長(zhǎng)狀態(tài)一致且到達(dá)種皮破裂階段的種子用于建庫(kù)(圖1b)。對(duì)于每個(gè)處理,在3個(gè)重復(fù)中各取5粒種子并將15粒種子混樣,使用液氮速凍后置于-80 ℃保存。

a)毛竹干種子dry seed;b)種皮破裂階段取樣的萌發(fā)種子sampled seed germinated at seed coat rupture stage。圖1 實(shí)驗(yàn)用毛竹種子Fig. 1 Moso bamboo seeds used in the experiment

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

使用植物多糖多酚RNA提取試劑盒(DP441,天根生物科技有限公司,中國(guó))提取樣品的總RNA,RNA質(zhì)檢合格后首先使用Ribo-ZeroTM試劑盒(MRZSR116,Epicentre,中國(guó))去除rRNA,隨后將RNA片段化并合成cDNA第1條鏈,在cDNA第2鏈合成時(shí)用dUTP取代dTTP。使用尿嘧啶特異性切除酶消化含有U的cDNA第2鏈,然后對(duì)第1鏈連接接頭后進(jìn)行PCR擴(kuò)增和富集得到最終的鏈特異性文庫(kù),并對(duì)其進(jìn)行初步定量。使用Illumina高通量測(cè)序平臺(tái)NovaSeq 6000(美國(guó))對(duì)質(zhì)檢合格的鏈特異性文庫(kù)進(jìn)行雙末端測(cè)序,測(cè)序長(zhǎng)度為150 bp。對(duì)測(cè)序后的原始數(shù)據(jù)進(jìn)行過濾,過濾掉接頭序列、含有≥1%的未知堿基的序列、低質(zhì)量(≤20)堿基數(shù)≥50%的序列。

1.3 circRNA鑒定

使用find_circ[22]和CIRI2[23]鑒定circRNAs,find_circ查找接頭序列與線性序列位置相反,且內(nèi)部存在GT-AG剪接位點(diǎn)的序列;CIRI2尋找具有目的信號(hào)和GT-AG剪切位點(diǎn)信息的連接點(diǎn)處的讀序(junction reads)序列。然后根據(jù)與參考基因組的比對(duì)結(jié)果[24-25],對(duì)鑒定結(jié)果取交集,并保留至少1個(gè)樣品中junction reads數(shù)≥5的circRNAs。

1.4 表達(dá)水平分析

使用具有circRNA剪接位點(diǎn)的junction reads數(shù)目估計(jì)circRNAs的表達(dá)量,將表達(dá)量進(jìn)行SRPBM(spliced reads per billion mapping)歸一化定量[26]。使用edgeR進(jìn)行差異表達(dá)分析[27],差異表達(dá)circRNAs篩選條件為差異倍數(shù)(foldchange)以2為底的對(duì)數(shù)絕對(duì)值≥1且P<0.05。使用TBtools繪制表達(dá)量熱圖[28]。

1.5 源基因及其功能富集分析

使用HISAT2將circRNAs比對(duì)至毛竹染色體水平參考基因組,鑒定其源基因[24-25]。對(duì)源基因進(jìn)行GO和KEGG注釋,P<0.05為富集結(jié)果顯著性閾值[29-30],使用ImageGP在線平臺(tái)展示富集結(jié)果[31]。

2 結(jié)果與分析

2.1 測(cè)序數(shù)據(jù)產(chǎn)量

本研究得到lncRNA原始數(shù)據(jù)如表1所示,過濾后共產(chǎn)出14.66～15.28 Gb數(shù)據(jù),clean reads的Q20值均≥98%,Q30值均≥95%,GC含量為52.4%～53.7%,將過濾后數(shù)據(jù)進(jìn)行比對(duì),77.70%～84.98%的數(shù)據(jù)可以比對(duì)到基因組,表明數(shù)據(jù)產(chǎn)量較高,且質(zhì)量較好,能夠滿足后續(xù)分析。測(cè)序原始數(shù)據(jù)提交至NCBI數(shù)據(jù)庫(kù)(PRJNA545401)。

表1 文庫(kù)測(cè)序數(shù)據(jù)統(tǒng)計(jì)

2.2 circRNA的鑒定

將兩個(gè)軟件的鑒定結(jié)果取交集,結(jié)果在本研究不同處理下萌發(fā)的毛竹種子中共鑒定到1 446個(gè)circRNAs,其長(zhǎng)度分布范圍為173～99 647 bp,每個(gè)樣本中長(zhǎng)度≤2 000 bp的circRNAs頻率均為最高(圖2a)。1 446個(gè)circRNAs主要分布在毛竹染色體scaffold 1—24,分布在scaffold 21和scaffold 17的數(shù)目最多,而分布在scaffold 1、scaffold 2、scaffold 19染色體的circRNAs最少,此外有少量(11個(gè))circRNAs分布在基因組未組裝片段中(圖2b)。1 446個(gè)circRNAs中有76.34%來(lái)源于基因組外顯子區(qū),19.23%來(lái)源于基因間隔區(qū),來(lái)源于內(nèi)含子區(qū)的circRNAs最少(4.42%)。每個(gè)樣本分布在基因組不同來(lái)源區(qū)域circRNAs的數(shù)量見圖2c,結(jié)果顯示每個(gè)樣本中來(lái)源于外顯子區(qū)的circRNAs數(shù)量最多,只有少量circRNAs來(lái)源于內(nèi)含子區(qū)。

a)5個(gè)樣本circRNAs長(zhǎng)度分布與頻率the length and frequency of circRNAs in five samples;b)circRNAs在毛竹染色體的分布distribution of circRNAs in moso bamboo chromosomes;c)circRNAs在基因組的分布distribution of circRNAs in genome。圖2 毛竹萌發(fā)種子中circRNAs的特征Fig. 2 Characteristics of circRNAs in germinated seeds of moso bamboo

2.3 源基因鑒定

對(duì)circRNAs的源基因進(jìn)行分析,結(jié)果表明,產(chǎn)自外顯子區(qū)和內(nèi)含子區(qū)的1 168個(gè)circRNAs共來(lái)源于867個(gè)基因,有698個(gè)源基因產(chǎn)生了1個(gè)circRNA,176個(gè)源基因產(chǎn)生了2個(gè)及以上circRNAs,14個(gè)源基因產(chǎn)生了5～13個(gè)circRNAs。如胼胝質(zhì)合成酶基因PH02Gene21577(基因編號(hào),后同)(callose synthase 9,CALS9)產(chǎn)生了13個(gè)circRNAs,核孔錨蛋白基因PH02Gene01074(nuclear-pore anchor,NUA)和CHD3型染色質(zhì)重塑因子PICKLE編碼基因 PH02Gene11133(CHD3-type chromatin-remodeling factor PICKLE,PKL)產(chǎn)生了9個(gè)circRNAs,胼胝質(zhì)合成酶基因PH02Gene00919(callose synthase 9,CALS9)產(chǎn)生了8個(gè)circRNAs。其中PKL能夠參與滲透脅迫和外源對(duì)種子萌發(fā)的調(diào)控,pkl突變體在脫落酸誘導(dǎo)下ABI3和ABI5過量表達(dá),其啟動(dòng)子區(qū)H3K9和H3K27甲基化水平比野生型大大降低,但不影響其在萌發(fā)后期幼苗的甲基化水平,同時(shí)PKL能夠沉默轉(zhuǎn)基因逆境誘導(dǎo)型RD29A啟動(dòng)子的轉(zhuǎn)錄沉默[32-33],但目前關(guān)于circRNAs在PKL調(diào)控中的功能仍是未知的,毛竹萌發(fā)期種子中PKL基因產(chǎn)生的circRNAs是否能夠參與PEG和NaCl脅迫下種子萌發(fā)的調(diào)控仍需進(jìn)一步研究。

2.4 表達(dá)特性分析

對(duì)照組(A1)毛竹萌發(fā)種子中表達(dá)量最高的circRNA為plant_circ_0001703,其源基因編號(hào)為PH02Gene18257,編碼油質(zhì)蛋白(olein)。10% PEG(B1)、15% PEG(C1)、50和100 mmol/L NaCl(E1)脅迫下毛竹萌發(fā)種子中表達(dá)量最高的circRNA分別為plant_circ_0001728、plant_circ_0000358、plant_circ_0001236、plant_circ_0001728,其源基因分別為球蛋白編碼基因PH02Gene28357、水通道蛋白TIP3-1編碼基因PH02Gene14831、功能未知基因PH02Gene04721和PH02Gene28357。對(duì)照組或脅迫處理下毛竹萌發(fā)種子中表達(dá)量前10的circRNAs熱圖及其源基因的注釋見圖3,其中包含了大量種子儲(chǔ)藏蛋白編碼基因,如球蛋白基因、vicilin-like種子儲(chǔ)藏蛋白基因、油質(zhì)蛋白基因等,表明在毛竹種子萌發(fā)過程中種子儲(chǔ)藏蛋白基因產(chǎn)生并積累了大量的circRNAs,但這部分circRNAs是否能夠調(diào)控毛竹種子萌發(fā)仍需要進(jìn)一步研究。

差異表達(dá)circRNAs(圖4)可見,與對(duì)照組相比,10%與15% PEG、50與100 mmol/L NaCl脅迫下差異表達(dá)circRNAs(differential expressed circRNAs,DECs)數(shù)量分別為524、505、467和474,其中有122個(gè)circRNAs在4種脅迫下均顯著差異表達(dá);此外,有503個(gè)circRNAs在PEG處理B1-C1比較組中差異表達(dá),465個(gè)circRNAs在不同濃度NaCl脅迫下差異表達(dá)(D1—E1);6個(gè)比較組中共有1 056個(gè)circRNAs差異表達(dá)。上述circRNAs在不同脅迫條件下萌發(fā)的毛竹種子中存在表達(dá)量的顯著變化,推測(cè)可能參與了脅迫下毛竹種子萌發(fā)的調(diào)控。

2.5 差異表達(dá)circRNAs源基因功能富集分析

圖5 不同比較組差異表達(dá)circRNAs源基因GO顯著富集前20結(jié)果Fig. 5 Top20 significantly enriched GO terms for original genes generating DECs in differential comparision groups

圖6 不同比較組差異表達(dá)circRNAs源基因KEGG顯著富集結(jié)果Fig. 6 Significantly enriched KEGG pathways of original genes generating DECs in different comparison groups

為了進(jìn)一步探究PEG和NaCl脅迫下毛竹萌發(fā)期種子中circRNAs的響應(yīng),對(duì)差異表達(dá)circRNAs(DECs)的源基因進(jìn)行GO和KEGG功能富集。6個(gè)比較組中顯著富集前20的GO條目見圖5,其中作用于酸酐的水解酶活性(hydrolase activity, acting on acid anhydrides,GO:0016817)在6個(gè)比較組中均被顯著富集到前20,推測(cè)其對(duì)毛竹萌發(fā)期種子中circRNAs對(duì)PEG和NaCl的響應(yīng)可能具有重要意義。KEGG富集結(jié)果(圖6)表明,6個(gè)比較組中顯著富集途徑數(shù)目分別為4(A1-B1)、3(A1-C1)、2(A1-D1)、5(A1-E1)、5(B1-C1)、1(D1-E1)。KEGG途徑中鞘脂代謝途徑(sphingolipid metabolism)在A1-B1、A1-C1、A1-E1 3個(gè)比較組均被顯著富集,推測(cè)鞘脂代謝途徑是circRNAs響應(yīng)PEG和NaCl的重要途徑。此外,過氧化酶體途徑,脂肪酸降解途徑,基本轉(zhuǎn)錄因子途徑在A1-C1和B1-C1比較組中被顯著富集,淀粉與蔗糖代謝在A1-C1比較組中被顯著富集,糖胺聚糖降解途徑在A1-E1、D1-E1比較組中均被顯著富集,ABC transporters在A1-E1比較組中被顯著富集,這些途徑與種子萌發(fā)密切相關(guān),推測(cè)在PEG或NaCl脅迫下毛竹種子萌發(fā)時(shí),上述途徑可能通過相關(guān)circRNAs響應(yīng)脅迫,并調(diào)控萌發(fā)進(jìn)程。

3 討 論

本研究在對(duì)照組、PEG和NaCl脅迫處理下萌發(fā)的毛竹種子中共鑒定1 446個(gè)circRNAs,對(duì)其特征分析表明這些circRNAs來(lái)源主要為外顯子區(qū),只有少部分circRNAs來(lái)源于內(nèi)含子區(qū),這與毛竹幼苗葉片中circRNAs的基因組來(lái)源區(qū)域是一致的,在擬南芥、玉米、水稻、小麥中也發(fā)現(xiàn)了相似的結(jié)果[34-37];同時(shí)本研究中circRNAs在染色體的分布是不均勻的,這與毛竹幼苗葉片中circRNAs的分布結(jié)果也是一致的。

已有研究[8,35,38-39]表明,植物中circRNAs參與了對(duì)干旱、鹽脅迫的響應(yīng)。在玉米368個(gè)自交系中,circRNAs的源基因富集在與circRNAs表達(dá)和耐旱性相關(guān)的SNP,表明了circRNAs在玉米干旱響應(yīng)中的重要作用[8]。擬南芥中過表達(dá)由保衛(wèi)細(xì)胞外向整流K+通道基因產(chǎn)生的circGORK導(dǎo)致轉(zhuǎn)基因植株對(duì)脫落酸高度敏感,但對(duì)干旱不敏感,表明circGORK在耐旱性中發(fā)揮了積極作用[35]。在紫花苜蓿(Medicagosativa)、鷹嘴豆(Cicerarietinum)和大豆(Glycinemax)中,大量的circRNAs能夠響應(yīng)干旱或鹽脅迫[38-39]。本研究在不同比較組中也鑒定了大量響應(yīng)PEG和NaCl脅迫處理的差異表達(dá)circRNAs,這些circRNAs對(duì)于干旱和鹽脅迫下的毛竹種子萌發(fā)可能具有重要調(diào)控作用。

對(duì)不同比較組差異表達(dá)circRNAs的源基因進(jìn)行GO富集分析,結(jié)果表明作用于酸酐的水解酶活性的相關(guān)基因在6個(gè)比較組中均顯著富集;涉及到的源基因包括CHD3型染色質(zhì)重塑因子PICKLE編碼基因(PKL)PH02Gene11692、PH02Gene11133、PH02Gene11133,ABC轉(zhuǎn)運(yùn)蛋白編碼基因(ABC transporter C,ABCC)PH02Gene23726、PH02Gene32020,鈣離子轉(zhuǎn)運(yùn)ATP酶基因(Calcium-transporting ATPase,ACA)PH02Gene33073、PH02Gene34105等,其中PH02Gene11133是本研究中產(chǎn)生circRNAs最多的源基因。前人的研究中PKL、ABCC和ACA均參與了種子萌發(fā)的調(diào)控[33, 40-41],推測(cè)作用于酸酐的水解酶活性相關(guān)基因可能對(duì)circRNAs調(diào)控干旱和鹽脅迫下的毛竹種子萌發(fā)具有重要意義。

根據(jù)差異表達(dá)circRNAs源基因的KEGG富集分析,鞘脂代謝途徑在A1-B1、A1-C1、A1-E1中均被顯著富集,推測(cè)其可能參與了PEG和NaCl脅迫下circRNAs對(duì)毛竹種子萌發(fā)的調(diào)控。鞘脂代謝途徑參與了玉米種子萌發(fā)過程中病原真菌誘導(dǎo)的細(xì)胞程序性死亡[42]和擬南芥種子萌發(fā)期對(duì)冷脅迫的響應(yīng)[43],但目前關(guān)于鞘脂代謝途徑參與種子萌發(fā)和萌發(fā)期抗逆的研究相對(duì)較少,鞘脂代謝途徑是否能夠通過circRNAs參與PEG和NaCl脅迫下的種子萌發(fā)仍需進(jìn)一步研究。本研究中DECs源基因中還涉及基礎(chǔ)轉(zhuǎn)錄因子,基礎(chǔ)轉(zhuǎn)錄因子在A1-C1,B1-C1中顯著富集,涉及的7個(gè)基因中6個(gè)屬于轉(zhuǎn)錄起始因子TFIID(transcription initiation factor, TFIID),推測(cè)TFIID可能參與了circRNA對(duì)毛竹種子萌發(fā)的調(diào)控。前人研究表明TFIID在干旱脅迫下對(duì)水稻生長(zhǎng)和光合作用的調(diào)節(jié)可能起著關(guān)鍵作用[44]。在擬南芥中TAF10參與構(gòu)成TFIID復(fù)合物,atTAF10的過表達(dá)提高了種子萌發(fā)期間對(duì)鹽脅迫的耐受性,而敲除突變體對(duì)鹽脅迫更為敏感,表明轉(zhuǎn)錄起始因子是植物鹽害毒性的生理靶標(biāo)復(fù)合體[45]。但目前關(guān)于TFIID通過circRNAs進(jìn)行基因調(diào)控的研究還未見報(bào)道,TFIID是否能夠通過產(chǎn)生circRNAs調(diào)控種子萌發(fā)還需要進(jìn)一步研究。

本研究首次對(duì)毛竹萌發(fā)期種子中的circRNAs進(jìn)行了初步鑒定和脅迫響應(yīng)分析,豐富了毛竹circRNAs數(shù)據(jù),并為脅迫下毛竹種子萌發(fā)調(diào)控因子的挖掘及相關(guān)研究提供參考;后續(xù)將進(jìn)一步開展circRNAs功能的驗(yàn)證工作,并對(duì)其在毛竹種子萌發(fā)中的調(diào)控機(jī)制進(jìn)行深入研究。