楊興云，喬丹丹，張雅潔，王少青，任俊才，李明陽(yáng)，屈明好，尚盼盼，楊成，黃琳凱，曾兵，*

（1. 西南大學(xué)動(dòng)物科學(xué)技術(shù)學(xué)院，重慶402460；2. 草食動(dòng)物科學(xué)重慶市重點(diǎn)實(shí)驗(yàn)室，重慶 400715；3. 福建傲農(nóng)生物科技集團(tuán)股份有限公司，福建 廈門 361111；4. 四川農(nóng)業(yè)大學(xué)草業(yè)科技學(xué)院，四川 成都 611130）

鴨茅（Dactylis glomerata）又名果園草（orchard grass）、雞腳草（cocksfoot），屬禾本科（Poaceae）早熟禾亞科（Festucoideae）鴨茅屬（Dactylis）［1］。鴨茅是一種高度異花授粉的冷季型多年生禾本科牧草，自然分布于歐洲及我國(guó)，其作為飼料作物有極高的營(yíng)養(yǎng)價(jià)值以及較強(qiáng)的適應(yīng)性和耐陰性［2?3］，但對(duì)水淹脅迫較敏感。在國(guó)內(nèi)鴨茅主產(chǎn)區(qū)，由于洪澇災(zāi)害頻繁，嚴(yán)重影響鴨茅的生長(zhǎng)發(fā)育、品質(zhì)和產(chǎn)量。但目前國(guó)內(nèi)外關(guān)于鴨茅抗逆性研究主要集中在干旱脅迫［4］、耐熱脅迫［5］、銹病脅迫［6］等方面，因此探究鴨茅在水淹脅迫下的抗逆分子機(jī)制有重要意義。

澇害作為植物正常生長(zhǎng)發(fā)育的非生物脅迫因子，對(duì)植物產(chǎn)生的危害主要是阻礙植物根系的供氧，抑制根系呼吸，導(dǎo)致植物重要代謝過程的根細(xì)胞能量供應(yīng)不足［7］，最終植物枯萎甚至壞死。水淹脅迫對(duì)植物生理生化的影響主要包括光合作用、根部抗氧化系統(tǒng)和根部激素合成三方面［8］。長(zhǎng)期水淹后，光合作用產(chǎn)生的光合產(chǎn)物不能通過韌皮部向下運(yùn)輸，導(dǎo)致植物的代謝循環(huán)受阻，葉片凈光合速率、光合色素、氣孔導(dǎo)度降低［9?10］，從而抑制植物正常生長(zhǎng)。此外，根系抗氧化系統(tǒng)也會(huì)發(fā)生相應(yīng)變化，水淹脅迫后植株體內(nèi)的抗壞血酸過氧化物酶（ascorbic acid peroxidase，APX）、谷胱甘肽還原酶（glutathione reductase，GR）、超氧化物歧化酶（superoxide dismutase，SOD）及過氧化物酶（peroxidase，POD）等含量急劇變化［11］。同時(shí)，水淹脅迫還影響植物體內(nèi)激素合成與運(yùn)輸動(dòng)態(tài)平衡，生 長(zhǎng) 素（indole-3-acetic acid，IAA）、細(xì) 胞 分 裂 素（cytokinin，CTK）、脫 落 酸（abscisic acid，ABA）、赤 霉 素（gibberellin，GA）、乙烯（ethylene，ETH）等激素在不同部位增加。其中，乙烯作為衡量植物耐澇性重要的指標(biāo)之一，在有氧部位積累，造成葉片老化脫落，抑制正常植物生長(zhǎng)［12］。

小RNA（microRNAs，miRNA）是一類內(nèi)源性的大小為21～22 個(gè)核苷酸的小分子單鏈RNA，位于基因組非編碼區(qū)［13］。通過靶向mRNA 的切割或抑制翻譯過程參與基因的表達(dá)調(diào)控［14］。miRNA 作為一種重要的調(diào)控因子，通過在轉(zhuǎn)錄后水平與其靶基因mRNA 的完美配對(duì)參與植物對(duì)環(huán)境脅迫的反應(yīng)［15］。研究表明，miRNA 在植物抗逆脅迫反應(yīng)中發(fā)揮著重要作用［16］，例如鹽脅迫［17］、溫度脅迫［18］、營(yíng)養(yǎng)脅迫［19］和水分脅迫［20］等。植物遭遇水淹時(shí)，氣體在水中的擴(kuò)散速度降低，植物組織會(huì)出現(xiàn)嚴(yán)重的缺氧表現(xiàn)，已有相關(guān)研究揭示低氧環(huán)境下miRNA 如何參與相關(guān)基因表達(dá)調(diào)控［21?23］。Licausi 等［24］將擬南芥（Arabidopsis thaliana）根進(jìn)行不同程度的缺氧處理，通過對(duì)1900個(gè)轉(zhuǎn)錄因子（transcription factors，TF）和180 個(gè)miRNA 差異表達(dá)分析后發(fā)現(xiàn)，低氧誘導(dǎo)基因的調(diào)控是通過各種TF 組合的同時(shí)相互作用來控制的，miRNA 在植物低氧基因調(diào)控的早期階段起次要作用［24］。劉智捷［25］對(duì)3 種不同耐澇表型的玉米（Zea mays）幼苗進(jìn)行短期漬水脅迫處理，分析各處理?xiàng)l件下根系miRNA 與其靶基因的表達(dá)動(dòng)態(tài)，結(jié)果表明miR159、miR164、miR167、miR393、miR408和miR528基因主要參與根系發(fā)育和水淹脅迫響應(yīng)。進(jìn)一步分析發(fā)現(xiàn)，約有1/2 目的基因與其miRNA 表達(dá)量呈負(fù)相關(guān)，其中大部分是參與低氧脅迫反應(yīng)的轉(zhuǎn)錄因子。在其他物種中，Jin 等［26］對(duì)水淹脅迫下的蓮花（Nelumbo nucifera）進(jìn)行高通量測(cè)序研究耐澇相關(guān)的miRNA，發(fā)現(xiàn)NNU_Far-miR159、NNU_GMA-miR393h 和NNU_Aly-miR319c-3p 在蓮花對(duì)淹水的響應(yīng)中起著重要的調(diào)節(jié)作用。

