甘蔗耐寒相關miRNA的生物信息學分析及其靶基因預測

朱鵬錦， 宋奇琦， 譚秦亮， 程 琴， 李佳慧， 龐新華， 周全光， 呂 平， 歐克緯， 盧業飛， 農澤梅， 唐桓偉， 龍盛風

(廣西壯族自治區亞熱帶作物研究所， 南寧 530001 )

miRNA是一類含有20～24個堿基的小分子內源性非編碼RNA，具有高度的保守性、時序性和組織特異性(Ambros, 2004)。miRNA主要在轉錄后水平上通過介導靶基因mRNA的切割或抑制翻譯來調節基因的表達，在生物體代謝過程中起到多種調控作用，如參與調控植物器官的形態建成(Sunkar et al., 2012)、生長發育(Thiebaut, et al.， 2012)、激素分泌、信號轉導以及對外界環境脅迫(Xiong & Zhu, 2003)的響應等過程。在擬南芥和水稻響應低溫的研究中發現，miR-167、miR-169、miR-319和miR-171等miRNA家族在低溫響應具有重要的生物學作用(Wang et al., 2010)，其中miR-169的靶基因為低溫誘導的重要基因CBF (Sunkar et al., 2007)；還在水稻響應低溫脅迫的研究中證明，miRNA-319和miRNA-171的靶基因屬于MYB類的轉錄因子，兩者的表達量互為消長關系，進而說明miRNA在水稻耐低溫途徑中所起的調控作用(Lü et al., 2010)。因此，準確高效分離和鑒定miRNA及其靶基因并分析它們的功能，更精確地掌握miRNA在植物抗逆脅迫過程中的調控機制，是目前植物miRNA研究領域的重要內容。

甘蔗是熱帶、亞熱帶地區重要的糖料作物，低溫是限制其擴大種植區域和實現高產穩產的重要因素之一(吳棉國等, 2010)。出現的偶發性大范圍冰凍雨雪或嚴重的霜凍天氣會導致我國甘蔗主產區遭受巨大的經濟損失(何燕等, 2009; 匡昭敏等, 2009)；廣西地理環境特殊，冬季持續低溫以及春季“倒春寒”形成的陰雨霜凍頻發導致大面積甘蔗受到冷害，具體表現在甘蔗葉片枯萎、莖桿壞死以致蔗糖含量、甘蔗產量下降(李楊瑞等, 2011)。古麗等(2011)采集甘蔗主產市縣40多年的氣象數據，結合甘蔗種植面積、甘蔗產量和蔗糖產量等指標進行分析發現，低溫霜凍災害是影響甘蔗種植、甘蔗生產和蔗糖產量的主要環境因子。因此，了解甘蔗的耐低溫調控機制是培育耐寒性強或適合我國熱帶北緣氣候的甘蔗品種的前提。

本研究在觀察不同甘蔗品種的田間農藝性狀和模擬低溫脅迫的生理生化研究基礎上，篩選出抗寒能力較強的甘蔗品種作為材料，利用高通量測序技術及生物信息學方法，獲得與甘蔗響應低溫脅迫相關的miRNA，分析miRNA的差異表達情況，明確miRNA與靶基因的作用關系，并對預測所得的靶基因進行基因本體(gene ontology, GO)分析，挖掘低溫脅迫應答基因，為選育耐寒性強的優良甘蔗新品種提供理論依據和技術支撐。

1 材料與方法

1.1 材料

以廣西壯族自治區農業科學院甘蔗研究所選育的‘桂糖28號’(GT28)、廣西蔗區主栽品種‘新臺糖22號’(ROC22)和廣西壯族自治區亞熱帶作物研究所新選育的‘桂熱2號’(GR2)為材料。

1.2 方法

