60Co-γ射線輻射水稻小RNA測序變異分析

張 莉 張 勇 雷星宇 張淵海 張逸妍 彭選明 楊 震

(湖南省農(nóng)業(yè)科學(xué)院核農(nóng)學(xué)與航天育種研究所，湖南 長沙 410125)

隨著人類日益增長的糧食需求，水稻作為重要的糧食作物，其增產(chǎn)、穩(wěn)產(chǎn)對國民經(jīng)濟(jì)發(fā)展具有重要影響。輻射誘變育種技術(shù)利用X、γ、α、β射線或中子等照射生物體，改變其遺傳物質(zhì)，使后代出現(xiàn)新的變異，再從中篩選和培育新種質(zhì)和新品種。根據(jù)照射方式，輻射誘變可分為內(nèi)照射和外照射。目前外照射方法被廣泛應(yīng)用，處理對象包括各種作物的種子、花粉、幼穗等[1]。輻射誘變育種技術(shù)在我國農(nóng)作物突變種質(zhì)資源創(chuàng)制、新品種培育等應(yīng)用領(lǐng)域取得了一批重要成果，在保障國家糧食安全、推進(jìn)農(nóng)業(yè)綠色發(fā)展方面發(fā)揮了獨(dú)特作用。我國育種家利用空間環(huán)境誘變水稻種子，培育出了航恢7號(hào)、航香10號(hào)、湘輻259、華航311號(hào)、Y兩優(yōu)1173等多個(gè)抗稻瘟病恢復(fù)系、新資源和新組合[2]。李強(qiáng)等[3]使用N離子束注入連麥6號(hào)，發(fā)現(xiàn)劑量為6×1017N+·cm-2時(shí)，第2代種子突變頻率最高，總突變率達(dá)11.02%，篩選到3個(gè)超親變異性狀(矮桿多蘗、早熟性、大穗)植株，豐富了種質(zhì)資源。李博等[4]采用5種不同劑量7Li離子束注入對大麥進(jìn)行誘變處理，發(fā)現(xiàn)在30 Gy劑量處理下，粒型變異幅度最大，50 Gy劑量下，品質(zhì)變異幅度最大，且30 Gy劑量處理獲得的極端突變體最多，并能獲得較為理想的變異。