近年來，隨著全球自然環(huán)境的持續(xù)發(fā)展惡化以及頻繁的極端天氣，導(dǎo)致某些地區(qū)暴雨頻發(fā)，洪水災(zāi)害越來越嚴(yán)重［27?28］。中國(guó)也是遭受洪澇災(zāi)害較為嚴(yán)重的國(guó)家，特別是長(zhǎng)江中下游地區(qū)，受災(zāi)面積占全國(guó)總受災(zāi)面積的3/4［29］。2020年8月，四川、重慶等地嚴(yán)重水災(zāi)造成了巨大社會(huì)經(jīng)濟(jì)損失。此外，持續(xù)的氣候變暖預(yù)示著世界所有區(qū)域的洪澇事件的頻率和強(qiáng)度還會(huì)增加［30］，這將導(dǎo)致更大的經(jīng)濟(jì)損失。鴨茅作為優(yōu)質(zhì)牧草，深度解析植物耐澇的分子響應(yīng)機(jī)制對(duì)加快作物耐澇品種選育具有重要意義。雖有人對(duì)水淹脅迫后“斯巴達(dá)”鴨茅進(jìn)行了轉(zhuǎn)錄組測(cè)序研究［31］，但對(duì)于鴨茅在水淹脅迫下miRNA 表達(dá)特征的研究尚未見報(bào)道。前期用人工模擬澇害脅迫，通過形態(tài)學(xué)、生理學(xué)根部顯微結(jié)構(gòu)以及組學(xué)研究［32?33］，篩選到了耐澇性較強(qiáng)的“滇北”鴨茅。

鑒于此，本研究以耐澇性較強(qiáng)的“滇北”鴨茅為試驗(yàn)材料，利用高通量測(cè)序與生物信息學(xué)技術(shù)，篩選與鴨茅水淹抗性相關(guān)的miRNA，分析它們?cè)谒兔{迫中的差異表達(dá)，旨在為今后研究鴨茅抗?jié)匙饔脵C(jī)制及抗?jié)尺z傳改良工作提供參考，同時(shí)為耐澇草本植物新品種的選育提供借鑒。

1 材料與方法

1.1 試驗(yàn)材料與處理

以鴨茅品種“滇北”（簡(jiǎn)稱“DB”）為試驗(yàn)材料，種子由四川農(nóng)業(yè)大學(xué)草業(yè)科技學(xué)院提供。試驗(yàn)于2018年9月開始，篩選均勻飽滿的鴨茅種子播種于育苗盤中，挑選長(zhǎng)勢(shì)大體一致的幼苗移栽至花盆（口徑16 cm，高14 cm），再放入人工培養(yǎng)箱，溫度設(shè)置為22 ℃/15 ℃（晝/夜），光周期14 h/10 h（晝/夜），光照10000 lx，土壤相對(duì)含水量保持在50%左右。根據(jù)植株生長(zhǎng)需要每周澆2 次1/2 Hoagland 溶液，待植株長(zhǎng)至5 cm 及葉片覆蓋率超過80%時(shí)進(jìn)行試驗(yàn)。

試驗(yàn)采用盆栽水淹脅迫，將3 盆長(zhǎng)勢(shì)良好且一致的鴨茅放入水箱（長(zhǎng)80 cm×寬57 cm×高50 cm）中，然后加水淹沒至植株葉片最高處，并始終保持水位淹沒鴨茅頂部。另外3 盆不水淹，做對(duì)照處理。于水淹脅迫處理0 h取樣對(duì)照組（CK），水淹脅迫處理8 和24 h 取樣處理組（D），取樣3 次，每次取樣3 個(gè)重復(fù)，共9 個(gè)樣品。采用混樣法剪取鴨茅葉片，液氮迅速冷凍，后置于?80 ℃低溫冰箱保存?zhèn)溆?。樣品送至武漢菲沙基因公司進(jìn)行測(cè)序。

1.2 小RNA 文庫(kù)的構(gòu)建與高通量測(cè)序

利用Trizol 試劑盒（Invitrogen，USA）分別提取9 個(gè)鴨茅葉片樣品中的總RNA。總RNA 質(zhì)量檢測(cè)：1）通過進(jìn)行瓊脂糖凝膠電泳分析RNA 降解和污染的程度；2）采用Nanodrop 測(cè)試RNA 濃度（OD260/280）；3）利用Qubit 對(duì)RNA 濃度進(jìn)行精確定量；4）采用Agilent 2100 檢測(cè)RNA 完整性。最后僅選擇質(zhì)量檢測(cè)合格的總RNA 樣品進(jìn)行小RNA（small RNA，sRNA）文庫(kù)構(gòu)建及測(cè)序工作（Illumina Hiseq TM2100）。

1.3 測(cè)序原始數(shù)據(jù)預(yù)處理

對(duì)高通量測(cè)序獲得的sRNA 進(jìn)行生物信息學(xué)分析前，需先對(duì)原始數(shù)據(jù)進(jìn)行雜質(zhì)處理。數(shù)據(jù)過濾的步驟如下：1）去除低質(zhì)量的reads，當(dāng)測(cè)序reads 中含有的低質(zhì)量（≤5）堿基數(shù)超過該條reads 長(zhǎng)度比例的50%時(shí)，過濾掉該reads；2）去除含N（任一個(gè)堿基）比例較高的reads，當(dāng)測(cè)序reads 中含有的N 的含量超過該條reads 長(zhǎng)度比例的10%時(shí)，過濾掉該reads；3）去除有5′接頭污染的reads；4）去除沒有3′接頭序列和插入片段的reads；5）去除3′接頭序列；6）去除含polyA/T/G/C 的reads；7）去除最終長(zhǎng)度小于17 bp 的reads；經(jīng)過上面的步驟過濾數(shù)據(jù)之后，得到了大量的sRNA 序列。