1.2.1 實驗設計與管理 選擇無病蟲害、蔗莖大小均勻的種莖切成單芽段，用清水沖洗種莖、干布擦洗干凈；用50%多菌靈可濕性粉劑1 000倍液浸12 h消毒；用蒸餾水浸洗種莖1 min后，再用蒸餾水浸洗擰干的棉布將種莖包好并做好標記；用橡皮筋扎好放進溫度為25 ℃的恒溫箱催芽。當種莖萌芽并長出幼根時，將其移植到裝有營養土的塑料盆并做好標記。每盆1段種莖，盆高17.5 cm，盆寬16 cm。育苗期間，每株施用完全營養液2次，每次10 mL。當幼苗兩葉一心時，選擇長勢、大小均一的幼苗進行低溫處理。低溫脅迫處理溫度為4 ℃，光照強度為5 000 lx，處理24 h；對照(CK)溫度為28 ℃，光照強度為5 000 lx。

1.2.2 RNA 提取及 Illumina 測序 采集低溫處理和對照樣本的葉片，及時用TRNzol Reagent試劑盒(天根生化科技有限公司, 北京)提取并檢測RNA。將低溫脅迫處理組與對照樣本組提取獲得的高質量RNA送往華大基因科技有限公司(深圳)進行文庫構建，每處理3個重復。采用Illumina HiSeqTM 2000平臺對質量合格的RNA文庫進行測序(Ali et al., 2008)。通過去接頭、去低質量、去污染等過程完成數據處理獲得干凈序列后，使用Bowtie 2軟件按照MiRbase > pirnabank > snoRNA(human/plant) > Rfam > other sRNA的優先級順序將小RNA(sRNA)遍歷注釋獲得未注釋RNA片段。

1.2.3已知miRNA的鑒定及新miRNA預測 用miRDeep 2軟件對所得的未注釋sRNA序列與參考基因組 [割手密(Saccharumspontaneum)基因組]進行序列比對分析 (Maurits et al., 2015)，鑒定出已知的甘蔗耐寒相關miRNA；將過濾獲得的未比對序列比對參考序列，通過堿基數目延伸、miRNA結構預測方法獲得新的miRNA。

1.2.4 差異表達miRNA的鑒定 對耐寒型甘蔗樣品的已知miRNA的讀數進行分析，判斷低溫脅迫前后不同耐寒型甘蔗樣品中miRNA的差異表達，以|log 2(FoldChange)|≥1，P<0.05為篩選標準，獲得低溫脅迫前后差異表達 miRNA。

1.2.5 miRNA靶基因預測 用psRNATarget(Wu et al., 2012)、TargetFinder(Fahlgren & Carrington, 2010)和Tapirhybrid(Peer，2010)軟件進行差異表達miRNA的靶基因預測，并在GO數據庫對預測靶基因進行同源性搜索，確定其參與的信號傳導及生物代謝途徑。

1.2.6 實時熒光定量PCR驗證 選取14個差異表達miRNA及其靶基因，先用Premier 5.0設計其成熟miRNA特異正向引物、通用反向引物(表4)，再用LightCycler?480 Instrument II進行 RT-PCR擴增，設陰性對照(不添加cDNA模板)以監控可能的污染。以甘蔗GAPDH為內參基因, 每個樣品均設置3次重復, 采用2-ΔΔCt法計算基因相對表達。

2 結果與分析

2.1 高通量測序數據分析

