鴨茅響應水淹脅迫的miRNA 差異表達分析

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

（1. 西南大學動物科學技術學院，重慶402460；2. 草食動物科學重慶市重點實驗室，重慶 400715；3. 福建傲農生物科技集團股份有限公司，福建 廈門 361111；4. 四川農業大學草業科技學院，四川 成都 611130）

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

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

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

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

鑒于此，本研究以耐澇性較強的“滇北”鴨茅為試驗材料，利用高通量測序與生物信息學技術，篩選與鴨茅水淹抗性相關的miRNA，分析它們在水淹脅迫中的差異表達，旨在為今后研究鴨茅抗澇作用機制及抗澇遺傳改良工作提供參考，同時為耐澇草本植物新品種的選育提供借鑒。

1 材料與方法

1.1 試驗材料與處理

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

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

1.2 小RNA 文庫的構建與高通量測序

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

1.3 測序原始數據預處理

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

1.4 響應水淹脅迫miRNA 的鑒定

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

miRNA 前體的標志性發夾結構能夠用來預測新的miRNA。基本原理是通過截取一定長度sRNA 比對上的參考序列，通過探尋其二級結構及Dicer 酶切位點信息、能量等特征進行分析，預測樣品中新miRNA，并進行各樣本中匹配上的sRNA 的序列、長度、出現的次數等信息，以及不同長度miRNA 的首位點堿基分布和所有miRNA的各位點堿基分布情況的統計。對于各個數據庫中均未注釋但可比對到基因組序列的sRNA，將作為候選的新miRNA。采用mirdeep2［35］軟件進行新miRNA 的分析，綜合考慮到sRNA 測序序列表達豐度，前體miRNA 最低自由能和復雜的二級結構，預測得到的候選新miRNA 前體的莖環結構圖，完成新miRNA 的預測。

1.5 miRNA 差異表達分析

基因差異表達分析的輸入數據為miRNA 表達水平分析中得到的miRNA 對應的read counts 數據。用R 語言包DESeq2 進行差異分析，篩選閾值為|log2 fold change|＞1&P＜0.05。對于之前得到的已知miRNA 和新miRNA，做miRNA 靶基因預測分析。

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

2 結果與分析

2.1 鴨茅小RNA 測序數據及質量控制

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

表1 不同鴨茅樣品測序數據統計Table 1 Statistical results of sequencing data of different D.glomerata samples

2.2 鴨茅小RNA 分類注釋

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

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

表2 小RNA 分類統計Table 2 Small RNA classification statistics（%）

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

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

表3 差異基因統計Table 3 Statistical of differentially expressed genes（DEGs）

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

2.4 鴨茅miRNA 的靶基因預測

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

表5 不同的靶基因對應同一個miRNATable 5 Different target genes correspond to the same miRNA

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

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

本研究對水淹脅迫下“滇北”鴨茅0 h vs 8 h 的70 個差異表達基因進行GO 功能富集分析，共注釋到431 個GO terms，其中生物過程（biological process）分類富集到256 個GO terms，細胞組成分類富集到75 個GO terms，分子功能分類富集到105 個GO terms。并挑選了富集程度最靠前的31 個GO terms（圖1），其中，顯著富集的GO terms 有10 個（Q-value＜0.05），生物過程類別中顯著富集的分別是：細胞過程（cellular process，GO：0009987）、代謝 過 程（metabolic process，GO：0008152）、生 物 調 控（biological regulation，GO：0065007）、生 物 過 程 調 控（regulation of biological process，GO：0050789）、發育過程（developmental process，GO：0032502）；細胞組成類別中顯著富集的分別是：細胞（cell，GO：0005623）、細胞部分（cell part，GO：0044464）、細胞器（organelle，GO：0043226）；分子功能類別中顯著富集的分別是：黏合物（binding，GO：0005488）、核酸結合轉錄因子活性（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：細胞過程Cellular process；2：代謝過程Metabolic process；3：生物調控Metabolic process；4：生物過程調控Regulation of biological process；5：單一有機體過程Single-organism process；6：發育過程Developmental process；7：對刺激的反應Response to stimulus；8：定位Localization；9：生物過程的正向調節Positive regulation of biological process；10：信號Signaling；11：多細胞生物的過程Multicellular organismal process；12：細胞組成組織或生物發生Cellular component organization or biogenesis；13：生物過程的負調控Negative regulation of biological process；14：繁殖Reproduction；15：生殖過程Reproductive process；16：細胞Cell；17：細胞部分Cell part；18：細胞器Organelle；19：膜Membrane；20：膜部分Membrane part；21：細胞器部分Organelle part；22：大分子復合體Macromolecular complex；23：細胞外區域Extracellular region；24：超分子絡合物Supramolecular complex；25：細胞連接Cell junction；26：合胞體Symplast；27：黏合物Binding；28：催化活性Catalytic activity；29：核酸結合轉錄因子活性Nucleic acid binding transcription factor activity；30：結構分子活性Structural molecule activity；31：轉運活力Transporter activity.