1.4 響應(yīng)水淹脅迫miRNA 的鑒定

將排除核糖體RNA（ribosomalRNA，rRNA）等的sRNA 序列在miRBase［34］數(shù)據(jù)庫(kù)中進(jìn)行比對(duì)分析，綜合成熟體相似度和前體莖環(huán)結(jié)構(gòu)兩方面比對(duì)情況得到鴨茅中miRNA 的表達(dá)情況。具體策略條件如下：1）成熟體相似度方面考慮到鴨茅與其他物種的差異，比對(duì)中允許兩個(gè)堿基錯(cuò)配；2）對(duì)條件1 中成熟體miRNA 進(jìn)行pre-miRNA 的預(yù)測(cè)鑒定，若出現(xiàn)基因組上無pre-miRNA 位點(diǎn)或不能形成二級(jí)莖環(huán)結(jié)構(gòu)的狀況，則排除其為鴨茅保守miRNA，獲得更為準(zhǔn)確和高可行性的驗(yàn)證。

miRNA 前體的標(biāo)志性發(fā)夾結(jié)構(gòu)能夠用來預(yù)測(cè)新的miRNA?；驹硎峭ㄟ^截取一定長(zhǎng)度sRNA 比對(duì)上的參考序列，通過探尋其二級(jí)結(jié)構(gòu)及Dicer 酶切位點(diǎn)信息、能量等特征進(jìn)行分析，預(yù)測(cè)樣品中新miRNA，并進(jìn)行各樣本中匹配上的sRNA 的序列、長(zhǎng)度、出現(xiàn)的次數(shù)等信息，以及不同長(zhǎng)度miRNA 的首位點(diǎn)堿基分布和所有miRNA的各位點(diǎn)堿基分布情況的統(tǒng)計(jì)。對(duì)于各個(gè)數(shù)據(jù)庫(kù)中均未注釋但可比對(duì)到基因組序列的sRNA，將作為候選的新miRNA。采用mirdeep2［35］軟件進(jìn)行新miRNA 的分析，綜合考慮到sRNA 測(cè)序序列表達(dá)豐度，前體miRNA 最低自由能和復(fù)雜的二級(jí)結(jié)構(gòu)，預(yù)測(cè)得到的候選新miRNA 前體的莖環(huán)結(jié)構(gòu)圖，完成新miRNA 的預(yù)測(cè)。

1.5 miRNA 差異表達(dá)分析

基因差異表達(dá)分析的輸入數(shù)據(jù)為miRNA 表達(dá)水平分析中得到的miRNA 對(duì)應(yīng)的read counts 數(shù)據(jù)。用R 語(yǔ)言包DESeq2 進(jìn)行差異分析，篩選閾值為|log2 fold change|＞1&P＜0.05。對(duì)于之前得到的已知miRNA 和新miRNA，做miRNA 靶基因預(yù)測(cè)分析。

由于單個(gè)軟件預(yù)測(cè)的miRNA 往往會(huì)得到特別多的靶基因，因此本研究選用TargetFinder 和psRobot 兩種軟件結(jié)果的交集作為最后的靶基因預(yù)測(cè)分析［36?37］，分別獲取各組差異表達(dá)的miRNA，根據(jù)miRNA 與其靶基因的關(guān)系，對(duì)每組差異表達(dá)miRNA 的靶基因集合分別進(jìn)行GO (gene ontology)功能分析和KEGG (kyoto encyclopedia of genes and genomes)富集分析。

2 結(jié)果與分析

2.1 鴨茅小RNA 測(cè)序數(shù)據(jù)及質(zhì)量控制

應(yīng)用高通量測(cè)序技術(shù)，對(duì)水淹脅迫0、8 和24 h 后的“滇北”鴨茅葉片進(jìn)行小RNA 測(cè)序。測(cè)序結(jié)果由表1 可知，水淹脅迫后“滇北”葉片中獲得109825656 原始測(cè)序數(shù)據(jù)（raw reads）；通過質(zhì)量控制后，樣品中未發(fā)現(xiàn)質(zhì)量值較低序列，經(jīng)過濾獲得高質(zhì)量序列（clean reads）數(shù)量是68789053；可從表中看出所有樣本的GC 含量為51.2%～52.9%（GC 含量：堿基G 和C 的數(shù)量總和占總的堿基數(shù)量的百分比）；Q20≥99.2%，Q30≥98.1%［Q20、Q30：堿基質(zhì)量值Qphred=?10 log10（error rate）數(shù)值大于20、30 的堿基占總體堿基的百分比］。這說明該測(cè)序結(jié)果可信且質(zhì)量較高，可用于進(jìn)一步的分析。

表1 不同鴨茅樣品測(cè)序數(shù)據(jù)統(tǒng)計(jì)Table 1 Statistical results of sequencing data of different D.glomerata samples

2.2 鴨茅小RNA 分類注釋

采用bowtie2［38］將小RNA 比對(duì)到mirBase（v21）數(shù)據(jù)庫(kù)指定范圍的miRNA 序列上，得到小RNA 與miRNA 比對(duì)情況。在水淹脅迫處理0、8 和24 h 的“滇北”鴨茅中鑒定到能比對(duì)上的miRNA 有3259、3309 和3188 個(gè)，分別占能比對(duì)上miRNA 的小RNA 總數(shù)（total sRNA）的4.22%、3.13%、2.59%。