選擇3個耐寒性不同甘蔗品種的幼苗進行低溫脅迫，每個品種的對照(CK)和低溫處理(T)3個重復，分別取葉片共18個樣本進行高通量測序，構建低溫脅迫前后sRNA文庫；對原始測序數據進行3′端去接頭、Trim低質量、片段大小選擇等質控處理，獲取高質量clean read序列(表1)。從表1可以看出，在GR2的CK和T葉片中分別挖掘出24 310 558、23 925 673條測序數據，最終處理分別得到21 343 194、21 576 594條，高質量clean read序列分別占總序列的87.83%、90.25%，且所構建的文庫與參考基因組進行比對，其CK和T分別有86.71%與84.32%的序列與參考基因組匹配；GT28的CK和T葉片中分別挖掘出23 915 244、23 528 103條測序數據，最終處理分別得到22 045 186、21 481 129條，高質量clean read序列分別占總序列的92.18%與91.28%，所構建的文庫與參考基因組進行比對，其CK和T分別有71.75%與84.10%的序列與參考基因組匹配；ROCC22的CK和T葉片中分別挖掘出27 069 867、23 631 208條測序數據，最終處理分別得到21 007 371、20 640 357條，高質量clean read序列占總序列的79.50%與87.35%，所構建的文庫與參考基因組進行比對，其CK和T分別有92.62%與85.67%的序列與參考基因組匹配(表1)。

表 1 樣品測序數據統計Table 1 Sequencing data statistics of samples

sRNA長度與不同功能有關，其中21～22 nt的主要與mRNA切割和轉錄后基因沉默相關，24 nt的主要與RNA導向的DNA甲基化和轉錄基因沉默相關。對total clean reads進行統計發現，sRNA主要集中在21～24 nt之間，但不同抗寒能力的材料之間有一定差異，且不同長度的測序頻率不同(圖1)。

圖 1 sRNA長度分布圖Fig. 1 Graphs sRNA length distribution

根據基因表達量進一步對每個樣品進行Pearson相關系數計算，并將這些系數以熱圖的形式反映出來(圖2)。從Pearson相關系數可以看出，同一品種相同處理的重復樣品間相關性較高(0.667～0.990)，表明所有樣品的表達量基本保持一致，說明其符合樣本重復性實驗標準，可以滿足后續的差異表達分析。

X、Y軸均代表每個樣品。顏色代表相關性系數，顏色越藍代表相關性越高，顏色越淺代表相關性越低。X and Y axes represent each sample. Color represents the correlation coefficient, the bluer the color, the higher the correlation, and the lighter the color, the lower the correlation.圖 2 樣品間相關性分析熱圖Fig. 2 Correlation heatmap of different samples

2.2 已知miRNA分析及新miRNA預測

用AASRA軟件將核苷酸序列比對到參考基因組及miRBase數據庫，鑒定得到分屬于84個已知miRNA家族的322個miRNA，其中對照組的297個miRNA分屬于69個家族，低溫處理組的305個miRNA分屬于74個家族，這些家族中擁有最多家族成員數量的為miR169(32個)，其次分別為miR166、miR171、miR167、miR156、miR396 (圖3)。

圖 3 miRNA家族數量分布前18位Fig. 3 Top 18 miRNA family members

根據成熟植物miRNA序列的高度保守型及miRNA前體擁有標志性發夾結構的特性，利用同源搜尋比對的方法在割手密miRBase數據庫中進行搜索，對具有同源性的miRNA序列進行提取，選取其中表達量最高的作為甘蔗響應低溫的保守miRNA的候選者，最終得到甘蔗響應低溫的110個新miRNA。

2.3 差異表達miRNA分析

通過比較對照組與低溫處理組miRNA的表達量變化情況，可以判斷低溫脅迫下抗寒能力不同的甘蔗品種miRNA的差異表達情況。在差異表達miRNA 檢測過程中，以|log 2(FC)|≥2，P<0.05作為篩選標準，低溫脅迫時差異表達miRNA 如表2所示，包括 100個已知miRNA (61個上調，39個下調)，37個新miRNA (15個上調，22個下調)。呈下調趨勢的miRNA家族包括miR8175、miR5564、miR444、miR166在內的21個家族，呈上調趨勢的miRNA家族包括miR156、 miR169、miR172、miR393、miR397、miR408等10個家族。這些具有差異性表達的miRNA在甘蔗響應低溫脅迫情況下可能發揮特定的功能。

表 2 甘蔗響應低溫差異表達miRNATable 2 Differential expressions of miRNAs in sugarcane responding to low temperature