高通量測序技術(shù)正向著高通量、低成本、長讀取長度的方向發(fā)展，其中轉(zhuǎn)錄組測序技術(shù)(ribonucleic acid-sequencing, RNA-seq)日漸成熟?；谛乱淮咄繙y序技術(shù)的小RNA(small RNA)測序，可一次性獲得數(shù)百萬條小RNA序列，能夠快速地鑒定某種組織在特定狀態(tài)下的所有已知小RNA并發(fā)現(xiàn)新的小RNA，為小RNA功能研究提供了有力工具。小RNA是生物體內(nèi)一類重要的功能分子，其主要功能是誘導(dǎo)基因沉默，調(diào)控細(xì)胞生長、發(fā)育、基因轉(zhuǎn)錄和翻譯等生物學(xué)過程。小RNA根據(jù)其起源、結(jié)構(gòu)，關(guān)聯(lián)蛋白和生物作用可分為三類：短干擾RNA(short-interferina RNA, siRNA)、微小RNA(micro RNA, miRNA)和PIWI蛋白互作RNA(PIWI-interacting RNA, piRNA)[5-6]。在這些小RNA中，miRNA是內(nèi)源性的20～24 nt長的非編碼實(shí)體，廣泛存在于真核生物中[7-8]。這些miRNA在轉(zhuǎn)錄和轉(zhuǎn)錄后水平上調(diào)控基因表達(dá)[9]，在植物生長發(fā)育中發(fā)揮重要作用[6,10-11]。miRBase(http://mirbase.org/)發(fā)布了24 521 個(gè)miRNA位點(diǎn)，共生產(chǎn)出30 424個(gè)成熟miRNA產(chǎn)物，屬于206個(gè)種屬[12]。

近年來，隨著現(xiàn)代分子生物學(xué)技術(shù)的進(jìn)步，輻射誘變突變體的遺傳及突變機(jī)制研究正在逐步開展，在細(xì)胞水平上主要闡述染色體畸形和突變的關(guān)系，在分子水平上圍繞DNA損傷、修復(fù)及其與突變的關(guān)系開展研究[13]，而損傷效應(yīng)在RNA以及蛋白質(zhì)水平上的研究報(bào)道相對較少。為探討不同劑量下60Co-γ射線輻射后的損傷效應(yīng)在水稻小RNA水平上的差異，本研究采用Illumina HiSeq高通量測序方法對不同劑量輻射后的高粳糯品種三葉期葉片進(jìn)行小RNA轉(zhuǎn)錄組測序及相關(guān)分析，以期從轉(zhuǎn)錄水平為水稻輻射誘變引起的苗期損傷效應(yīng)的分子機(jī)理研究提供參考。

1 材料與方法

1.1 植物材料與誘變方法

水稻高粳糯品種是湖南湘西自治州農(nóng)家種植的糯稻品種，植株過高，成熟后期易倒伏。為獲得矮化植株變異突變體，將其干種子在湖南省農(nóng)業(yè)科學(xué)院核農(nóng)學(xué)與航天育種研究所湖南輻照中心進(jìn)行60Co-γ射線輻射處理，輻射劑量為0(CK)、300(Gy3)、400 Gy(Gy4)，平均劑量率為8.37 Gy·min-1。將輻射后的種子浸種催芽，播種到水稻田土盆(35 cm×25 cm×5 cm)中，每盆6株×6行，共3盆。光溫生長箱中生長(10 h/14 h，光/暗，光照強(qiáng)度10 000 Lx)，定期澆水。取三葉期葉片，20株為1個(gè)重復(fù)，共3個(gè)生物學(xué)重復(fù)，采集后的葉片浸泡于液氮中保存。

1.2 試驗(yàn)方法

1.2.1 總RNA提取和純度檢測、小RNA文庫構(gòu)建、Illumina測序 根據(jù)RNAiso Plus[寶生物工程(大連)有限公司]提取試劑盒的說明，從組織樣本中分離出總RNA，部分步驟稍作修改。在1%瓊脂糖凝膠上監(jiān)測RNA降解和污染，送至北京諾禾致源科技股份有限公司檢測純度。RNA濃度由Qubit 2.0熒光計(jì)(美國Life Technologies公司)測量，完整性由2100安捷倫生物分析儀(美國安捷倫有限公司)評估。純化后的RNA用于構(gòu)建小RNA文庫，并在Illumina HiSeqTM2500/MiSeq測序平臺(tái)(北京諾禾致源科技股份有限公司)進(jìn)行測序。

1.2.2 測序數(shù)據(jù)處理分析 測序得到的原始讀取序列(raw reads)去除含有帶接頭的、低質(zhì)量的讀取序列(reads)，得到過濾序列(clean reads)。篩選出序列長度為18～40 nt的小RNA用于進(jìn)一步分析。使用Bowtie軟件[14]將篩選后的小RNA定位到參考序列上，分析小RNA在參考序列上的分布。將上述比對到參考序列上的reads，用miRBase20.0搜索已知的miRNA[15]。整合miREvo[16]和mirdeep2[17]這些miRNA預(yù)測軟件來進(jìn)行新miRNA的分析；用水稻非編碼RNA(non-coding RNA，ncRNA)序列注釋測得的小RNA分析ncRNA，用水稻重復(fù)序列信息注釋測得小RNA分析重復(fù)序列，用PlantNATsDB(http://bis.zju.edu.cn/pnatdb/)檢測天然反義轉(zhuǎn)錄段干擾RNA(Natural Antisense Transcripts-siRNA，NAT-siRNA.)，反式短干擾RNA(Trans Acting Short-interfering RNA，TAS-iRNA)分析采用水稻和擬南芥已知TAS基因?yàn)閿?shù)據(jù)庫進(jìn)行同源比對，運(yùn)用順式元件短干擾RNA(ta-siRNA，trans-acting siRNA)識(shí)別軟件UEA sRNA tools[18]來進(jìn)行TAS預(yù)測；對各樣本中已知和新miRNA進(jìn)行表達(dá)量統(tǒng)計(jì)，并用TPM[19]進(jìn)行表達(dá)量歸一化處理；采用基于負(fù)二項(xiàng)分布的DESeq2[20]進(jìn)行差異表達(dá)的miRNAs分析。用火山圖可以推斷差異miRNA的整體分布情況，從差異倍數(shù)(fold change)和校正后的顯著水平(Padj，q值)兩個(gè)水平進(jìn)行評估。默認(rèn)差異miRNA的篩選條件為：Padj<0.05；將psRobot.tar運(yùn)用到psRobot軟件[21]中預(yù)測miRNAs的靶基因。

1.2.3 GO富集和KEGG通路分析 GO(Gene Ontology)富集分析方法為GOseq[22]；KEGG[23](Kyoto Encyclopedia of Genes and Genomes)是有關(guān)Pathway的主要公共數(shù)據(jù)庫，以候選靶基因KEGG富集散點(diǎn)圖表示KEGG富集分析結(jié)果，通過Rich factor、q值和富集到此通路上的基因個(gè)數(shù)來衡量，Rich factor越大，表示富集的程度越大。q值是做過多重假設(shè)檢驗(yàn)校正之后的P值，q值的取值范圍為[0，1]，越接近于零，表示富集越顯著。

2 結(jié)果與分析

2.1 測序數(shù)據(jù)總結(jié)

Gy3、Gy4及對照組樣品clean data分別達(dá)到98.42%、94.98%以及96.67%，堿基GC含量均在50%以上，Q30堿基百分比均在97.3%及以上，說明測序質(zhì)量良好，可進(jìn)行后續(xù)分析(表1)。從clean data中選取長度在18～30 nt的小RNA序列作進(jìn)一步研究，對照組小RNA的表達(dá)頻率分布前四位分別為18.98%(24 nt)、9.99%(29 nt)、9.71%(30 nt)、8.73%(28 nt)；Gy3處理組中表達(dá)頻率分布前四位分別為17.94%(24 nt)、9.33%(20 nt)、8.90%(25 nt)、7.91%(23 nt)；Gy4處理組中表達(dá)頻率分布前四位分別為19.55%(24 nt)、9.39%(20 nt)、8.19%(29 nt)、7.78%(21 nt)(圖1)，24 nt的小RNA明顯高于其他長度的小RNA，經(jīng)γ射線處理后高粳糯20-25 nt小RNA所占比例有所增加。

圖1 不同劑量下高粳糯的小RNA長度分布Fig.1 Length distribution of small RNA under different doses treated Gaogengnuo

表1 測序數(shù)據(jù)小結(jié)Table 1 Summary of resequencing data

一般情況下，質(zhì)量較好的植物樣品中rRNA總量所占比例應(yīng)低于60%。由表2可知，3個(gè)處理rRNA總量均遠(yuǎn)低于60%。天然反義轉(zhuǎn)錄物(natural antisense transcripts，NAT)是指可以跟其他轉(zhuǎn)錄本互補(bǔ)形成RNA雙鏈的編碼或非編碼RNA序列，對照組中NAT所占比例為17.20%，Gy3和Gy4組NAT占比分別為20.66%和19.84%，兩種劑量輻射處理后NAT比例都有所增加。對照組新的miRNAs有863個(gè)，Gy3有1 203個(gè)，Gy4有867個(gè)，300 Gy處理下較對照增加了340個(gè)新miRNAs。未被注釋序列占比隨著輻射劑量的增加而增加。外顯子中的小RNA多位于正鏈，內(nèi)含子的小RNA則多位于負(fù)鏈。

表2 小RNA種類和數(shù)量匯總Table 2 Summary of small RNA types and amount

2.2 已知和新miRNA鑒定

對小RNA文庫進(jìn)行搜索來識(shí)別已知的miRNA，經(jīng)過BLASTN搜索和序列分析，Gy3處理中共有312個(gè)已知前體，265個(gè)miRNA成熟序列；Gy4中共有315個(gè)已知前體，261個(gè)miRNA成熟序列。miRNA前體，能夠用來預(yù)測新的miRNA，在對照中檢測到31個(gè)已知前體，26個(gè)miRNA成熟序列，Gy3中共有34個(gè)已知前體，29個(gè)miRNA成熟序列，Gy4中共有34個(gè)已知前體，29個(gè)miRNA成熟序列(表3)。其中novel_209、novel_279、novel_91這3個(gè)新miRNA只在輻射處理之后的高粳糯中被檢測到，novel_261在3種處理中均是數(shù)目最多的新miRNA，其次是novel_34。與對照相比，Gy3處理使得相同類型新miRNA數(shù)量增多(表4)。

表3 已知和新預(yù)測的miRNATable 3 Known and novel miRNA

表4 新miRNA序列和數(shù)目Table 4 Sequence and number of novel miRNA

2.3 新預(yù)測miRNA的首位堿基偏好性

miRNA在由前體發(fā)育為成熟體時(shí)，其過程是由Dicer酶切完成的，酶切位點(diǎn)的特異性使得miRNA成熟體序列首位堿基具有很強(qiáng)的偏向性。通過分析不同長度miRNA的首位點(diǎn)堿基分布，發(fā)現(xiàn)在新miRNA中3種處理組24 nt的miRNA占絕大多數(shù)，其首位堿基位點(diǎn)是A，18 nt和19 nt長度的miRNA中100%都是5′U為首位酶切位點(diǎn)。對照組中24～27 nt長度的新miRNA以5′A首位酶切位點(diǎn)占主要比例；Gy3處理組中20 nt、28～29 nt與18～19 nt長度的新miRNA的首位堿基酶切位點(diǎn)相同，100%都是U，在30 nt時(shí)則都是G；Gy4處理組中28 nt與18～19 nt長度的新miRNA的首位堿基酶切位點(diǎn)相同，100%都是U，在29 nt時(shí)則都是A(圖2)。新miRNA中24 nt的miRNA在Gy3處理組中較對照增加了290個(gè)，在Gy4處理組中較對照減少了81個(gè)。

圖2 新預(yù)測miRNA的首位堿基偏好性Fig.2 The first nucleotide bias of novel miRNA

2.4 差異miRNA分析鑒定

兩兩比較3個(gè)處理組文庫中miRNA的標(biāo)準(zhǔn)化表達(dá)水平，以識(shí)別差異表達(dá)miRNA。差異miRNA篩選條件為：q值<0.01，| log2(foldchange)|>1。維恩圖分析結(jié)果顯示Gy3 vs CK比較組有68個(gè)差異miRNA，包括61個(gè)已知miRNA和7個(gè)新的miRNA(novel_149、novel_205、novel_232、novel_238、novel_288、novel_43、novel_95)； Gy4 vs CK比較組有65個(gè)差異miRNA，包括60個(gè)已知miRNA和5個(gè)新的miRNA(novel_100、novel_102、novel_150、novel_261、novel_63)；Gy3 vs CK與Gy4 vs CK再進(jìn)行差異比較，有86個(gè)差異miRNA，包括79個(gè)已知miRNA和7個(gè)新miRNA(novel_184、novel_195、novel_209、novel_240、novel_279、novel_90、novel_91)(圖3-a)。火山圖顯示Gy3 vs CK比較組有154個(gè)表達(dá)差異，其中Gy3中有88個(gè)上調(diào)，66個(gè)下調(diào)；Gy4 vs CK比較組有151個(gè)表達(dá)差異，其中Gy4中有75個(gè)上調(diào)，76個(gè)下調(diào)；Gy3 vs Gy4再進(jìn)行差異比較，有86個(gè)上調(diào)，68個(gè)下調(diào)表達(dá)(圖3-b、c、d)。

注：a：Gy3 vs CK和Gy4 vs CK差異miRNA維恩圖，每個(gè)圈里的數(shù)字代表相應(yīng)樣本間的差異miRNAs，重疊部分代表共同的差異miRNAs；b、c、d代表Gy3 vs CK、Gy4 vs CK和Gy3 vs Gy4的差異miRNA火山圖，正三角和倒三角分別代不同比較組中上調(diào)和下調(diào)的miRNAs，黑點(diǎn)代表無顯著差異。Note: a:The Veen diagram of differentially expressed miRNAs of Gy3 vs CK, Gy4 vs CK. The number in each circlerepresented the expressed miRNAs in the corresponding sample, and in the overlapping part of the circles represented the co-expressed miRNAs. b, c, d: Dare volcanic diagram of differentially expressed miRNA in Gy3 vs CK, Gy4 vs CK and Gy3 vs Gy4. The regular triangle and inverted triangle represented significantly up-regulated and down-regulated miRNAs, the black dots represented no significant difference.圖3 3個(gè)處理組比較差異miRNAsFig.3 The differentially expressed miRNAs in three treatments

2.5 差異miRNA候選靶基因GO分析

將CK、Gy3和Gy4分別兩兩比較，得到差異表達(dá)miRNA，根據(jù)miRNA與其靶基因間的對應(yīng)關(guān)系，對每組差異表達(dá)miRNA的靶基因的集合進(jìn)行GO分析。Gy3和CK相比較沒有顯著差異的miRNA，無需進(jìn)行GO分析。Gy4和CK的差異miRNA的靶基因進(jìn)行GO分析，如圖4所示，差異顯著的靶基因均涉及生物過程、細(xì)胞組分、分子功能三大類，但分子功能所占基因數(shù)目最多。其中生物過程中共有20個(gè)小類型，前三類分別是羧酸代謝過程(51)、酮酸代謝過程(51)和有機(jī)酸代謝過程(51)。連續(xù)的分析表明，靶基因在功能上高度多樣化。在細(xì)胞組件方面，將基因分成18個(gè)小類，其中前三類分別是細(xì)胞質(zhì)(96)、細(xì)胞器部分(75)和細(xì)胞內(nèi)細(xì)胞器部分(74)。在分子功能方面共有20個(gè)小類型，前三類分別是分子功能(619)、離子結(jié)合(253)和蛋白質(zhì)結(jié)合(184)?？偨Y(jié)三種類型，分子功能的基因最多，其次是結(jié)合和代謝，可能這些差異miRNA與代謝密切相關(guān)。Gy3和Gy4兩種劑量處理組比較差異miRNA所涉及的靶基因富集如圖5所示，只富集到生物過程和分子功能，其中生物過程又分為2個(gè)小類共27個(gè)基因，DNA組裝(13)和DNA構(gòu)象變化(14)。分子功能分為5個(gè)小類共527個(gè)基因，離子結(jié)合(252)、過渡金屬離子結(jié)合(137)、氧化還原酶活性(82)、電子載體活性(25)、銅離子結(jié)合(31)?？赡苓@些差異miRNA涉及生物體內(nèi)氧化還原反應(yīng)。

注：1：氨?；z氨酰tRNA；2：蛋白質(zhì)翻譯的tRNA氨?；?；3：氨基酸活化作用；4：tRNA氨?；?；5：細(xì)胞氨基酸代謝過程；6：tRNA代謝過程；7：DNA組裝；8：羧酸代謝過程；9：酮酸代謝過程；10：有機(jī)酸代謝過程；11：ncRNA代謝過程；12：DNA構(gòu)象變化；13：正調(diào)控分解代謝過程；14：正調(diào)控自噬；15：細(xì)胞分解代謝過程正調(diào)控；16：細(xì)胞分解代謝過程調(diào)控；17：自噬調(diào)控；18：分解代謝過程調(diào)控；19：發(fā)病機(jī)理；20：自噬；21：中間絲狀體；22：中間纖維細(xì)胞骨架；23：膜的外在成分；24：部分細(xì)胞骨架；25：細(xì)胞骨架；26：轉(zhuǎn)運(yùn)質(zhì)子的v型ATP酶，V0結(jié)構(gòu)域；27：轉(zhuǎn)運(yùn)質(zhì)子的v型ATP酶復(fù)合物；28：細(xì)胞質(zhì)；29：細(xì)胞內(nèi)無膜邊界細(xì)胞器；30：無膜邊界細(xì)胞器；31：基膜；32：細(xì)胞外基質(zhì)成分；33：細(xì)胞器部分；34：細(xì)胞內(nèi)細(xì)胞器部分；35：蛋白質(zhì)的細(xì)胞外基質(zhì)；36：二扇形ATP酶復(fù)合物，質(zhì)子轉(zhuǎn)運(yùn)域；37：內(nèi)膜系統(tǒng)；38：細(xì)胞外基質(zhì)；39：動(dòng)力蛋白結(jié)合；40：蛋白質(zhì)同源二聚化活性；41：同源蛋白結(jié)合；42：絲氨酸t(yī)RNA連接酶活性；43：轉(zhuǎn)錄因子結(jié)合；44：氨酰連接酶活性；45：連接酶活性，形成碳氧鍵；46：連接酶活性，形成氨酰tRNA及相關(guān)化合物；47：蛋白質(zhì)二聚化活性；48：未折疊蛋白結(jié)合；49：連接酶活性；50：銅離子結(jié)合；51：電子載體活性；52：序列特異性DNA結(jié)合；53：核酸結(jié)合轉(zhuǎn)錄因子活性；54：轉(zhuǎn)錄因子活性，序列特異性DNA結(jié)合；55：分子功能；56：細(xì) 胞骨架蛋白結(jié)合；57：離子結(jié)合；58：蛋白結(jié)合。BP：生物過程；CC：細(xì)胞組分；MF：分子功能。Note: 1: Seryl-tRNA aminoacylation. 2: tRNA aminoacylation for protein translation. 3: Amino acid activation. 4: tRNA aminoacylation. 5: Cellular amino acid metabolic process. 6: tRNA metabolic process. 7: DNA packaging. 8: Carboxylic acid metabolic process. 9: Oxoacid metabolic process. 10: Organic acid metabolic process. 11: ncRNA metabolic process. 12: DNA conformation change. 13: Positive regulation of catabolic process. 14: Positive regulation of autophagy. 15: Positive regulation of cellular catabolic process. 16: Regulation of cellular catabolic process. 17: Regulation of autophagy. 18: Regulation of catabolic process. 19: Pathogenesis. 20: Autophagy. 21: Intermediate filament. 22: Intermediate filament cytoskeleton. 23: Extrinsic component of membrane. 24: Cytoskeletal part. 25: Cytoskeleton. 26: Proton-transporting Ⅴ-type ATPase, V0 domain. 27: Proton-transporting Ⅴ-type ATPase complex. 28: Cytoplasm. 29: Intracellular non-membrane-bounded organelle. 30: Non-membrane-bounded organelle. 31: Basement membrane. 32: Extracellular matrix component. 33: organelle part. 34: Intracellular organelle part. 35: Proteinaceous extracellular matrix. 36: Proton-transporting two-sector ATPase complex, proton-transporting domain. 37: Endomembrane system. 38: Extracellular matrix. 39: Dynein binding. 40: Protein homodimerization activity. 41: Identical protein binding. 42: Serine-tRNA ligase activity. 43: Transcription factor binding. 44: Aminoacyl-tRNA ligase activity. 45: Ligase activity, forming carbon-oxygen bonds. 46: Ligase activity, forming aminoacyl-tRNA and related compounds. 47: Protein dimerization activity. 48: Unfolded protein binding. 49: Ligase activity. 50: copper ion binding. 51: Electron carrier activity. 52: Sequence-specific DNA binding. 53: Nucleic acid binding transcription factor activity. 54: Transcription factor activity, sequence-specific DNA binding. 55: Molecular function. 56: Cytoskeletal protein binding. 57: Ion binding. 58: Protein binding. BP: Biological Process. CC: CellularComponent. MF: Molecular Function.圖4 Gy4 vs CK的GO富集圖Fig.4 Gene ontology terms for differentially expressed miRNAs in Gy4 vs CK

注：1：DNA組裝；2：DNA構(gòu)象變化；3：銅離子結(jié)合；4：電子載體活性；5：氧化還原酶活性；6：過渡離子結(jié)合；7：離子結(jié)合。Note: 1: DNA packaging. 2: DNA conformation change. 3: Copper ion binding. 4: Electron carrier activity. 5: Oxidoreductase activity. 6: Transition metal ion binding. 7: Ion binding.圖5 Gy3 vs Gy4差異miRNA的GO富集圖Fig.5 Gene ontology terms for differentially expressed miRNAs in Gy3 vs Gy4

2.6 差異miRNA候選靶基因KEGG富集分析

生物體內(nèi)不同基因相互協(xié)調(diào)行使其生物學(xué)功能，通過途徑顯著性富集能確定候選靶基因參與的最主要生化代謝途徑和信號(hào)轉(zhuǎn)導(dǎo)途徑。將Gy3 vs CK、Gy4 vs CK、 Gy3 vs Gy4 的差異miRNA對應(yīng)的靶基因進(jìn)行KEGG富集分析，結(jié)果以候選靶基因KEGG富集散點(diǎn)圖方式展示(圖6)，Gy3 vs CK組富集到的前20條途徑中植物病原互作在基因數(shù)量方面富集程度最高(8個(gè)基因)，其次是過氧物酶體(7個(gè)基因)和內(nèi)吞作用(7個(gè)基因)；富集最顯著的則是植物晝夜節(jié)律，q值為0.49，有5個(gè)候選基因。Gy4 vs CK組富集到的前20條途徑中代謝途徑在基因數(shù)量方面富集程度最高(57個(gè)基因)，其次是甘油磷脂代謝(8個(gè)基因)、精氨酸和脯氨酸代謝(6個(gè)基因)、真核生物的核糖體發(fā)生(6個(gè)基因)和內(nèi)吞作用(6個(gè)基因)；富集最顯著的則是甘油磷脂代謝和醚脂類代謝，q值均為0.12，其中醚脂類代謝有5個(gè)候選基因。Gy3 vs Gy4組富集到的前20條途徑中甘油磷脂代謝、真核生物的核糖體發(fā)生和內(nèi)吞作用在基因數(shù)量方面富集程度最高(7個(gè)基因)，其次是醚脂類代謝(6個(gè)基因)；富集最顯著的是醚脂類代謝，q值為0.02，其次是甘油磷脂代謝和氨腈酸代謝，q值為0.35。

注：q值的顏色表示富集因子的顯著性，圓圈表示所涉及的靶基因，大小為與基因數(shù)量成正比。Note: The coloring of the q-values indicates thesignificance of the rich factor; the circle indicates the target genes that are involved, and the size isproportional to the gene number.圖6 差異miRNA的KEGG分類圖Fig.6 KEGG pathway of the differentially expressed miRNAs

3 討論

目前，轉(zhuǎn)錄組測序和從頭組裝是獲取高通量轉(zhuǎn)錄本序列最快速有效的方法，并成功應(yīng)用于新基因發(fā)現(xiàn)、功能基因定位和分子標(biāo)記開發(fā)[24-25]。水稻常規(guī)品種高粳糯株高約150 cm，為創(chuàng)制矮化株型的突變體，對高粳糯干種子進(jìn)行不同劑量γ射線輻射處理，輻射后在M1代的苗期都出現(xiàn)了不同程度的損傷效應(yīng)。目前基于植物輻射誘變引起當(dāng)代苗期損傷效應(yīng)的分子機(jī)理研究報(bào)道相對較少，本研究在RNA轉(zhuǎn)錄水平上探討和分析誘變引起的損傷效應(yīng)，以期為輻射誘變育種技術(shù)提供理論支撐。

基于高通量測序技術(shù)的小RNA測序已經(jīng)在許多物種中展開，Xu等[26]對嗜水氣單胞菌敏感(susceptible grass carp, SGC)和耐水草魚(resistant grass carp, RGC)進(jìn)行sRNA測序，共鑒定出61個(gè)保守miRNA和124個(gè)候選新miRNA以及21個(gè)差異miRNA。He等[27]分別構(gòu)建了白肉甘薯(xusshu-18)和紫肉甘薯(xuzshu-3)小RNA和降解組文庫。小RNA測序共鑒定出191個(gè)已知miRNAs和33個(gè)新miRNA，降解組測序鑒定出115個(gè)已知i-miRNAs和5個(gè)新i-miRNAs切割的180個(gè)靶基因，有121個(gè)差異miRNA，綜合分析發(fā)現(xiàn)有26個(gè)差異表達(dá)的miRNA和36個(gè)相應(yīng)的靶點(diǎn)可能參與花青素的生物合成。水稻根系鹽脅迫響應(yīng)miRNA和tRF的鑒定：孟淑君等[28]共檢測到12種水稻根系中鹽脅迫響應(yīng)miRNA，靶基因多為轉(zhuǎn)錄因子編碼基因，推測其通過對其靶基因轉(zhuǎn)錄因子的轉(zhuǎn)錄后調(diào)控參與了鹽脅迫響應(yīng)的表達(dá)調(diào)控。本研究中Gy3 處理共鑒定出265個(gè)已知miRNA成熟序列，29個(gè)新的miRNA成熟序列；Gy4處理共鑒定出261個(gè)已知miRNA成熟序列，29個(gè)新miRNA成熟序列，兩種劑量處理后對已知miRNA和新的miRNA的貢獻(xiàn)均等，與對照組相比均略增加。推測300 Gy 輻射劑量更適合發(fā)掘新的miRNA，400 Gy輻射劑量更適合發(fā)掘已知miRNA。Gy3 vs CK組有154個(gè)差異miRNA，Gy4 vs CK組有151個(gè)差異miRNA，在這些差異表達(dá)基因中，有很多miRNA家族在其他物種中也被報(bào)道如oa-miR397和osa-miR168。如紅花miR397a的表達(dá)可增強(qiáng)植物對NaCl的敏感性[29]；mir168b突變體表現(xiàn)出脫落酸低敏性和干旱超敏性[30]。輻射后的高粳糯可能對某些非生物脅迫表現(xiàn)出一定的敏感性。

γ射線輻射生物體后，生物體產(chǎn)生過量自由基不能及時(shí)清除，從而引起輻射損傷[31]。何穎等[32]發(fā)現(xiàn)低劑量γ射線輻射后的大鼠血清代謝產(chǎn)物發(fā)生變化，主要涉及免疫功能、能量代謝、糖代謝和脂代謝。水稻在經(jīng)γ射線輻射后在當(dāng)代苗期出現(xiàn)不同程度損傷效應(yīng)，生長后期損傷效應(yīng)逐漸恢復(fù)，可能涉及某些基因通過調(diào)節(jié)表達(dá)量來迅速應(yīng)對外界脅迫效應(yīng)，采用生物信息學(xué)對輻射后的高粳糯差異基因涉及途徑分析，Gy4和CK的差異miRNA靶基因富集到的分子功能的基因最多，其次是結(jié)合和代謝，可能這些差異miRNA與代謝密切相關(guān)。通過KEGG富集分析，Gy4 vs CK組代謝途徑富集程度最高，富集最顯著的是甘油磷脂代謝和醚脂類代謝。Gy3 vs Gy4組甘油磷脂代謝、真核生物的核糖體發(fā)生和內(nèi)吞作用在基因數(shù)量方面富集程度最高，富集最顯著的是醚脂類代謝，可見輻射處理相對更多影響高粳糯參與代謝活動(dòng)的基因。60Co-γ輻射可能通過影響水稻幼苗甘油磷酸酯代謝和醚脂類代謝抑制了其生長發(fā)育。下一步工作將結(jié)合相關(guān)突變體表型開展候選基因的表達(dá)量驗(yàn)證，以期克隆發(fā)現(xiàn)相關(guān)基因。

4 結(jié)論

利用高通量測序平臺(tái)對不同劑量γ射線輻射后的水稻品種高粳糯進(jìn)行小RNA測序。300 Gy和400 Gy劑量處理后對已知miRNA和新miRNA的貢獻(xiàn)均等，與對照組相比均有少量增加。GO聚類分析結(jié)果顯示Gy4 vs CK的差異靶基因在分子功能的基因最多。Gy3 vs Gy4組差異靶基因以離子結(jié)合占主要比例。KEGG顯示Gy3 vs CK組集中在植物病原互作和植物晝夜節(jié)律途徑；Gy4 vs CK組集中在代謝途徑；Gy3 vs Gy4組集中在甘油磷脂代謝和醚脂類代謝途徑。本研究從RNA水平變異的角度分析輻射誘變對水稻的影響，以期為進(jìn)一步研究輻射誘變引起水稻當(dāng)代苗期的損傷效應(yīng)提供一種新途徑。