根據(jù)與Rfam 數(shù)據(jù)庫(kù)的比對(duì)作為非編碼RNA（non-coding RNA，ncRNA）的注釋結(jié)果，“滇北”鴨茅非編碼RNA 序列可以注釋為核仁小RNA（small nucleolar RNA，snoRNA）、核糖體RNA（ribosomalRNA，rRNA）、小核RNA（small nuclear RNA，snRNA）、轉(zhuǎn)運(yùn)RNA（transfer ribonucleic acid，tRNA）4 種。結(jié)合基因組上外 顯子（exon），內(nèi)含子（intron）和重復(fù)序列的信息，以及與Rfam，mirBase（v21）數(shù)據(jù)庫(kù)的比對(duì)結(jié)果，比對(duì)上的小RNA 做分類統(tǒng)計(jì)。由表2 可知，比對(duì)到的miRNA 序列從11.00%上升到15.80%，再下降到14.00%，說明水淹脅迫處理8 h 時(shí)，鴨茅可能積極響應(yīng)水淹脅迫。

表2 小RNA 分類統(tǒng)計(jì)Table 2 Small RNA classification statistics（%）

2.3 鴨茅保守miRNA 鑒定及其家族分析

植物中成熟的miRNA 序列高度保守，生物信息學(xué)方法用于比較基因組學(xué)研究從而探索不同物種中miRNA 基 因 家 族的 保 守 性［39］?；?于miRNA 比 對(duì) 分析結(jié)果，對(duì)miRNA 在鴨茅水淹脅迫葉片中的表達(dá)進(jìn)行統(tǒng)計(jì)。通過表3 可知，推測(cè)8 h 是植物遭受水淹脅迫的代謝活躍期。Ren 等［40］對(duì)水淹脅迫處理8 h 后的毛白楊（Populus tomentosa）進(jìn)行小RNA 測(cè)序分析，結(jié)果表明7 個(gè)保守的miRNA 家族和5 個(gè)新的miRNAs 在水淹脅迫下的表達(dá)發(fā)生了顯著變化，對(duì)差異表達(dá)miRNAs 潛在靶點(diǎn)的注釋表明，許多編碼轉(zhuǎn)錄因子、酶和信號(hào)轉(zhuǎn)導(dǎo)成分的基因參與了非生物脅迫反應(yīng)。與水淹脅迫24 h 相比，8 h“滇北”鴨茅有72 個(gè)miRNA 發(fā)生差異表達(dá)，有38 個(gè)上調(diào)，34 個(gè)下調(diào)。由此可推斷，水淹脅迫處理8 h 時(shí)，鴨茅通過調(diào)節(jié)生理代謝活動(dòng)響應(yīng)水淹脅迫。對(duì)檢測(cè)到的已知miRNA 和新miRNA進(jìn)行家族分析，探索其所屬的miRNA 家族在其他物種中的存在情況。如果來源于同一個(gè)家族的miRNA，它們?cè)?～7 這些位點(diǎn)（這幾個(gè)位點(diǎn)是作用于靶基因的種子位點(diǎn)）表現(xiàn)為一致的保守型，就會(huì)被定義為一個(gè)miRNA 家族［41］，其中序列相似性比對(duì)，是最為簡(jiǎn)單直接的方法。對(duì)鑒定到的miRNA 進(jìn)行DESeq2 差異分析，“滇北”鴨茅中總共鑒定到228 個(gè)差異表達(dá)的miRNA?？梢园l(fā)現(xiàn)鴨茅在水淹脅迫下發(fā)生差異表達(dá)的miRNA 主要屬于miR166（78 個(gè)）、miR159（35 個(gè)）、miR167（29 個(gè)）、miR396（17個(gè)）、miR156（21 個(gè)）這5 個(gè)miRNA 家族（表4），表明這幾個(gè)家族的miRNA 在植物響應(yīng)水淹脅迫過程中發(fā)揮了重要作用。

表3 差異基因統(tǒng)計(jì)Table 3 Statistical of differentially expressed genes（DEGs）

表4 鴨茅中響應(yīng)水淹脅迫的差異表達(dá)miRNA 家族Table 4 Differentially expressed miRNA families in D.glomerata in response to flooding stress

2.4 鴨茅miRNA 的靶基因預(yù)測(cè)

miRNA 主要通過與靶mRNA 的結(jié)合，或促使mRNA 降解，或阻礙其翻譯，從而抑制目的基因的表達(dá)。用生物信息學(xué)方法準(zhǔn)確快速地預(yù)測(cè)miRNA 的靶基因，可以為研究miRNA 功能提供線索。對(duì)不同水淹脅迫處理下的鴨茅材料中得到的已知miRNA 和新miRNA 做miRNA 靶基因預(yù)測(cè)分析，獲得了miRNA 與靶基因的對(duì)應(yīng)關(guān)系。利用TargetFinder 和psRobot，進(jìn)一步預(yù)測(cè)了上述miRNA 的靶基因，共獲得9136 個(gè)靶基因，5 個(gè)主要的miRNA 家族中分別鑒定到miR166（1817 個(gè)）、miR167（320 個(gè)）、miR159（835 個(gè)）、miR396（924 個(gè)）、miR156（5204 個(gè)）。再選取5 個(gè)基因家族中差異顯著的miRNA 及其對(duì)應(yīng)的靶基因進(jìn)行分析，通過表5 可知，在研究中發(fā)現(xiàn)某些miRNA 有多個(gè)靶基因與之對(duì)應(yīng)，表明miRNA 在響應(yīng)水淹脅迫過程中并不是單基因調(diào)控，而是一個(gè)miRNA 可以參與調(diào)控多個(gè)生物代謝過程，這也證實(shí)了鴨茅在水淹脅迫中響應(yīng)過程復(fù)雜。

表5 不同的靶基因?qū)?yīng)同一個(gè)miRNATable 5 Different target genes correspond to the same miRNA

2.5 差異miRNA 靶基因GO 富集分析

