王業(yè)建，梁曉玲，阿布來提·阿布拉，韓登旭，楊 杰，郗浩江，劉 俊，李銘東

(新疆農(nóng)業(yè)科學院糧食作物研究所， 烏魯木齊 830091)

0 引 言

【研究意義】干旱能抑制植物的生長和發(fā)育，并導致植物的品質(zhì)和產(chǎn)量下降。了解植物對非生物脅迫的反應對提高作物生產(chǎn)力有重要意義。【前人研究進展】在長期的進化過程中，植物已經(jīng)開發(fā)出一系列調(diào)節(jié)機制，以應對不同層面的這些不利條件，包括細胞、生理、生物化學和分子過程[1-2]。已確定植物中激素介導的信號交叉參與對干旱和鹽脅迫的反應，如脫落酸(ABA)、水楊酸和乙烯。基因表達調(diào)控是植物在轉錄后水平對抗干旱的關鍵策略[3-4]?；虮磉_調(diào)控最重要的參與者之一是microRNA(miRNA)。miRNA是一種長度約為21個核苷酸的非編碼小RNA分子。研究表明，miRNA具有與其靶mRNA的完美或接近完美的互補核苷酸結合的特點，通過mRNA切割或通過翻譯抑制負向調(diào)節(jié)靶基因的表達[5]。除了miRNA在發(fā)育和代謝中的作用外，miRNA還參與非生物和生物應激反應。例如，miR394是一種保守的miRNA，已在一系列植物物種中發(fā)現(xiàn)，如擬南芥、水稻、棉花、柳枝稷和甘藍型油菜[6-7]。目前，有多達40個與非生物脅迫相關的植物miRNA家族，其中多種與干旱脅迫反應有關[8]。玉米(Zeamays)基因組中miRNA的表征和分析取得了重大進展[9-10]。雖然玉米是一種相對耐旱和耐鹽的作物，但玉米高鹽度和過度缺水可導致一系列代謝紊亂，包括滲透作用(脫水)，營養(yǎng)不平衡和離子毒性，都會產(chǎn)生極大的負面影響[11]?；谵D錄組和轉基因分析，已顯示大量基因在玉米或其他物種中響應鹽度和干旱脅迫而顯示出持久的表達。例如，轉基因擬南芥中玉米CBL相互作用蛋白激酶基因(GhCIPK6)的過表達導致對鹽，干旱和ABA脅迫的耐受性提高，GlyCIP6可能是對抗玉米中鹽和干旱脅迫的正向調(diào)節(jié)因子[12]。轉基因煙草表達玉米C群絲裂素活化蛋白激酶(MAPK)基因(GhMPK2)具有較低的水分損失率，表現(xiàn)出對鹽和干旱的耐受性增強，GhMPK2可能正向調(diào)節(jié)煙草和玉米的鹽和耐旱性[13]?；谖㈥嚵械霓D錄組分析揭示了玉米中某些鹽/干旱介導的信號轉導途徑，其中許多候選基因表達差異，可能是耐受干旱脅迫的潛在標記，如WRKY，ERF，跨膜轉運蛋白，丙酮酸脫羧酶和蔗糖合成酶等[14]。【本研究切入點】與親本相比，近交親本系的F1自交后代表現(xiàn)出優(yōu)越的性能和脅迫耐受性[15-17]。這種效應被稱為雜種優(yōu)勢，在植物育種中被廣泛利用。研究運用現(xiàn)代分子生物學及高通量測序技術，測序鑒定玉米雄穗花器官分化期響應干旱脅迫的miRNA及其靶基因?！緮M解決的關鍵問題】通過miRNA測序，篩選出耐旱自交系“PHBA6”和干旱敏感自交系“吉63”之間，以及PHBA6干旱脅迫和對照處理之間差異表達的miRNA，挖掘調(diào)控耐旱自交系“PHBA6”雄穗花器官分化期干旱適應性的miRNA。