對水淹脅迫下“滇北”鴨茅0 h vs 24 h 的66 個差異表達基因進行GO 功能富集分析，共注釋到404 個GO terms 上，其中生物過程分類富集到242 個GO terms，細胞組成分類富集到73 個GO terms，分子功能分類富集到94 個GO terms。并挑選了富集程度最靠前的29 個GO terms（圖2），其中，顯著富集的GO terms 有11 個（Qvalue ＜0.05），其中生物過程類別中最顯著富集的分別是：細胞過程（cellular process，GO：0009987）、代謝過程（metabolic process，GO：0008152）、生物調節（biological regulation，GO：0065007）、生物過程調控（regulation of biological process，GO：0050789）、發育過程（developmental process，GO：0032502）；細胞組成類別中最顯著富集的分別是：細胞（cell，GO：0005623）、細胞部分（cell part，GO：0044464）、細胞器（organelle，GO：0043226）、超分子絡合物（supramolecular complex，GO：0099080）；分子功能類別中最顯著富集的分別是：黏合物（binding，GO：0005488）、核酸結合轉錄因子活性（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：細胞過程Cellular process；2：代謝過程Metabolic process；3：生物調節Biological regulation；4：生物過程調控Regulation of biological process；5：單一有機體過程Single-organism process；6：發育過程Developmental process；7：對刺激的反應Response to stimulus；8：定位Localization；9：信號Signaling；10：生物過程的正向調節Positive regulation of biological process；11：多細胞生物的過程Multicellular organismal process；12：細胞組成組織或生物發生Cellular component organization or biogenesis；13：生物過程的負調控Negative regulation of biological process；14：繁殖Reproduction；15：生殖過程Reproductive process；16：細胞Cell；17：細胞部分Cell part；18：細胞器Organelle；19：膜Membrane；20：膜部分Membrane part；21：細胞器部分Organelle part；22：高分子復合物Macromolecular complex；23：超分子絡合物Supramolecular complex；24：細胞外區域Extracellular region；25：黏合物Binding；26：催化活性Catalytic activity；27：核酸結合轉錄因子活性Nucleic acid binding transcription factor activity；28：結構分子活性Structural molecule activity；29：轉運活力Transporter activity.

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

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

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

“滇北”在水淹脅迫0 h vs 8 h 有73 個差異表達基因注釋到19 條pathway（圖3），發現差異miRNA 的靶基因主要參與IL-17 信號通路（IL-17 signaling pathway）、Th17 細胞分化（Th17 cell differentiation）、雌激素信號通路（estrogen signaling pathway）、NOD-樣 受 體 信 號 通 路（NOD-like receptor signaling pathway）、壞 死 性 凋 亡（necroptosis）、PI3K-Akt 信號 通 路（PI3K-Akt signaling pathway）、植物?病原體相 互作用（plant?pathogen interaction）、抗原處理和呈遞（antigen processing and presentation）、內質網中的蛋白質加工（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 個差異表達基因注釋到48 條pathway，前20 條如圖4 所示，發現差異miRNA 的靶基因主要參與Th17 細胞分化（Th17 cell differentiation）、IL-17 信號通路（IL-17 signaling pathway）、雌激素信號通路（estrogen signaling pathway）、NOD-樣受體信號通路（NOD-like receptor signaling pathway）、壞死性凋亡（necroptosis）、PI3K-Akt 信號通路（PI3K-Akt signaling pathway）、植物?病原體相互作用（plant?pathogen interaction）、抗原處理和呈遞（antigen processing and presentation）、內質網中的蛋白質加工（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 討論

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

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

miR167 可通過調控不同的靶基因參與到植物發育和逆境脅迫適應的不同階段中。miR167 的靶基因是生長激素反應因子ARF基因家族中的ARF6和ARF8，生物學功能是增強生長素反應活性以及側根生長，miR167 與其靶基因ARF6和ARF8之間存在反饋調控，ARF6激活miR167 的表達，ARF8抑制miR167 的表達，它們之間的相互激活或抑制調節著不定根的形成［52］。白羽扇豆（Lupinus albus）在低磷脅迫條件下，miR167 靶向調節轉錄因子增強生長素反應活性，促使側根生長上調表達［53］。張丹鳳［54］對水淹脅迫下玉米幼苗的適應機制進行研究發現，水淹脅迫影響玉米幼苗根系細胞中生長素（IAA）的動態平衡。IAA 含量的變化會誘導zma-miR167 的表達，zmamiR167 通過與AGO1相互作用，負調控靶基因ARF17的表達。ARF17抑制基因GH3的表達，進一步誘導下游基因ra2的表達，最終導致不定根的產生。“滇北”鴨茅在不同的處理時間段下，miR167 差異表達基因的數量為8 h少于24 h，這表明miR167 受水淹脅迫后，不定根的形成與鴨茅水淹脅迫密切相關。

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

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

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

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

4 結論

本研究利用小RNA 測序方法，對水淹脅迫處理后“滇北”鴨茅miRNA 在不同時間點的表達模式及其靶基因進行了分析。初步發現，miRNA 在水淹脅迫不同時間點表現出一定的動態變化。同時，根據miRNA 序列及靶標功能預測，推測鴨茅miR166、miR159、miR167、miR396 和miR156 很可能參與了鴨茅的水淹脅迫應答過程。這些結果為鴨茅中與水淹脅迫相關miRNA 的進一步研究以及探究miRNA 在鴨茅抗澇過程中的潛在功能和作用機理提供了參考。但有關響應水淹脅迫的miRNA 如何調控其靶基因，預測的耐澇相關靶基因又如何發揮其作用等問題還需要進一步探究。