GO（gene ontology）是基因功能國(guó)際標(biāo)準(zhǔn)分類體系。在預(yù)測(cè)水淹脅迫下響應(yīng)miRNA 的靶基因后，研究這些靶基因在GO 中所映射到terms 的情況并對(duì)其進(jìn)行生物學(xué)富集分析，可以對(duì)水淹響應(yīng)miRNA 靶基因進(jìn)行功能分類集合分析和具有生物學(xué)意義的基因功能描述。GO 為進(jìn)一步解析miRNA 的功能和相關(guān)靶基因研究提供了思路，該數(shù)據(jù)庫(kù)主要分為三大類：生物過程（biological process，BP）、細(xì)胞組成（cellular component，CC）、分子功能（molecular function，MF）［42?43］。

本研究對(duì)水淹脅迫下“滇北”鴨茅0 h vs 8 h 的70 個(gè)差異表達(dá)基因進(jìn)行GO 功能富集分析，共注釋到431 個(gè)GO terms，其中生物過程（biological process）分類富集到256 個(gè)GO terms，細(xì)胞組成分類富集到75 個(gè)GO terms，分子功能分類富集到105 個(gè)GO terms。并挑選了富集程度最靠前的31 個(gè)GO terms（圖1），其中，顯著富集的GO terms 有10 個(gè)（Q-value＜0.05），生物過程類別中顯著富集的分別是：細(xì)胞過程（cellular process，GO：0009987）、代謝 過 程（metabolic process，GO：0008152）、生 物 調(diào) 控（biological regulation，GO：0065007）、生 物 過 程 調(diào) 控（regulation of biological process，GO：0050789）、發(fā)育過程（developmental process，GO：0032502）；細(xì)胞組成類別中顯著富集的分別是：細(xì)胞（cell，GO：0005623）、細(xì)胞部分（cell part，GO：0044464）、細(xì)胞器（organelle，GO：0043226）；分子功能類別中顯著富集的分別是：黏合物（binding，GO：0005488）、核酸結(jié)合轉(zhuǎn)錄因子活性（nucleic acid binding transcription factor activity，GO：0001071）。

圖1 “滇北”淹水0 h vs 8 h 差異miRNA 靶基因GO 分類Fig.1 Histogram of GO classification of differential miRNA target genes for 0 h vs 8 h flooding in“Dianbei”1：細(xì)胞過程Cellular process；2：代謝過程Metabolic process；3：生物調(diào)控Metabolic process；4：生物過程調(diào)控Regulation of biological process；5：?jiǎn)我挥袡C(jī)體過程Single-organism process；6：發(fā)育過程Developmental process；7：對(duì)刺激的反應(yīng)Response to stimulus；8：定位Localization；9：生物過程的正向調(diào)節(jié)Positive regulation of biological process；10：信號(hào)Signaling；11：多細(xì)胞生物的過程Multicellular organismal process；12：細(xì)胞組成組織或生物發(fā)生Cellular component organization or biogenesis；13：生物過程的負(fù)調(diào)控Negative regulation of biological process；14：繁殖Reproduction；15：生殖過程Reproductive process；16：細(xì)胞Cell；17：細(xì)胞部分Cell part；18：細(xì)胞器Organelle；19：膜Membrane；20：膜部分Membrane part；21：細(xì)胞器部分Organelle part；22：大分子復(fù)合體Macromolecular complex；23：細(xì)胞外區(qū)域Extracellular region；24：超分子絡(luò)合物Supramolecular complex；25：細(xì)胞連接Cell junction；26：合胞體Symplast；27：黏合物Binding；28：催化活性Catalytic activity；29：核酸結(jié)合轉(zhuǎn)錄因子活性Nucleic acid binding transcription factor activity；30：結(jié)構(gòu)分子活性Structural molecule activity；31：轉(zhuǎn)運(yùn)活力Transporter activity.

對(duì)水淹脅迫下“滇北”鴨茅0 h vs 24 h 的66 個(gè)差異表達(dá)基因進(jìn)行GO 功能富集分析，共注釋到404 個(gè)GO terms 上，其中生物過程分類富集到242 個(gè)GO terms，細(xì)胞組成分類富集到73 個(gè)GO terms，分子功能分類富集到94 個(gè)GO terms。并挑選了富集程度最靠前的29 個(gè)GO terms（圖2），其中，顯著富集的GO terms 有11 個(gè)（Qvalue ＜0.05），其中生物過程類別中最顯著富集的分別是：細(xì)胞過程（cellular process，GO：0009987）、代謝過程（metabolic process，GO：0008152）、生物調(diào)節(jié)（biological regulation，GO：0065007）、生物過程調(diào)控（regulation of biological process，GO：0050789）、發(fā)育過程（developmental process，GO：0032502）；細(xì)胞組成類別中最顯著富集的分別是：細(xì)胞（cell，GO：0005623）、細(xì)胞部分（cell part，GO：0044464）、細(xì)胞器（organelle，GO：0043226）、超分子絡(luò)合物（supramolecular complex，GO：0099080）；分子功能類別中最顯著富集的分別是：黏合物（binding，GO：0005488）、核酸結(jié)合轉(zhuǎn)錄因子活性（nucleic acid binding transcription factor activity，GO：0001071）。

圖2 “滇北”淹水0 h vs 24 h 差異miRNA 靶基因GO 分類Fig.2 Histogram of GO classification of differential miRNA target genes for 0 h vs 24 h flooding in“Dianbei”1：細(xì)胞過程Cellular process；2：代謝過程Metabolic process；3：生物調(diào)節(jié)Biological regulation；4：生物過程調(diào)控Regulation of biological process；5：?jiǎn)我挥袡C(jī)體過程Single-organism process；6：發(fā)育過程Developmental process；7：對(duì)刺激的反應(yīng)Response to stimulus；8：定位Localization；9：信號(hào)Signaling；10：生物過程的正向調(diào)節(jié)Positive regulation of biological process；11：多細(xì)胞生物的過程Multicellular organismal process；12：細(xì)胞組成組織或生物發(fā)生Cellular component organization or biogenesis；13：生物過程的負(fù)調(diào)控Negative regulation of biological process；14：繁殖Reproduction；15：生殖過程Reproductive process；16：細(xì)胞Cell；17：細(xì)胞部分Cell part；18：細(xì)胞器Organelle；19：膜Membrane；20：膜部分Membrane part；21：細(xì)胞器部分Organelle part；22：高分子復(fù)合物Macromolecular complex；23：超分子絡(luò)合物Supramolecular complex；24：細(xì)胞外區(qū)域Extracellular region；25：黏合物Binding；26：催化活性Catalytic activity；27：核酸結(jié)合轉(zhuǎn)錄因子活性Nucleic acid binding transcription factor activity；28：結(jié)構(gòu)分子活性Structural molecule activity；29：轉(zhuǎn)運(yùn)活力Transporter activity.