續表2

續表2

2.4 miRNA作用靶基因預測與分析

根據篩選差異表達的miRNA與對應物種的基因序列信息，使用psRNATarget、TargetFinder和Tapurhybrid軟件進行靶基因預測，結果見圖4。3種預測方法共預測了1 844個相同的潛在靶基因，其中1 696個屬于已知miRNA靶基因，148個屬于新miRNA靶基因。為進一步分析預測所得靶基因的潛在生物學功能，對預測所得的靶基因進行GO分析，共鑒定出甘蔗響應低溫脅迫miRNA靶基因：在生物過程(biological process)有13個功能亞類，在細胞組分(cellular component)有11個功能亞類，在分子功能(molecular function)有7個功能亞類(圖4)。在生物過程類別中，預測靶基因的主要生物學功能富集于細胞過程和代謝過程；在細胞組分類別中，預測靶基因生物學功能主要富集于細胞、細胞器和生物膜；在分子功能類別中，預測靶基因生物學功能主要富集于結合和催化活性及轉運活性。大多靶基因功能均與這些結合功能及其他相近的結合功能相關。大多miRNA通過直接或間接介導靶基因的表達調控相關代謝途徑響應低溫脅迫，這些miRNA所調控的靶基因對甘蔗的耐寒性起關鍵的調控作用(表 3)。

1. 細胞過程; 2. 代謝過程; 3.生物調節; 4. 發育過程; 5. 定位; 6. 響應刺激; 7. 多細胞有機體過程; 8. 繁殖; 9. 繁殖過程; 10. 信號傳導; 11. 細胞成分組織或生物發生; 12. 多有機體過程; 13. 生長; 14. 細胞; 15. 細胞器; 16. 細胞膜; 17. 細胞膜部分; 18. 蛋白復合體; 19. 細胞器部分; 20. 細胞連接; 21. 共質體; 22. 膜包管腔; 23. 超分子復合體; 24. 胞外區; 25. 結合; 26. 轉錄因子活性; 27. 催化活性; 28. 轉運活性; 29. 結構分子活性; 30. 分子功能因子; 31. 抗氧化活性.1. Cellular process; 2. Metabolic process; 3. Biological regulation; 4. Developmental process; 5. Localization; 6. Response to stimulus; 7. Multicellular organismal process; 8. Reproduction; 9. Reproductive process; 10. Signaling; 11. Cellular component organization or biogenesis; 12. Multi-organism process; 13. Growth; 14. Cell; 15. Organelle; 16. Membrane; 17. Membrane part; 18. Protein-containing complex; 19. Organelle part; 20. Cell junction; 21. Symplast; 22. Membrane-enclosed lumen; 23. Supramolecular complex; 24. Extracellular region; 25. Binding; 26. Transcription regulator activity; 27. Catalytic activity; 28. Transporter activity; 29. Structural molecule activity; 30. Molecular function regulator; 31. Antioxidant activity.圖 4 差異表達miRNA靶基因的GO功能分析Fig. 4 GO function analysis of miRNAs target genes for differentially expressed

表 3 響應低溫的部分甘蔗miRNA及其靶基因Table 3 Some miRNAs and their target genes in sugarcane responding to low temperature

2.6 qRT-PCR驗證

利用qRT-PCR技術對測序所得到的miRNA及靶基因的差異豐度進行驗證，miRNA及其靶基因的引物序列如表4所示。在對照組與低溫處理組中共篩選出14個差異表達的miRNA及其靶基因，并對它們進行qRT-PCR驗證，包括miR156、miR160g_1、miR167d、miR169 a-3p_3、miR171b-3p_3、miR172d-5p_4、miR393-3p_1、miR396b、 miR397a_3、miR398b、miR399k_1、miR408d、novel_mir36、novel_mir89(圖5)。根據表4和圖5的結果，除novel-miR36外，選擇的其余13個miRNA在qRT-PCR實驗中的表達模式與高通量測序檢測到的表達模式一致(圖5：A)，這表明大部分高通量測序結果均能被熒光定量 PCR 技術所驗證，此次研究數據為真實可信的。所篩選的14個miRNA，除novel-miR36外，其余13個miRNA均與其靶基因呈負調控關系(圖5：B)，證明miRNA通常以負調控或沉默靶基因的方式來調控植物的生長、發育及對環境的應激反應。