1 材料與方法

1.1 材 料

材料為耐旱自交“PHBA6”(PHZ51×PHG47)和干旱敏感自交系“吉63”[(127-32×鐵84)(W24×W20)輻)] , 由新疆農(nóng)業(yè)科學院糧食作物研究所提供。

材料于2017年在新疆農(nóng)科院海南三亞科技示范園基地大棚分期播種，前期進行正常的滴灌澆水處理，2個自交系的雄穗處于花器官分花期時進行干旱脅迫處理，耐旱種質(zhì)“PHBA6”進行正常滴灌澆水處理作為對照，干旱處理15 d后，分別取干旱處理的“PHBA6”和“吉63”的雄穗(分別命名為PDST15和JDST15)，以及對照處理的“PHBA6”的雄穗(命名為PNTT15)，迅速包好放入液氮中, 于-80℃保存。

1.2 方 法

1.2.1 總RNA提取及質(zhì)量控制

分別合并來自6個對照和PEG處理的玉米雄穗，使用miRNA Mini Kit(Qiagen)分離小RNA。 使用Agilent 2100生物分析儀(Agilent)測試RNA樣品的質(zhì)量和濃度。 將小RNA與3'和5'適配子連接，并使用TruSeq TM Small RNAkit(Illumina)根據(jù)說明書進行RT-PCR。 從凝膠中分離預期的最終PCR產(chǎn)物，純化并使用Hiseq2000測序儀(Illumina)測序。sRNA測序和基本生物信息學分析分別在MAGrogen(韓國)和LC Sciences(Texas, USA)上進行。

1.2.2 文庫制備和miRNA測序