綜上可知，“滇北”鴨茅水淹0 h vs 8 h 和0 h vs 24 h 中，GO 3 個(gè)組分中富集到較多靶基因的條目基本相同。但注釋到的GO term 數(shù)目以及富集在各條目上的靶基因，8 h 多于24 h，這表明“滇北”在水淹脅迫8 h 響應(yīng)比24 h更復(fù)雜；推測(cè)這些鴨茅差異miRNA 可能通過調(diào)控參與以上過程的靶基因的表達(dá)來影響鴨茅葉片發(fā)育及響應(yīng)環(huán)境。

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

在生物體內(nèi)，不同基因相互協(xié)調(diào)行使其生物學(xué)功能，通過Pathway 顯著性富集能確定差異表達(dá)基因參與的最主要生化代謝途徑和信號(hào)轉(zhuǎn)導(dǎo)途徑。KEGG（kyoto encyclopedia of genes and genomes）是有關(guān)Pathway 的主要公共數(shù)據(jù)庫(kù)［44］，該數(shù)據(jù)庫(kù)提供的整合代謝途徑（pathway）有助于進(jìn)一步解讀基因的功能，以Q-value＜0.05 的pathway 為差異表達(dá)顯著富集的pathway。在富集散點(diǎn)圖中Rich factor 指該pathway 中富集到的差異基因個(gè)數(shù)與注釋基因個(gè)數(shù)的比值；Rich factor 越大，表示富集的程度越大；Q-value 是進(jìn)行多重假設(shè)檢驗(yàn)校正之后的P-value。

“滇北”在水淹脅迫0 h vs 8 h 有73 個(gè)差異表達(dá)基因注釋到19 條pathway（圖3），發(fā)現(xiàn)差異miRNA 的靶基因主要參與IL-17 信號(hào)通路（IL-17 signaling pathway）、Th17 細(xì)胞分化（Th17 cell differentiation）、雌激素信號(hào)通路（estrogen signaling pathway）、NOD-樣 受 體 信 號(hào) 通 路（NOD-like receptor signaling pathway）、壞 死 性 凋 亡（necroptosis）、PI3K-Akt 信號(hào) 通 路（PI3K-Akt signaling pathway）、植物?病原體相 互作用（plant?pathogen interaction）、抗原處理和呈遞（antigen processing and presentation）、內(nèi)質(zhì)網(wǎng)中的蛋白質(zhì)加工（protein processing in endoplasmic reticulum）。

圖3 “滇北”水淹0 h vs 8 h 差異miRNA 靶基因KEGG 富集Fig.3 KEGG enrichment scatter plot of differential miRNA gene KEGG after 0 h vs 8 h flooding in“Dianbei”

“滇北”在水淹脅迫0 h vs 24 h 有69 個(gè)差異表達(dá)基因注釋到48 條pathway，前20 條如圖4 所示，發(fā)現(xiàn)差異miRNA 的靶基因主要參與Th17 細(xì)胞分化（Th17 cell differentiation）、IL-17 信號(hào)通路（IL-17 signaling pathway）、雌激素信號(hào)通路（estrogen signaling pathway）、NOD-樣受體信號(hào)通路（NOD-like receptor signaling pathway）、壞死性凋亡（necroptosis）、PI3K-Akt 信號(hào)通路（PI3K-Akt signaling pathway）、植物?病原體相互作用（plant?pathogen interaction）、抗原處理和呈遞（antigen processing and presentation）、內(nèi)質(zhì)網(wǎng)中的蛋白質(zhì)加工（protein processing in endoplasmic reticulum）。

圖4 “滇北”水淹0 h vs 24 h 差異miRNA 靶基因KEGG 富集Fig.4 KEGG enrichment scatter plot of differential miRNA gene KEGG after 0 h vs 24 h flooding in“Dianbei”

3 討論

水淹脅迫會(huì)導(dǎo)致植物根系缺氧，并引起代謝從有氧呼吸向厭氧發(fā)酵轉(zhuǎn)變，從而導(dǎo)致基因轉(zhuǎn)錄和蛋白質(zhì)合成的快速變化［39］。前期研究結(jié)果表明，在鴨茅對(duì)淹水脅迫的響應(yīng)過程中，多個(gè)差異基因和生物學(xué)過程參與了水淹的響應(yīng)，特別是活性氧（reactive oxygen species，ROS）的清除過程［33］，在本研究中，鴨茅水淹處理0、8 和24 h 共鑒定到208 個(gè) 差 異 表 達(dá)的miRNA，差 異 表 達(dá)的miRNA 主 要屬 于miR166、miR159、miR167、miR396 和miR156 這5 個(gè)miRNA 家族，這些結(jié)果表明“滇北”鴨茅中廣泛存在各類保守miRNA，可能對(duì)“滇北”適應(yīng)水淹脅迫后的缺氧環(huán)境起著重要的調(diào)節(jié)作用。在植物對(duì)外界環(huán)境響應(yīng)過程中，miRNA 通過對(duì)靶基因的調(diào)控來影響植物的生理功能［45］。基于對(duì)差異miRNA 進(jìn)行靶基因預(yù)測(cè)及靶基因的GO 和KEGG 功能分析，發(fā)現(xiàn)這些靶基因主要參與細(xì)胞生理過程（cellular process）、代謝過程（metabolic process）、IL-17 信號(hào)通路（IL-17 signaling pathway）、Th17 細(xì)胞分化（Th17 cell differentiation）等植物逆境響應(yīng)過程，表明“滇北”鴨茅差異表達(dá)miRNA 通過調(diào)節(jié)對(duì)環(huán)境適應(yīng)相關(guān)靶基因的表達(dá)，影響不同時(shí)間點(diǎn)葉片的抗逆性。相關(guān)研究發(fā)現(xiàn)水淹脅迫能夠影響植物miRNA 的形成，部分響應(yīng)水淹脅迫的miRNA 在植物中被鑒定出來，這些miRNA 在水淹脅迫中的調(diào)控功能也得到證實(shí)［40］。