表 4 miRNA和靶基因qRT-PCR引物設計Table 4 Design of qRT-PCR primer for miRNAs and target genes

A. miRNA qRT-PCR 檢驗; B. 靶基因qRT-PCR 檢驗。A. Results of miRNA qRT-PCR; B. Results of target genes.圖 5 miRNA及其靶基因qRT-PCR分析Fig. 5 qRT-PCR analysis of miRNA and target genes

3 討論與結論

甘蔗熱帶種原產熱帶地區，喜溫，現代甘蔗栽培品種以熱帶種的種質為主體。作為廣西傳統的主導產業，而低溫不僅是限制其擴大種植區域和實現高產穩產的重要因素之一，還影響蔗農收益和糖業穩定發展(蘇永秀等, 2006)。為了解甘蔗響應低溫的內在分子機制，挖掘其與耐寒相關的miRNA及相關靶基因，本研究對不同基因型甘蔗進行低溫脅迫處理，通過高通量測序技術及生物信息學方法，系統分析不同抗寒能力的甘蔗對低溫脅迫的響應。

通過高通量測序技術及生物信息學分析發現，sRNA主要集中在21～24 nt之間，不同抗寒能力的材料之間有一定差異且不同長度的測序頻率不同。sRNA在不同物種間的長度分布會有所區別，如擬南芥(Pasquinelli et al., 2000)、小麥(Meng et al., 2013)、棉花(Sripathi et al., 2014)的sRNA長度最多分布在24 nt，楊樹(李明娜等, 2014)、 大豆 (Turner et al., 2012)和番茄 (Pilcheret al., 2007)的sRNA長度最多分布在21 nt，這些結果與本研究結果高度相似。本研究通過生物信息學手段共挖掘137個miRNA在低溫脅迫前后進行差異性表達，其中100個已知miRNA (61個上調，39個下調)，37個新miRNA(15個上調, 22個下調)。在低溫脅迫下植物通過調節miRNA的表達水平，進而調節對應靶基因的表達，從而引起相關代謝與信號轉導途徑的變化來實現對逆境的響應，其信號轉導途徑主要包括胞外信號途徑、胞內第二信使、轉錄因子以及功能基因等。Wu和Poethig (2006)的研究發現，低溫脅迫導致miR156下調而其靶基因表達量增加， 進而調控擬南芥營和養期延長、生長代謝變緩，以應對不良環境，這與本研究低溫脅迫后miR156下調的結果一致。王麗麗等(2017)在擬南芥、小麥和水稻研究中發現miR160響應低溫脅迫。而本研究中，在低溫脅迫后miR160g_1表達下調，靶基因表達受抑制且作用于生長素信號通路從而在抵抗低溫中發揮作用。Pourcel等( 2005)研究表明miR397的靶基因與漆酶有關，與細胞壁木質素的合成、抗病、對環境的適應過程密切相關。而本研究發現，低溫脅迫后甘蔗葉片miR397a_3呈上調表達，其靶基因與抗壞血酸氧化酶(L-ascorbate oxidase)有關，可能參與調控抗壞血酸氧化酶的表達而增強低溫響應能力。Li等(2014)研究旱芹miRNA對低溫脅迫的響應發現， miR160、 miR164、miR394、 miR395、miR408具有表達差異。Sunkar和Zhu (2004)在研究分析擬南芥響應脅迫miRNA中，發現一些miRNA能被多種脅迫因素誘導，如miR393受到低溫、干旱、高鹽的誘導表達。Gupta等(2014)研究小麥miRNA對低溫脅迫、鹽脅迫、滲透脅迫的響應發現，在鹽脅迫和低溫脅迫下，miR168、miR397均表達下調，而miR172表達上調；miR393在滲透和鹽脅迫下表達量上升，在低溫脅迫下表達量下降。Sun等(2015) 研究發現，在低溫脅迫下對葡萄的miR169的表達量上調，但在擬南芥、扁桃中研究發現miR169的表達量下調。結合前人的研究結果，本研究認為特定miRNA對低溫脅迫的響應可能因植物種類、同一植物的不同基因型、不同的組織類型、脅迫的時間等而有所差異。