用SolexaQA軟件包(Solexa)計算質(zhì)量統(tǒng)計數(shù)據(jù)，采用Illumina HiSeq 2500儀器進行深度測序得到50 bp 的單端測序讀數(shù)(reads)。來自4個庫中每個庫的生成數(shù)據(jù)文件用于后續(xù)分析。從miRNA讀數(shù)中去除適應序列后，建立小RNA的總讀數(shù)和獨特讀數(shù)。 從讀數(shù)中，映射到rRNA，tRNA，snoRNA和snRNA的小RNA被丟棄，剩余的讀數(shù)被映射到玉米基因組和轉錄本中。這些基因組合轉錄本數(shù)據(jù)可使用Bowtie軟件在美國國家醫(yī)學圖書館(NLM)的Genbank(http://ftp.ncbi.nlm. nih.gov)和法國的膜生物學國家重點實驗室獲得(http://www.biomemb. cnrs.fr/ date_ contigs.doc)。為了使用Bowtie進行對齊，需要去除缺乏3’的序列、5’污染的序列和小于18個核苷酸的序列。

1.2.3 保守和新miRNA的預測

為了鑒定保守的miRNA，使用Bowtie v 0.12.7將sRNA與NCBI中發(fā)表在miRBase 20.0(http：/ / www.miRBase.org)植物成熟miRNA比對。具有完整匹配序列的miRNA被認為是保守序列。具有1～3個錯配的miRNA被認為是miRNA變體，并且與miiRNA數(shù)據(jù)庫中沒有相似性的miRNA被認為是新的miRNA。 在該研究中鑒定的miRNA前體序列的長度范圍為57～264 bp。 使用UAEsmall RNA workbench軟件驗證了二級發(fā)夾結構的形成和推定的新miRNA序列在相應前體上的比對。

1.2.4 miRNA表達分析比較和功能

按照先前描述的方法，用miRNA Read/Clear×Total Read Count來進行表達量歸一化處理，即計算每個miRNA中的序列計數(shù)的標準化。文庫大小參數(shù)被認為是特定樣本和偽參考樣本的比率之間的中值。在去除了文庫中未顯示讀數(shù)的miRNA序列后，在文庫中技術miRNA獨特序列。使用數(shù)字基因表達數(shù)據(jù)軟件(degseq包)進行差異表達分析。將變異的生物系數(shù)(BCV)值設定為0.2。如果多個測試校正的P值并且錯誤發(fā)現(xiàn)率(FDR)<0.05，則預測樣品池間miRNA表達是差異顯著的。使用PermutMatrix軟件和KOBAS軟件進行層次聚類分析和KEGG 通路的富集分析。

1.2.5 miRNA靶基因的計算預測

使用psRNAtarget軟件在RNA-seq重疊群上搜索miRNA來完成對從4個文庫測序的保守的，變體和新的miRNA的靶基因的預測。參數(shù)設置為：最大期望值(ME)為3，完成度評分(hspsize)為20，目標位點周圍的側翼長度(上游17 bp、下游13 bp和中心錯配區(qū)9～11 nt的翻譯抑制為靶基因可分析范圍。

注釋潛在的mRNA靶標，并使用blast2GO v2.3.5軟件注釋功能，基于推斷的蛋白質(zhì)相似性，使用BLASTx在Swiss-Prot / Uniprot蛋白質(zhì)數(shù)據(jù)庫中獲得的其他序列和保守結構域。每個序列被注釋一個基因本體(GO)或更的細胞組分，分子功能和生物過程。

2 結果與分析

2.1 測序數(shù)據(jù)對比

研究表明，分別從JDST15、PDST15和PNTT15的miRNA文庫中獲得12133269、14669833和1009941 reads。去除低質(zhì)量的數(shù)據(jù)后，共剩下9632540、550036和6868703條有效讀取(valid reads)被用miRNA比對鑒定與預測分析。樣本中測序序列的分布概況，其中一些低比重占有率的大片段序列，例如rRNA或mRNA，樣本的總RNA質(zhì)量較高未被降解。3個文庫都顯示出與其他RNA家族相似的分布，包括rRNA(Uniq～1.34%)，sn RNA(Uniq～0.02%)，snoRNA(Uniq～0.01%)和tRNA(Uniq～0.13%)。表1

選取Rfam數(shù)據(jù)庫來注釋測序得到的小RNA序列，盡可能地發(fā)現(xiàn)并去除其中可能的rRNA、 scRNA、snoRNA、snRNA、tRNA等非miRNA序列。分別對總的測序數(shù)據(jù)(總數(shù))以及Unique數(shù)據(jù)(種數(shù))進行了統(tǒng)計。選取Rfam數(shù)據(jù)庫來注釋測序得到的小RNA序列，盡可能地發(fā)現(xiàn)并去除其中可能的rRNA、 scRNA、 snoRNA、 snRNA、tRNA等非miRNA序列。分別對總的測序數(shù)據(jù)(總數(shù))以及Unique數(shù)據(jù)(種數(shù))進行了統(tǒng)計。PNTT15組特異性rRNA、tRNA、snoRNA和snRNA所占比例分別為：54.50%、26.74%、6.39%和6.62%；PDST15組分別為為42.09%、26.09%、18.43%和10.97；JDST15組分別為53.12%、24.49%、15.58%和7.50%。圖1

圖1 Rfam數(shù)據(jù)比對Fig.1 Rfam data comparison results

2.2 Valid小RNA長度分布

研究表明，統(tǒng)計過濾后的Valid數(shù)據(jù)的總數(shù)及種數(shù)長度，大部分數(shù)據(jù)分布在20～24 nt，符合Dicer 酶切割的典型特征。圖2

圖2 Valid數(shù)據(jù)長度分布統(tǒng)計Fig.2 Valid data length distribution statistics

2.3 PDST15 vs PNTT15和JDST15 vs PDST15差異miRNA的檢出和發(fā)現(xiàn)

研究表明，在玉米中已經(jīng)報道的前體miRNA檢出176個。其中在JDST15、PDST15和PNTT15組中分別檢出160、151和145個。成熟體miRNA檢出199個在3組中分別檢出156、159和146個。在JDST15、PDST15和PNTT15組分別新發(fā)現(xiàn)前體miRNA 83、81和50個，成熟體miRNA7、76和42個。

根據(jù)miRBase(Release 20)(Kozomara和Griffiths-Jones，2011)中所有已知植物miRNA的3個文庫的所有清潔讀數(shù)的比對發(fā)現(xiàn)，miRNA在物種間高度保守，共有709個已知植物miRNA家族。這些miRNA家族占總獨特讀數(shù)序列的約0.47%和平均總冗余閱讀序列的23.59%。在這些miRNA家族中，71個miRNA在2個不同物種中響應干旱有特異型。47個miRNA家族僅在PHBA6自交系對照處理中發(fā)現(xiàn)。例如，miR1868和miR2099分別僅在干旱處理的樣品中表達。干旱處理的文庫共享65個miRNA家族，這些家族沒有出現(xiàn)在對照組文庫中。在3個文庫中共鑒定出709個miRNA家族中的357個，表明在維持正常生物活性中的關鍵作用，例如miR156 / 157、miR159、miR168和miR172。所有3個文庫都存在相似的最常見miRNA家族，包括miR156，miR157、miR166、miR167和miR3954。 Pearsonχ2檢驗顯示，709個(79.69%)miRNA家族中有565個在3個文庫中差異表達(P值≤0.05)。共有443個(61.37%)miRNA家族在三種處理的成對比較中表現(xiàn)出顯著差異，包括miR157，miR159，miR2948和miR3694。 miR1854和miR1148在對照與干旱和鹽與干旱之間具有最大的倍數(shù)變化，高達≥10倍，干旱脅迫強烈抑制其表達。表1～3

表1 3個文庫中總miRNA標簽分類Table 1 The classification of total small RNA tags in 3 library

2.4 miRNA靶基因預測和富集性

研究表明，共預測到4 658個基因，利用Fisher 精確假設檢驗對差異miRNA進行分析，得出miRNA和基因富集的KEGG通路和GO信號通路。在PDST15 vs PNTT15組中這些基因主要參與生物學過程中的蛋白磷酸化，轉錄本的調(diào)節(jié)，ARF蛋白信號轉導的調(diào)控；細胞組成中的維持細胞膜和核完整性；分子過程中的蛋白結合，ATP結合。在JDST15 vs PDST15組中，差異基因參與的生物學過程主要為：轉錄的調(diào)節(jié)，脂肪酸生物合成過程和蛋白磷酸化；細胞組成中的維持細胞膜和核完整性；分子過程中的蛋白結合，ATP結合以及DNA結合。

差異表達miRNAs靶向基因富集的通路主要為：碳水化合物形成，脂質(zhì)形成，氨基酸代謝產(chǎn)物，輔助因子和維生素生成，其他次生代謝的生物合成，聚糖生物合成和環(huán)境適應等。圖3～5

注：PHBA6干旱脅迫和對照前50種最豐富的保守miRNA(左)，不同品種玉米干旱脅迫后前50種最豐富的保守miRNA的熱圖。紅色，上調(diào); 綠色，下調(diào)

圖4 通過差異表達的miRNAs預測的靶基因注釋的部分基因GOFig.4 GO analysis of partial genes annotated by target genes predicted by differentially expressed miRNAs

圖5 預測差異表達miRNAs靶向的前10個途徑Fig.5 Top 10 pathways for predicting differentially expressed miRNAs

3 討 論

在干旱脅迫下，許多保守的和新的miRNA表達存在差異；一些miRNA甚至在干旱處理中特異性表達。在3個文庫中獲得了大量的miRNA讀數(shù)，同時在3個文庫中觀察到冗余讀數(shù)，獨特讀數(shù)和匹配的獨特讀數(shù)的類似大小分布，其中長度為24個核苷酸的讀數(shù)占大多數(shù)。匹配的冗余讀數(shù)在21個核苷酸類中具有最多的讀數(shù)，其次是24個核苷酸。玉米中的小RNA豐度和大小與擬南芥[18]和水稻[19]中報道的結果大致一致。

在所有3種處理中，miR156/157家族是最豐富的小RNA，占總小RNA含量的60%。與對照相比，miR156和miR157在干旱處理中分別下調(diào)0.44和1.23倍。據(jù)報道，玉米根和葉中的miR156和miR157均在高濃度鹽(> 2.5%)下表達量均下降，并且隨著PEG的增加以劑量依賴性方式下調(diào)。目前的在干旱脅迫下雄穗中的結果也與前人報道基本一致。研究表明，miR156 / 157在植物中負向地靶向SPL轉錄因子，并且miR156 / 157過表達導致擬南芥中成年性狀和開花的延遲[20]。miR156(Corngrass1，Cg1)的過量產(chǎn)生導致玉米幼年營養(yǎng)期的延長。目前關于miR156 / 157的研究主要集中在其在形態(tài)變化和開花調(diào)控中的作用[21]。研究中，小RNA測序提供了干旱脅迫干擾miR156 / 157表達的證據(jù)，miR156 / 157在干旱脅迫下的新作用。

大豆中NF-Y復合物的NF-YA(GmNFYA3)可被ABA和非生物脅迫誘導，包括干旱，NaCl和冷刺激[22]。擬南芥中GmNFYA3的過量表達導致葉片水分流失減少和干旱脅迫增加，并提高其對高鹽度和外源ABA的敏感性。煙草中的體內(nèi)試驗顯示miR169指導GmNFYA3 mRNA切割[23]。在玉米中，預測NF-YA3(contig16841和contig22907)是ghr-miR169的靶標。此外，與對照相比，干旱和處理中的ghr-miR169a / b / c表達顯著下調(diào)0.04至0.93倍。相反，ghr-miR169d / e / f / g在干旱處理中顯著上調(diào)0.04～0.85倍。Ghr-miR169i / j表達在干旱處理中受到抑制。推測至少ghr-miR169a / b / c可能通過作用于玉米中的NF-YA3而在抗旱中發(fā)揮積極作用。

植物水通道蛋白是大型主要內(nèi)在蛋白家族的一類，眾所周知，通過不同的小分子(包括水和其他小營養(yǎng)素)的生物膜發(fā)揮作用[24]。除了在水和養(yǎng)分的吸收和運輸中的作用外，還發(fā)現(xiàn)水通道蛋白涉及一系列非生物脅迫，例如干旱和冷脅迫。例如，轉基因煙草中質(zhì)膜水溶蛋白的過表達在有利的生長條件下提高了植物活力，但在干旱或鹽脅迫下沒有，因為通過質(zhì)膜水通道蛋白的共轉運水在干旱或鹽脅迫期間具有有害作用[25]。已發(fā)現(xiàn)從玫瑰中分離出的TIP型水通道蛋白基因Rh-TIP1受到乙烯和水缺乏處理的抑制[26]。確定20個miRNA可能靶向玉米中的22個水通道蛋白，其中兩個miRNA-tar-get對(ghr-miR4371和contig780;以及ghr-miR4371和BK007054.1)也通過降解組序列分析檢測到。水通道蛋白可能通過作為miRNA靶標(例如ghr-miR4371)參與對玉米中的干旱脅迫的響應。

過量表達miR393的轉基因水稻植物對鹽和堿處理更敏感，miR393是通過靶向非生物相關基因對鹽和堿脅迫響應的負調(diào)節(jié)因子[27]。在該研究中，獲得了相反的結果，其中ghr-miR393a / b / c / d / e表達通過干旱和鹽度處理上調(diào)。然而，玉米中預測的ghr-miR393靶標也是應激相關基因，以及幾種激素應答基因，包括NADPH：細胞色素P450還原酶(CPR1)，III類過氧化物酶(POX4)，蛋白質(zhì)AUXIN SIGNALING F-BOX 3，和TIR1(運輸抑制器響應1)。如果玉米中的miR393也是水稻中的負壓脅迫調(diào)節(jié)因子，一種可能的解釋是，玉米幼苗中miR393s在干旱脅迫下的上調(diào)有助于擴大玉米的脅迫信號，觸發(fā)通過信號轉導串擾來對抗干旱脅迫的更有效或更強大的途徑，可能是通過生長素相關途徑。

使用基于CitationRank的算法，在植物的干旱和鹽度脅迫的背景下編碼相關的編碼基因?；贕O分析允許確定基因屬于哪些GO術語(生物過程，分子功能和細胞成分)[28]。基于GO的分析可以提供更多關于理解miRNA功能的想法。研究中總共274個miRNA(256個保守的miRNA和18個新的miRNA)和其1驗364個靶標被分類為542個分子功能，679個生物學過程和138個細胞組分。至少137種miRNA及其與應激反應相關。229種靶標能夠分為1.2種分子功能，118種生物過程和37種細胞成分。35個miRNA-靶標對屬于響應于干燥的生物學過程(GO：0009269)和對水剝奪的反應(GO：0009414)，例如ghr-miR159b、ghr-miR166h、ghr-miR399i、 ghr-n26和ghr-miR399f。許多分類的生物過程與信號轉導有關，如生長素代謝(GO：0009850)和生物合成論(GO：0009851)，乙烯介導的信號通路(GO：0009873)，對生物刺激的反應(GO) ：0009607)，細胞分裂素代謝(GO：0009690)和生物合成(GO：0009691)，以及茉莉酸代謝(GO：0009694)和生物合成。其中，MYB3R被評為干旱環(huán)境中最重要的基因。最近，越來越多的研究報道MYB與植物的耐旱性有關。例如，煙草中的R2R3型MYB基因NbPhan被認為在解決干旱脅迫方面發(fā)揮了作用，其煙草沉默導致煙葉嚴重萎縮和增加水分流失率[29]。MYB對纖維發(fā)育也至關重要。玉米中miRNA和MYB轉錄因子之間的相關調(diào)節(jié)機制有望用于雙重目的，即纖維發(fā)育和耐旱性。對于鹽度反應，包括顆粒蛋白重復半胱氨酸蛋白酶，醇脫氫酶1，富含甘氨酸的RNA結合蛋白3和赤霉素2-β-雙加氧酶7的基因被確定為在鹽度響應的重要性方面排名較高。有許多報道顯示實際上顆粒蛋白重復半胱氨酸蛋白酶與干旱反應無關，而與植物中的程序性細胞死亡無關[30-31]。干旱脅迫可能更容易通過半胱氨酸蛋白酶觸發(fā)植物中的程序性細胞死亡。

4 結 論

鑒定了337前體miRNA，其中包含289種已知的miRNA和48種新的miRNA。在3個文庫，兩個分組中共有155種差異表達miRNA。目標預測，基于GO功能分類和基于遺傳和基因組(KEGG)的功能富集顯示，miRNA可能通過靶向一系列與脅迫相關的基因而在干旱脅迫中發(fā)揮作用。至少55個預測的靶基因進一步被60個miRNA調(diào)控。NAC、MYB和MAPK基因家族在干旱脅迫下評分最高，在植物抗旱中重要作用。根據(jù)目標基因預測，一系列玉米miRNA與這些排名靠前的基因相關，包括miR164、miR172、miR1520、miR6158、ghr-n24和ghr-n56等。miRNAs可能在玉米雄穗花器官分化期耐旱中發(fā)揮重要作用，miRNAs的篩選為分子輔助育種和轉基因育種將提供新的靶標。