miR166 通過調(diào)控靶基因HD-ZipⅢ參與調(diào)控植物的生長(zhǎng)發(fā)育過程，也可以作用于HD-ZipⅢmRNA 參與調(diào)節(jié)擬南芥逆境脅迫響應(yīng)［46?47］。在玉米中，水淹脅迫下zam-miR166 的積累，影響根分生組織細(xì)胞的分化從而導(dǎo)致維管系統(tǒng)的形成，木質(zhì)部組織擴(kuò)張［48］。有研究表明，苜蓿（Medicago truncatula）中過表達(dá)的miR166 會(huì)導(dǎo)致側(cè)根數(shù)量形成減少和根維管束的變化，并且miR166 是少數(shù)在低氧脅迫下降低的家族之一［49］。擬南芥在遭遇澇漬脅迫后，miR166 依賴于線粒體功能，在基因表達(dá)調(diào)控和低氧環(huán)境反應(yīng)中發(fā)揮作用［50］。也有關(guān)于miR166 在其他逆境脅迫中的報(bào)道，如經(jīng)高溫處理后的不同品種蕹菜（Ipomoea aquatica），表現(xiàn)出截然相反的表達(dá)模式，是由于其在長(zhǎng)時(shí)間高溫后莖頂端分生組織受損造成［51］。水淹脅迫后“滇北”鴨茅差異表達(dá)的miRNA 中，miR166 數(shù)目最多，推測(cè)miR166 可能在植物遭受水淹脅迫中整合關(guān)鍵的缺氧信號(hào)途徑，以此抵御該非生物脅迫。

miR167 可通過調(diào)控不同的靶基因參與到植物發(fā)育和逆境脅迫適應(yīng)的不同階段中。miR167 的靶基因是生長(zhǎng)激素反應(yīng)因子ARF基因家族中的ARF6和ARF8，生物學(xué)功能是增強(qiáng)生長(zhǎng)素反應(yīng)活性以及側(cè)根生長(zhǎng)，miR167 與其靶基因ARF6和ARF8之間存在反饋調(diào)控，ARF6激活miR167 的表達(dá)，ARF8抑制miR167 的表達(dá)，它們之間的相互激活或抑制調(diào)節(jié)著不定根的形成［52］。白羽扇豆（Lupinus albus）在低磷脅迫條件下，miR167 靶向調(diào)節(jié)轉(zhuǎn)錄因子增強(qiáng)生長(zhǎng)素反應(yīng)活性，促使側(cè)根生長(zhǎng)上調(diào)表達(dá)［53］。張丹鳳［54］對(duì)水淹脅迫下玉米幼苗的適應(yīng)機(jī)制進(jìn)行研究發(fā)現(xiàn)，水淹脅迫影響玉米幼苗根系細(xì)胞中生長(zhǎng)素（IAA）的動(dòng)態(tài)平衡。IAA 含量的變化會(huì)誘導(dǎo)zma-miR167 的表達(dá)，zmamiR167 通過與AGO1相互作用，負(fù)調(diào)控靶基因ARF17的表達(dá)。ARF17抑制基因GH3的表達(dá)，進(jìn)一步誘導(dǎo)下游基因ra2的表達(dá)，最終導(dǎo)致不定根的產(chǎn)生?！暗岜薄兵喢┰诓煌奶幚頃r(shí)間段下，miR167 差異表達(dá)基因的數(shù)量為8 h少于24 h，這表明miR167 受水淹脅迫后，不定根的形成與鴨茅水淹脅迫密切相關(guān)。

miR156，miR159 和miR396 這3 個(gè)miRNA 是植物發(fā)育階段過渡的關(guān)鍵調(diào)控因子［55］，在非生物脅迫中也有重要作用。miR156 是在植物中發(fā)現(xiàn)的第一個(gè)miRNA，其調(diào)節(jié)發(fā)育階段轉(zhuǎn)換的功能最初在擬南芥中得到證實(shí)，它通過靶向表達(dá)SPL基因來實(shí)現(xiàn)［56］。SPL是一種植物特異性轉(zhuǎn)錄因子，具有保守的SBP 結(jié)構(gòu)域，可識(shí)別并特異性結(jié)合下游靶基因的啟動(dòng)子區(qū)域，從而調(diào)節(jié)植物的生長(zhǎng)發(fā)育［57］。在擬南芥中，miR156b 的過表達(dá)會(huì)產(chǎn)生更多的葉片并且迫使開花延遲，同時(shí)伴隨著多個(gè)SPL基因（SPL2/3/4/5/6/9/10/11/13/15）表達(dá)的下調(diào)［56］。miR156 在植物暴露于各種環(huán)境脅迫時(shí)異常表達(dá)，研究表明miR156 是植物熱應(yīng)激記憶所必需的，并且植物表現(xiàn)出對(duì)熱應(yīng)激的耐受性增強(qiáng)［58］。在水淹脅迫方面，CSA-miR156G-3p 在漬水的黃瓜（Cucumis sativus）下胚軸中下調(diào)，參與植物的非生物脅迫響應(yīng)［59］。本研究中，8 h 差異表達(dá)的miR156 多于24 h，說明8 h 時(shí)鴨茅為降低水淹脅迫對(duì)植物生長(zhǎng)發(fā)育的影響，通過表達(dá)更多相關(guān)基因來適應(yīng)該非生物脅迫。