本研究利用psRNATarget、TargetFinder 及Tapirhybrid三種分析方式預測靶基因，對預測所得的靶基因進行GO分析發現，在生物過程類別中預測靶基因的主要生物學功能富集于細胞過程和代謝過程；在細胞組分類別中預測靶基因生物學功能主要富集于細胞、細胞器和生物膜；在分子功能類別中預測靶基因生物學功能主要富集于結合和催化活性及轉運活性。在逆境脅迫條件下植物通過調節miRNA的表達水平，進而調節對應靶基因的表達，從而引起相關代謝與信號轉導途徑的變化來實現對逆境的響應。因此，進一步對參與調控植物激素信號傳導、光合色素合成、抗氧化酶系統、泛素介導蛋白水解、淀粉與蔗糖代謝等通路的靶基因進行qRT-PCR驗證。本研究發現，miR156靶向調控SBP轉錄因子，低溫脅迫前后，miR156下調表達，負向調控SBP轉錄因子，使甘蔗生長代謝變緩而耐寒能力增強。梅琳(2007)研究冬小麥發現，低溫脅迫下miR160的靶基因ARF17在過表達miR160f擬南芥植株中表達量明顯下降。在擬南芥、水稻、玉米研究中證明miR160、miR167的靶基因是ARF(AuxinResponseFactors), 主要通過調控生長素信號通路應對逆境。而本研究也發現甘蔗在低溫脅迫后miR160g_1上調表達且其作用于轉錄因子ARF。張譯云等(2012)研究發現毛白楊在低溫脅迫下miR169ac表達下調，對其靶基因轉錄因子NAC負向調控。黨春艷(2013)在研究高山離子芥低溫脅迫響應中發現，miR169a在冷脅迫下無顯著表達，而miR169在冷脅迫下表達下調，表明miR169家族之間響應低溫存在差異。本研究中，miR169 a-3p_3在低溫脅迫后表達上調，作用于細胞核轉錄因子使其表達下調，推測miR169 a-3p_3與調控甘蔗低溫耐受能力密切相關。

此外，本研究通過生物信息學挖掘甘蔗響應低溫的miRNA，不僅參與抗氧化酶系統、植物激素信號傳導、遺傳信息等，而且參與類胡蘿卜素代謝、卟啉與葉綠素代謝、淀粉與蔗糖代謝途徑相關基因的表達。有研究表明葉綠體內的類胡蘿卜素既可作為吸收光能的輔助色素將能量傳遞給葉綠素a，又可在持續脅迫條件下將過剩光能安全耗散而保護光合機構(Liu et al., 2004; Holt et al., 2005)。我們前期的研究也發現，低溫脅迫時甘蔗能夠通過增加光合機構過剩激發能的耗散和調整其葉片光合色素含量與構成進行有效的光能利用和分配。本研究對部分差異表達miRNA及其靶基因進行qRT-PCR檢驗發現，參與類胡蘿卜素代謝的miR167d對靶基因具有負向調控作用。因此，后續我們將對低溫條件下參與光合生理過程的miRNA分子調控機制進行深入研究，為抗寒育種提供重要理論依據和實驗基礎。