miR159 的靶基因編碼MYB，Ma 等［55］研究發(fā)現(xiàn)miR159 的上調(diào)可以抑制MYB29 的表達(dá)，降低脂肪族脂肪硫代葡萄糖苷（giucosinolates，GSLS）水平，從而調(diào)節(jié)莖向根莖轉(zhuǎn)化過程中的生長(zhǎng)，以響應(yīng)地上莖到根莖的轉(zhuǎn)變。在玉米的根中，miR159 被水澇誘導(dǎo)上調(diào)，并負(fù)責(zé)沉默編碼GAMYBs、MYB33 和MYB101 同系物的兩個(gè)mRNA［22］。有研究表明，miR159 對(duì)各種環(huán)境脅迫有反應(yīng)，與野生型相比，過量表達(dá)miR159 的轉(zhuǎn)基因水稻（Oryza sativa）對(duì)熱脅迫更敏感，這表明miR159 的下調(diào)可能有助于增強(qiáng)水稻對(duì)熱脅迫耐受性［60］。Zhang 等［48］發(fā)現(xiàn)Zam-miR159 在淹水初期被抑制，24 h 后被誘導(dǎo)。miRNA 的預(yù)測(cè)涉及碳水化合物和能量代謝，包括淀粉合成酶、轉(zhuǎn)化酶、蘋果酸酶和ATP 酶，表明淹水響應(yīng)miRNA 在轉(zhuǎn)錄后水平上參與調(diào)節(jié)玉米根系的代謝、生理和形態(tài)適應(yīng)。水淹脅迫處理后miR159 明顯變化且顯著上調(diào)，說明其與植物耐澇性的發(fā)揮密切相關(guān)。

miR396 通過調(diào)控其靶基因生長(zhǎng)調(diào)節(jié)因子（growth regulating factor，GRF）在植物不同發(fā)育階段發(fā)揮多種作用。Zhang 等［61］通過對(duì)發(fā)育中的玉米籽粒進(jìn)行RNA 測(cè)序，分析miR396 和GRF基因在玉米有效灌漿過程中的動(dòng)態(tài)變化，發(fā)現(xiàn)在該過程中隨著玉米灌漿進(jìn)程的推進(jìn)，miR396 高表達(dá)并逐漸下降，與其靶基因GRF呈顯著負(fù)相關(guān)，且大多持續(xù)增加。miR396 是一種已知的應(yīng)激反應(yīng)miRNA，在鹽和堿脅迫處理期間，結(jié)構(gòu)性過表達(dá)miR396 的轉(zhuǎn)基因擬南芥和水稻植株顯著阻礙了根的生長(zhǎng)，并影響植株的生長(zhǎng)和發(fā)育［62］。在本研究中，miR156、miR159 和miR396 在水淹脅迫處理后高表達(dá)，這可能是由于水淹脅迫限制了鴨茅正常生長(zhǎng)發(fā)育的關(guān)鍵調(diào)控因子，意味著“滇北”鴨茅需要更多更強(qiáng)的抗性基因來維持植物的存活。

植物在遭遇水淹脅迫時(shí)會(huì)導(dǎo)致基因的表達(dá)變化，可通過生物信息學(xué)分析對(duì)DEGs 進(jìn)行功能預(yù)測(cè)，探究植物耐澇機(jī)理和分子調(diào)控網(wǎng)絡(luò)。Xu 等［59］通過研究黃瓜不定根的形成對(duì)內(nèi)澇脅迫的影響，GO 富集分析表明，在細(xì)胞氧化還原穩(wěn)態(tài)、細(xì)胞骨架、光合作用和細(xì)胞生長(zhǎng)等方面功能的靶基因被過度表達(dá)，說明細(xì)胞相關(guān)生物過程在響應(yīng)水淹脅迫中發(fā)揮重要作用。林延慧等［63］通過轉(zhuǎn)錄組和蛋白組的關(guān)聯(lián)分析，在大豆（Glycine max）中篩選出3 個(gè)與水淹脅迫絲裂原活化蛋白激酶（mitogen-activated protein kinase，MAPK）信號(hào)通路相關(guān)的差異表達(dá)基因，分別是Glyma. 07G127100、Glyma. 05G124000 和Glyma. 05G123700，發(fā)現(xiàn)水淹脅迫引發(fā)了植物類似于病原菌侵染的基本防御反應(yīng)，傳遞信號(hào)通過FLS2 蛋白進(jìn)一步激活了MAPK 信號(hào)通路來抵御澇害脅迫。

4 結(jié)論

本研究利用小RNA 測(cè)序方法，對(duì)水淹脅迫處理后“滇北”鴨茅miRNA 在不同時(shí)間點(diǎn)的表達(dá)模式及其靶基因進(jìn)行了分析。初步發(fā)現(xiàn)，miRNA 在水淹脅迫不同時(shí)間點(diǎn)表現(xiàn)出一定的動(dòng)態(tài)變化。同時(shí)，根據(jù)miRNA 序列及靶標(biāo)功能預(yù)測(cè)，推測(cè)鴨茅miR166、miR159、miR167、miR396 和miR156 很可能參與了鴨茅的水淹脅迫應(yīng)答過程。這些結(jié)果為鴨茅中與水淹脅迫相關(guān)miRNA 的進(jìn)一步研究以及探究miRNA 在鴨茅抗?jié)尺^程中的潛在功能和作用機(jī)理提供了參考。但有關(guān)響應(yīng)水淹脅迫的miRNA 如何調(diào)控其靶基因，預(yù)測(cè)的耐澇相關(guān)靶基因又如何發(fā)揮其作用等問題還需要進(jìn)一步探究。