小麥響應干旱脅迫環(huán)狀RNA的鑒定

李寧，柳坤，劉彤彤，史雨剛，王曙光，楊進文，孫黛珍

小麥響應干旱脅迫環(huán)狀RNA的鑒定

李寧，柳坤，劉彤彤，史雨剛，王曙光，楊進文，孫黛珍

山西農(nóng)業(yè)大學農(nóng)學院，山西太谷 030801

【目的】干旱是限制全球小麥生產(chǎn)的主要非生物脅迫之一，探索小麥應對干旱的分子調(diào)控機制對小麥分子育種具有重要意義。環(huán)狀RNA（circRNA）已被證實在植物抵御外界環(huán)境脅迫的過程中扮演著重要角色。鑒定小麥響應干旱脅迫的circRNA，有助于構(gòu)建小麥干旱脅迫響應的調(diào)控網(wǎng)絡，為解析小麥的抗旱性機制奠定基礎。【方法】以2個抗旱性差異顯著的小麥品種（周麥13和冀麥38）為試驗材料，對其在干旱及對照條件下的根部樣本進行circRNA測序。鑒定小麥circRNA并對其進行特征分析，篩選與干旱脅迫響應相關的差異表達circRNA，并對其靶向microRNA（miRNA）進行預測。進一步根據(jù)miRNA及其靶基因在干旱脅迫下的表達模式，構(gòu)建小麥響應干旱脅迫的潛在circRNA-miRNA-mRNA調(diào)控模塊。【結(jié)果】共鑒定獲得1 409個小麥circRNA，其中，多數(shù)（68.91%）為外顯子circRNA，且僅有133個circRNA在2個品種中被同時鑒定獲得。在干旱脅迫下共鑒定獲得239個差異表達circRNA，其中138個circRNA在抗旱型品種周麥13（ZM13）中特異性差異表達，19個circRNA在2個品種中同時差異表達。共預測到34個靶向miRNA以及1 408個miRNA靶基因。根據(jù)這些差異表達circRNA、靶向miRNA以及miRNA靶基因在干旱脅迫后的表達模式，共篩選出5個分別以tae-miR9664-3p、tae-miR1122b-3p、tae-miR9662a-3p、tae-miR6197-5p和tae-miR1120c-5p為中心的小麥響應干旱脅迫的潛在circRNA-miRNA-mRNA調(diào)控模塊。【結(jié)論】小麥circRNA具有明顯的品種特異性，且在不同抗旱型小麥品種之間具有不同的表達模式。共鑒定到239個響應干旱脅迫的小麥circRNA以及5個潛在的circRNA-miRNA-mRNA調(diào)控模塊。

小麥；circRNA；干旱脅迫；miRNA

0 引言

【研究意義】干旱是全球糧食生產(chǎn)所面臨的嚴重自然災害之一[1]，提高作物的抗旱能力是廣大育種工作者的主要研究方向。小麥（L.）是世界三大糧食作物之一，挖掘其響應干旱脅迫的調(diào)控網(wǎng)絡并分析其遺傳基礎，對小麥耐旱性分子設計育種和耐旱性機制解析具有重要意義。【前人研究進展】環(huán)狀RNA（circular RNA，circRNA）是一類閉合的非線性RNA分子，主要由pre-mRNA通過可變剪切加工產(chǎn)生，大量存在于真核細胞轉(zhuǎn)錄組中[2]。由于高通量測序和高效生物信息學分析技術的不斷突破，研究人員已經(jīng)在小鼠、人類、古細菌和秀麗隱桿線蟲等物種中相繼鑒定出了大量的circRNA[3]。對其特征分析結(jié)果表明，circRNA可以來自外顯子、內(nèi)含子和基因間區(qū)域，且其表達模式通常在不同組織和發(fā)育階段存在特異性[4]。此外，有證據(jù)表明circRNA比線性RNA要更加穩(wěn)定，因此，它們可能通過一些特殊的方式參與生命活動過程[4]。近年來，circRNA在動物中研究最多的功能之一是它可以與microRNA（miRNA）競爭性結(jié)合，來調(diào)節(jié)其靶基因的功能。如，Zheng等[5]研究表明，人類中的有18個潛在的miRNA結(jié)合位點，可以與9個miRNA相結(jié)合。此外，Wang等[6]研究證實，與心臟相關的circRNA（）可作為競爭性內(nèi)源RNA（ceRNA）與miR-223相結(jié)合并抑制其活性。此外，在人類中的研究證實，一些circRNA可以通過特定的相互作用關系來促進其宿主基因在順式作用中的表達[7]。與動物相比，人們對植物中circRNA鑒定的研究相對較少。2014年，Wang等[8]首次在植物中進行了circRNA鑒定，并且發(fā)現(xiàn)了擬南芥中存在circRNA的明確證據(jù)。之后，Ye等[9]通過全基因組水平的檢測，在水稻（）和擬南芥中分別鑒定出12 037和6 012個circRNA，證實circRNA在單子葉植物和雙子葉植物中均廣泛存在。后來，研究人員相繼在番茄[10]、馬鈴薯[11]、小麥[12]、大豆[13]、茶[14]、大麥[15]和玉米[16]中鑒定到circRNA。盡管人們已經(jīng)在植物中鑒定到大量的circRNA，但迄今為止，關于circRNA在植物當中的功能研究十分有限。Tan等[17]在番茄中過表達一個來自類胡蘿卜素合成關鍵基因（）的circRNA，可導致番茄紅素合成酶基因的表達降低以及番茄紅素和β-胡蘿卜素的積累減少。Cheng等[18]過量表達了一個來自AT5G37720第一個內(nèi)含子的circRNA，可以改變約800個基因的表達并最終影響了擬南芥的發(fā)育。此外，越來越多研究表明，circRNA在調(diào)節(jié)植物對環(huán)境脅迫的響應方面發(fā)揮著重要作用。在水稻中鑒定到27個在磷酸鹽充足和饑餓條件下差異表達的外顯子circRNA[9]。Pan等[19]在熱脅迫條件下的擬南芥中鑒定到的circRNA數(shù)量比正常條件下要多，推測其可能通過由circRNA介導的ceRNA網(wǎng)絡參與植物對熱脅迫的反應。Wang等[12]在小麥葉片中鑒定到88個circRNA，其中66個在脫水脅迫后差異表達，并且從這些差異表達circRNA中預測到26個靶向miRNA。【本研究切入點】盡管在小麥中已進行了相關circRNA的鑒定，但鑒于circRNA在不同品種間的特異性，仍需對其進行大量的鑒定工作來完善小麥circRNA信息。此外，本研究在鑒定小麥circRNA的基礎上，充分利用本研究組前期在相同小麥品種中已獲得的miRNA及mRNA的轉(zhuǎn)錄組數(shù)據(jù)，快速構(gòu)建了小麥響應干旱脅迫的circRNA-miRNA-mRNA調(diào)控模塊。【擬解決的關鍵問題】本研究利用2個不同抗旱型小麥品種為試驗材料，對小麥circRNA進行鑒定和特征分析，并進一步鑒定在干旱脅迫后差異表達的circRNA，篩選小麥干旱脅迫響應相關的潛在circRNA-miRNA-mRNA模塊，為探索小麥響應干旱脅迫的調(diào)控網(wǎng)絡及其circRNA的作用機制奠定基礎。

1 材料與方法

1.1 試驗材料和干旱處理方法

2個小麥（L.）品種分別為周麥13（ZM13）和冀麥38（JM38）。ZM13是一個抗旱型品種，JM38是一個干旱敏感型品種[20]。小麥種子經(jīng)75%乙醇溶液消毒2 min，再用2% H2O2消毒30 min，無菌水洗滌3次[21]。室溫條件下，將種子置于蒸餾水中發(fā)芽3 d，然后轉(zhuǎn)移至96孔塑料水培盒（13 cm×8.5 cm×11 cm），光照培養(yǎng)箱（22℃光照14 h/20℃黑暗10 h）中培養(yǎng)[22]。待幼苗培養(yǎng)7 d（一葉一心期）后，對其進行干旱處理，干旱溶液為20% PEG-6000的Hoagland營養(yǎng)液，設置對照，3次試驗重復。

1.2 總RNA提取和circRNA測序

干旱脅迫4 d時，收集2個品種CK和干旱處理的幼苗根部，每個樣本分別由3株幼苗全部根系混合而成，每個品種每個處理3次重復，共計12個樣本，液氮冷凍，-80℃保存。按照Yin等[10]方法進行總RNA的提取和濃度的測定。利用Ribo-Zero rRNA Removal Kit（Epicentre, Madison, WI, USA）去除總RNA中的核糖體RNA，利用RNase R（Epicentre, Madison, WI, USA）去除線性RNA，構(gòu)建circRNA-seq文庫[10]。委托北京百邁客生物科技有限公司利用Illumina nova-seq 6000平臺對文庫進行雙末端測序。

1.3 circRNA的鑒定和表達分析

首先去除原始數(shù)據(jù)（raw data）中未知（N）堿基超過5%、包含接頭序列以及低質(zhì)量堿基（Q≤20）超過50%的低質(zhì)量read。使用HISAT軟件[23]將高質(zhì)量的read（clean read）與小麥參考基因組（IWGSC_ RefSeq_v1.1.Triticum_aestivum.）進行序列比對。利用find_circ軟件[24]進行circRNA鑒定。參照Zhu等[25]方法，將每個circRNA在樣本中的表達水平進行標準化。利用DESeq2[26]軟件進行circRNA的差異表達分析，差異表達circRNA的篩選標準為：|log2(Fold change)|＞1且＜0.05。

1.4 circRNA-miRNA靶向關系預測、miRNA靶基因預測和基因的功能注釋

使用Target Finder軟件[27]預測差異表達circRNA的miRNA結(jié)合位點并對miRNA的靶基因進行預測。使用KEGG（Kyoto Encyclopedia of Genes and Genomes）數(shù)據(jù)庫對差異表達circRNA的宿主基因進行KEGG注釋。利用小麥多組學網(wǎng)站W(wǎng)heatOmics 1.0（WheatOmics sdau.edu.cn）[28]鑒定宿主基因在水稻中的同源基因，并利用國家水稻數(shù)據(jù)中心（www.ricedata. com）查詢同源基因的已知功能。使用BLAST2GO軟件[29]對宿主基因進行GO注釋。使用Swiss-Prot數(shù)據(jù)庫對miRNA的靶基因進行功能注釋。

1.5 circRNA-miRNA-mRNA調(diào)控模塊預測方法

利用本研究組已鑒定到的JM38和ZM13在干旱脅迫下的miRNA數(shù)據(jù)[20]（NCBI，SRA數(shù)據(jù)庫，BioProject PRJNA837867）分析本研究中差異表達circRNA的靶向miRNA在干旱脅迫下的表達模式，從中篩選出差異表達的靶向miRNA。之后同樣利用本研究組已鑒定到的JM38和ZM13在干旱脅迫下的mRNA數(shù)據(jù)[20]（NCBI，SRA數(shù)據(jù)庫，BioProject PRJNA838787）分析本研究中靶向miRNA的靶基因在干旱脅迫下的表達模式，從中篩選出差異表達的靶基因。然后根據(jù)靶基因的注釋功能及同源基因的已知功能對其進行篩選。最終構(gòu)建以差異表達circRNA為基礎，與其存在靶向關系且差異表達的miRNA為中心，以及根據(jù)注釋功能選擇的差異表達靶基因為目標的circRNA-miRNA-mRNA調(diào)控模塊，作為小麥響應干旱脅迫的潛在circRNA-miRNA- mRNA調(diào)控模塊。

1.6 差異表達circRNA的驗證

利用實時熒光定量PCR（qRT-PCR）方法對隨機選擇的9個差異表達circRNA的表達模式進行驗證。按照Trizol試劑說明書提取總RNA，按照RevertAid First Strand cDNA Synthesis Kit（Thermo Scientific, USA）反轉(zhuǎn)錄試劑盒說明書進行cDNA的合成。使用SYBR Green Master Mix在CFX 96實時PCR系統(tǒng)（Bio-Rad）上進行qRT-PCR分析。采用2?ΔΔCt方法計算每個circRNA的相對表達量。以小麥管家基因作為內(nèi)參基因[12]，所用引物的信息見表1。

2 結(jié)果

2.1 小麥中circ RNA的鑒定

為鑒定小麥中的circRNA，提取這兩個品種在對照及干旱脅迫4 d后的根部樣本進行測序。結(jié)果顯示，從12個樣本中共獲得約210.36 Gb的高質(zhì)量數(shù)據(jù)（clean data），單個樣品的clean data在14.50 Gb以上，且每個樣本高質(zhì)量read（clean read）的Q30均能達到93.01%以上，進一步說明該測序結(jié)果的可靠性。與小麥參考基因組的比對結(jié)果顯示，各樣品clean read的比對效率從98.74%到99.79%不等（表2）。

表1 熒光定量PCR引物序列

表2 circRNA測序數(shù)據(jù)結(jié)果統(tǒng)計

a：CK：對照條件下的樣本；T：干旱條件下的樣本；1—3分別代表3個不同重復

a: CK: samples under control conditions; T: samples under drought conditions; 1-3 represent three different repetitions

利用circRNA預測工具find_circ在12個測序樣本中共鑒定到1 409個circRNA，其中，在干旱敏感品種JM38中鑒定到722個circRNA，在抗旱品種ZM13中鑒定到820個circRNA（圖1-A）。此外，在1 409個circRNA中，有589和687個circRNA分別在JM38和ZM13中特異性表達，僅有133個circRNA在2個品種中均被鑒定到（圖1-A），表明circRNA在小麥中具有較高的品種特異性。

2.2 小麥circRNA的特征分析

基因組起源分析顯示，在1 409個已鑒定到的circRNA中，多數(shù)（971個，68.91%）是外顯子circRNA，其次，有299個（21.22%）是基因間circRNA，其余139個（9.87%）為內(nèi)含子circRNA（圖1-B）。在染色體分布方面，circRNA主要分布在小麥5B、1B、2A、2B以及3B等染色體上（圖1-C）。關于circRNA長度的變異情況，大多數(shù)circRNA（1 079個，76.58%）的長度低于800 bp（圖1-D）。此外，大部分（78.69%）長度大于2 600 bp的長circRNA屬于基因間circRNA（圖1-D）。

A：2個品種中circRNA數(shù)量的韋恩圖；B：由外顯子、內(nèi)含子和基因間區(qū)域產(chǎn)生的circRNA的數(shù)量和百分比；C：circRNA在小麥染色體上的分布；D：小麥circRNA的長度分布

2.3 干旱脅迫誘導的差異表達circRNA

為鑒定響應干旱脅迫的小麥circRNA，比較了2個品種的對照和干旱處理之間circRNA的表達譜。在所有鑒定到的1 409個circRNA中，有239個circRNA在干旱脅迫下存在顯著差異表達，其中，在干旱敏感型品種JM38中共鑒定到91個差異表達的circRNA，在抗旱品種ZM13中共鑒定到157個差異表達的circRNA，其中138個circRNA在ZM13中特異性差異表達，19個circRNA在2個品種中同時差異表達（圖2-A）。此外，值得注意的是，在JM38中更多的circRNA下調(diào)表達（55個，60.44%），而在ZM13中，有88%以上（139個，88.54%）的circRNA均為上調(diào)表達（圖2-B）。說明circRNA在不同抗旱型小麥品種中存在不同的表達模式。

為驗證circRNA測序結(jié)果的準確性和可靠性，利用qRT-PCR方法從239個差異表達circRNA中隨機選擇9個circRNA進行差異表達驗證。結(jié)果顯示，除了2個circRNA（chr3D:456098361|456098885和chr3B:438537516|438537814）以外，其余7個circRNA在2種試驗方法中表現(xiàn)出相似的表達模式（圖3），這可能是由于其絕對表達量較低或2種不同的表達模式計算方法所造成的結(jié)果。

2.4 差異表達circRNA宿主基因的生物學功能分析

circRNA可通過對其宿主基因的順式調(diào)節(jié)在轉(zhuǎn)錄控制中發(fā)揮重要作用[7]。在2個品種中差異表達的239個circRNA中共鑒定出182個宿主基因。為了解這些差異表達circRNA的生物學功能，首先對宿主基因進行了GO注釋分析。結(jié)果顯示，共有114個宿主基因被注釋到30個不同的GO分類上（圖4-A）。分別按照生物過程、細胞成分和分子功能對這些GO類別進行分類。注釋到較多宿主基因的生物過程有：代謝過程（GO：0008152）、細胞過程（GO：0009987）、單有機體過程（GO：0044699）和刺激響應（GO：0050896）。注釋到較多宿主基因的細胞成分類別有：膜（GO：0016020）、細胞（GO：0005623）、膜組分（GO：0044425）、細胞器（GO：0043226）和大分子復合物（GO：0032991）。注釋到較多宿主基因的分子功能類別有：結(jié)合（GO：0005488）、催化活性（GO：0003824）、核酸結(jié)合轉(zhuǎn)錄因子活性（GO：0001071）、轉(zhuǎn)運蛋白活性（GO：0005215）和抗氧化活性（GO：0016209）（圖4-A）。

A：2個品種中差異表達circRNA數(shù)量的韋恩圖；B：2個品種中上調(diào)表達circRNA和下調(diào)表達circRNA的數(shù)量；C：差異表達circRNA在2個品種中的表達熱圖，數(shù)字代表log2(fold change)值

其次，分析了在抗旱品種ZM13中特異性差異表達circRNA的宿主基因的GO注釋結(jié)果，共鑒定出110個宿主基因，其中，共有69個宿主基因被注釋到28個不同的GO類別上。和所有差異表達circRNA宿主基因的GO注釋結(jié)果相似，注釋到較多宿主基因的GO類別有：結(jié)合、催化活性、代謝過程、細胞過程、膜、細胞和膜組分等（圖4-B）。

此外，對2個品種中共同差異表達的19個circRNA的宿主基因進行GO注釋分析，共鑒定出14個宿主基因，其中有8個宿主基因被注釋到20個不同的GO類別上（表3）。除了上文中提到的主要GO類別以外，注釋到宿主基因數(shù)量在2個及2個以上的GO類別還有生物調(diào)節(jié)（GO：0065007）以及細胞外區(qū)域（GO：0005576）（表3）。

通過利用KEGG數(shù)據(jù)庫對所有的差異表達circRNA宿主基因的功能進行進一步注釋。結(jié)果顯示，共有67個宿主基因被注釋到36個不同的KEGG通路上。包括植物-病原菌互作（ko04626）、內(nèi)吞作用（ko04144）、真核生物中的核糖體生物發(fā)生（ko03008）、植物激素信號轉(zhuǎn)導（ko04075）、苯丙烷生物合成（ko00940）和過氧化物酶體（ko00190）等（表4）。通過對ZM13中特異性差異表達circRNA的宿主基因進行了KEGG注釋，結(jié)果顯示，共有38個宿主基因被注釋到86個不同的KEGG通路上（表5）。其中，富集到較多宿主基因的通路包括：植物-病原菌互作、內(nèi)吞作用和過氧化物酶體等。此外同樣對2個品種中共同差異表達circRNA的宿主基因進行KEGG注釋。結(jié)果顯示，共有7個宿主基因被注釋到7個不同的KEGG通路上（表6）。

圖3 差異表達circRNA的表達模式驗證

表3 2個品種中共同差異表達circRNA的宿主基因的GO注釋

僅列出基因數(shù)量大于等于2的GO term Only GO terms with a No. of genes greater than or equal to 2 are listed

表4 全部差異表達circRNA宿主基因的KEGG注釋

僅列出基因數(shù)量大于等于2的通路 Only pathways with a No. of genes greater than or equal to 2 are listed

A：兩品種中所有差異表達circRNA宿主基因的GO注釋結(jié)果；B：ZM13中特異性差異表達circRNA宿主基因的GO注釋結(jié)果

A: GO annotation of host genes of all differentially expressed circRNAs; B: GO annotation of host genes of differentially expressed circRNAs specifically in ZM13

圖4 差異表達circRNA宿主基因的GO注釋

Fig. 4 GO annotation of host genes of differentially expressed circRNAs

此外，對所有宿主基因在水稻中的同源基因進行比對。通過查詢這些同源基因的生物學功能發(fā)現(xiàn)，其中，來自9個宿主基因的10個水稻同源基因已被證實在水稻抵御多種非生物脅迫的過程中發(fā)揮著重要的作用（表7）。

2.5 差異表達circRNA的靶向miRNA預測

circRNA可以通過與miRNA競爭性結(jié)合，來阻止miRNA調(diào)節(jié)其靶基因的表達[6]。通過對circRNA的靶向miRNA進行預測，發(fā)現(xiàn)共有24個差異表達circRNA（其中，12個circRNA在ZM13中特異性差異表達，3個circRNA在2個品種中同時差異表達）預測到miRNA結(jié)合位點，并鑒定出34個不同的靶向miRNA（表8）。在所有差異表達的circRNA中，chr1A: 589450964|589476548具有最多的miRNA結(jié)合位點，為7個；其次是chr5B:528574134|528588868，具有6個miRNA結(jié)合位點；chr3B:778404617|778442815、chr6D: 454657953|454748495、chr7B:1308928|1352258和chr4D: 65120787|65147651均各具有5個miRNA結(jié)合位點（表8）。此外，還發(fā)現(xiàn)miRNA的靶向circRNA也不是唯一的，也就是說多個circRNA可能會具有同一個靶向miRNA。如tae-miR1117和tae-miR1133都具有7個靶向circRNA。在2個品種中同時差異表達的3個circRNA（chr5A:482358772|482427161、chr2B:787486469| 787538256和chr4A:668278340|668323672）共預測到5個靶向miRNA（表8）。

表5 ZM13中特異性差異表達circRNA的宿主基因的KEGG注釋

僅列出基因數(shù)量大于等于2的通路 Only pathways with a No. of genes greater than or equal to 2 are listed

表6 2個品種中共同差異表達circRNA的宿主基因的KEGG注釋

表7 差異表達circRNA宿主基因在水稻中的同源基因及其抗逆性功能

表8 差異表達circRNA及其靶向miRNA

circRNA ID斜體代表其在ZM13中特異性差異表達，circRNA ID加粗代表在2個品種中同時差異表達

CircRNA IDs in italics represent their specific differential expression in ZM13, and circRNA IDs in bold represent differentially expressed in both varieties

2.6 circRNA-miRNA-mRNA調(diào)控模塊預測

為了探索小麥響應干旱脅迫的circRNA-miRNA- mRNA調(diào)控網(wǎng)絡。首先利用本研究組已鑒定到的小麥干旱脅迫響應相關的miRNA數(shù)據(jù)分析34個靶向miRNA在干旱脅迫下的表達情況。結(jié)果表明，共有6個miRNA在干旱脅迫后差異表達，且均為上調(diào)表達（圖5-A—B）。其中tae-miR9664-3p和tae-miR7757-5p在2個小麥品種均上調(diào)表達（圖5-C）。然后利用Target Finder軟件對這6個miRNA的靶基因進行預測，共預測到1 408個靶基因。之后，同樣利用已鑒定到的小麥干旱脅迫響應相關的mRNA數(shù)據(jù)分析這1 408個靶基因干旱脅迫下的表達情況。結(jié)果表明，共有261個靶基因在干旱脅迫后差異表達（圖5-D—F）。進一步查詢這261個基因的注釋功能，從中篩選出32個可能與小麥脅迫響應相關候選靶基因。它們的注釋功能有：過氧化物酶、轉(zhuǎn)錄因子、細胞色素P450、鈣依賴性蛋白激酶、絲氨酸/蘇氨酸蛋白激酶、生長素反應蛋白和水通道蛋白等（表9）。

綜上所述，最終構(gòu)建了5個小麥干旱脅迫響應相關的潛在circRNA-miRNA-mRNA調(diào)控模塊，其中，包括8個差異表達circRNA、5個差異表達的靶向miRNA以及32個差異表達的miRNA靶基因（表9）。其中模塊3中的circRNA（chr5D:543694782| 543696753）以及其的靶向miRNA（tae-miR9662a- 3p）均在抗旱型品種ZM13中特異性差異表達（表9）。

A、B：差異表達的靶向miRNA在2個品種中的分布；C：差異表達的靶向miRNA在2個品種中的表達模式；D、E：差異表達靶基因在2個品種中的分布；F：差異表達靶基因在2個品種中的表達模式

3 討論

3.1 小麥circRNA的鑒定與特征分析

circRNA是由下游剪接供體和上游剪接受體之間的非線性反向剪接事件所產(chǎn)生一類特殊的RNA。近年來，隨著circRNA-seq技術的發(fā)展，研究人員在動物和植物中鑒定到了大量的circRNA。這些廣泛表達且高度保守的circRNA的發(fā)現(xiàn)增加了人們對于非編碼RNA潛在功能的認識[40]。然而，與動物相比，人們對植物中circRNA的生物發(fā)生和調(diào)控功能在很大程度上都是未知的[11]。Ye等[9]對水稻和擬南芥中的circRNA進行了特征分析，并且發(fā)現(xiàn)植物circRNA的生物發(fā)生機制可能與動物中的有所不同。此外，最近在水稻[9]、小麥[12]和馬鈴薯[11]上的研究表明，circRNA在調(diào)節(jié)植株對外界脅迫的反應中發(fā)揮著重要的作用。

表9 circRNA-miRNA-mRNA調(diào)控模塊

本研究分別對2個在抗旱性方面存在顯著差異的小麥品種進行了circRNA測序，共鑒定到1 409個circRNA，其中，在干旱敏感品種JM38中鑒定到722個circRNA，在抗旱品種ZM13中鑒定到820個circRNA（圖1-A）。通過對比2個品種中鑒定到的circRNA，發(fā)現(xiàn)僅有133個（9.44%）circRNA在2個品種中被同時鑒定到，說明circRNA在小麥中具有明顯的品種特異性，這與前人研究的結(jié)果一致。本研究在小麥中鑒定到的circRNA數(shù)量要少于前人在水稻（26 160個）[9]、擬南芥（31 079個）[9]和大豆（5 372個）[13]中的鑒定結(jié)果，但要高于前人同樣在小麥（88個）[12]以及玉米（496個）[16]、茶（342個）[14]和番茄（854個）[10]中的鑒定結(jié)果，這可能是由于不同的研究中在物種、組織（例如葉、根、果實、莖、芽等）以及circRNA預測工具（例如CIRI2、fnd_circ、CIRCexplorer）上存在的差異所導致的結(jié)果[41]。此外，前人研究結(jié)果表明，circRNA可以在外顯子、內(nèi)含子和基因之間的區(qū)域產(chǎn)生，在水稻、擬南芥、茶和番茄的研究中，大多數(shù)已鑒定的circRNA是位于外顯子的circRNA[9, 14, 42]。本研究鑒定到的大多數(shù)小麥circRNA（971個，68.91%）同樣是外顯子circRNA（圖1-B）。

3.2 小麥干旱脅迫響應相關的circRNA

植物對干旱脅迫的反應是一個非常復雜的過程，其中，涉及到許多干旱誘導基因和信號轉(zhuǎn)導途徑的參與。為了探索小麥circRNA在干旱脅迫響應中的作用，本研究共鑒定到239個在干旱脅迫下存在顯著差異表達的circRNA。在這些差異表達的circRNA中，有超過一半（138個，57.74%）的circRNA僅來自于抗旱型品種ZM13中，還有19個circRNA同時來自于2個品種，只有72個（30.13%）circRNA在干旱敏感型品種JM38中特異性差異表達（圖2-A）。這可能是由于circRNA在小麥中的品種特異性所導致的結(jié)果。此外，值得的注意的是，在JM38中更多的circRNA下調(diào)表達（55個，60.44%），而在ZM13中，有88%以上（139個，88.54%）的circRNA均為上調(diào)表達（圖2-B）。這說明circRNA在不同抗旱型小麥品種中存在不同的表達模式。

本研究通過對這些差異表達circRNA的宿主基因進行了GO和KEGG注釋分析。結(jié)果表明，包括刺激響應、細胞膜組分、轉(zhuǎn)運蛋白活性和抗氧化活性等一些與干旱脅迫反應相關的GO term都注釋到了較多的宿主基因（圖4）。此外，部分宿主基因也注釋到了一些重要的脅迫應答通路：植物激素信號轉(zhuǎn)導、苯丙烷生物合成和過氧化物酶體等（表4）。

檢索circRNA宿主基因在其他作物中的同源基因的功能，可以為預測這些差異表達circRNA的生物學功能提供一些幫助。通過查詢這些差異表達circRNA的宿主基因在水稻中的同源基因及其功能發(fā)現(xiàn)，其中，部分水稻同源基因已被證實在抵御多種非生物脅迫的過程中發(fā)揮著重要的作用（表7）。如chr7A:222784992| 222785194的宿主基因在水稻中的同源基因和可以提高植株對干旱以及鹽脅迫的耐受性[38-39]。差異表達circRNA chr2A:52664278|52664553、chr3A:638334965|638335214和chr6A:608632503|608632694宿主基因的同源基因均可提高水稻的耐鹽性[30, 32, 37]。此外，其余同源基因在植株抵御高溫、冷和營養(yǎng)脅迫等過程扮演著重要的角色[31-36]（表7）。這些circRNA表達模式及功能預測結(jié)果，為挖掘小麥干旱脅迫響應相關的circRNA提供了重要依據(jù)。

3.3 小麥干旱脅迫響應相關的潛在circRNA-miRNA- mRNA調(diào)控模塊

研究表明，circRNA可以通過與miRNA競爭性相合，來降低miRNA對靶基因的沉默效應，從而調(diào)節(jié)miRNA靶基因的表達[43]。為了探索與小麥干旱脅迫相關的circRNA-miRNA-mRNA調(diào)控網(wǎng)絡。本研究構(gòu)建了5個可能與小麥干旱脅迫響應相關的circRNA- miRNA-mRNA調(diào)控模塊（表9）。通過分析這些差異表達靶基因注釋功能，可以對于了解這些circRNA- miRNA-mRNA調(diào)控模塊的潛在功能提供一定的幫助。以tae-miR9664-3p為中心的第1個circRNA-miRNA- mRNA調(diào)控模塊中共篩選出10個差異表達的靶基因，其中的注釋功能為過氧化物酶；的注釋功能為乙烯反應轉(zhuǎn)錄因子；此外還包括3個細胞色素P450基因，3個鈣依賴性蛋白激酶和2個ABC轉(zhuǎn)運體A家族成員（表9）。眾所周知，過氧化物酶可通過催化H2O2的氧化還原來調(diào)節(jié)生物體內(nèi)活性氧的平衡，在植物生長發(fā)育和抵御非生物脅迫中發(fā)揮著重要作用[44]。乙烯反應轉(zhuǎn)錄因子與細胞發(fā)育、激素和逆境信號傳遞有關，在植物抗逆信號轉(zhuǎn)導中具有重要的調(diào)控作用[45]。細胞色素P450同樣被證實參與了植物多種生理代謝途徑，其次生代謝產(chǎn)物在植物體信號傳導以及應對生物和非生物脅迫中發(fā)揮著重要作用[46]。鈣依賴性蛋白激酶是植物中重要的一類Ca2+傳感器，在植物發(fā)育及多種脅迫響應中均扮演重要角色[47-48]。以tae-miR1122b-3p為中心的第2個circRNA-miRNA-mRNA調(diào)控模塊中共篩選出9個差異表達的靶基因，的注釋功能為絲氨酸/蘇氨酸蛋白激酶，同樣注釋為乙烯反應轉(zhuǎn)錄因子；為生長素反應蛋白，是一個水通道蛋白（表9）。植物水通道蛋白是介導植物體內(nèi)水分跨膜運輸?shù)闹饕ǖ溃鋸V泛存在于植物細胞膜上，具有調(diào)節(jié)植物開花、氣孔運動、生殖生長、響應多種逆境脅迫等功能[49]。此外，在水稻中的同源基因為，Cao等[50]研究證實，過量表達可以提高煙草對鹽及干旱脅迫的抵抗能力。在水稻中的同源基因也被證實可以影響水稻的籽粒產(chǎn)量、耐鹽性、根部滲透系數(shù)及種子萌發(fā)率[51]。此外，值得注意的是，以tae-miR9662a-3p為中心的模塊3中的circRNA以及其靶向miRNA均只在抗旱型品種ZM13中特異性差異表達，其靶基因包括3個鋅指蛋白和1個富含半胱氨酸的受體蛋白激酶（表9）。該模塊對ZM13的抗旱性研究提供了一定的研究基礎。

除了上述基因外，這5個模塊中所涉及到的miRNA靶基因還包括1個富含半胱氨酸的受體蛋白激酶、3個轉(zhuǎn)錄因子、3個絲氨酸/蘇氨酸蛋白激酶、2個肉桂酰輔酶A還原酶和1個普遍脅迫蛋白（表9）。這些基因均可作為小麥脅迫響應相關的候選靶基因。以上研究結(jié)果為探索小麥響應干旱脅迫的circRNA- miRNA-mRNA調(diào)控網(wǎng)絡奠定了基礎。

4 結(jié)論

共鑒定到1 409個小麥circRNA，其中多數(shù)為外顯子circRNA，且表現(xiàn)出明顯的品種特異性。共鑒定到239個受干旱脅迫誘導的差異表達circRNA，并預測到34個與之相關的靶向miRNA。根據(jù)這些差異表達circRNA、靶向miRNA以及miRNA靶基因的表達模式，共篩選出5個與小麥干旱脅迫相關的潛在circRNA-miRNA-mRNA調(diào)控模塊。

[1] Akram N A, Waseem M, Ameen R, Ashraf M. Trehalose pretreatment induces drought tolerance in radish (L.) plants: some key physio-biochemical traits. Acta Physiologiae Plantarum, 2016, 38(1): 3.

[2] Jeck W R, Sharpless N E. Detecting and characterizing circular RNAs. Nature biotechnology, 2014, 32(5): 453-461.

[3] Salzman J. Circular RNA expression: its potential regulation and function. Trends in Genetics, 2016, 32(5): 309-316.

[4] Li L, Guo J, Chen Y, Chang C, Xu C. Comprehensive CircRNA expression profile and selection of key circRNAs during priming phase of rat liver regeneration. BMC Genomics, 2017, 18(1): 80.

[5] Zheng Q, Bao C, Guo W, Li S, Chen J, Chen B, Luo Y, Lyu D, Li Y, Shi G. Circular RNA profiling reveals an abundant circHIPK3 that regulates cell growth by sponging multiple miRNAs. Nature Communications, 2016, 7: 11215.

[6] Wang K, Bo L, Fang L, Wang J X, Liu C Y, Bing Z, Zhou L Y, Teng S, Man W, Tao Y.A circular RNA protects the heart from pathological hypertrophy and heart failure by targeting miR-223. European Heart Journal, 2016, 37(33): 2602-2611.

[7] Li Z, Huang C, Bao C, Chen L, Lin M, Wang X, Zhong G, Yu B, Hu W, Dai L. Exon-intron circular RNAs regulate transcription in the nucleus. Nature Structural and Molecular Biology, 2015, 22: 256.

[8] Wang P L, Bao Y, Yee M C, Barrett S P, Hogan G J, Olsen M N, Dinneny J R, Brown P O, Salzman J. Circular RNA is expressed across the eukaryotic tree of life. PLoS One, 2014, 9(3): e90859.

[9] Ye C Y, Chen L, Liu C, Zhu Q H, Fan L J.Widespread noncoding circular RNAs in plants. New Phytologist, 2015, 208(1): 88-95.

[10] Yin J L, Liu M Y, Ma D F, Wu J W, Li S L, Zhu Y X, Han B.Identification of circular RNAs and their targets during tomato fruit ripening. Postharvest Biology and Technology, 2018, 136: 90-98.

[11] Zhou R, Zhu Y X, Zhao J, Fang Z W, Wang S P, Yin J L, Chu Z H, Ma D F.Transcriptome-wide identification and characterization of potato circularRNAs in response to Pectobacterium carotovorumsubspecies Brasilienseinfection.International Journal of Molecular Sciences, 2018, 19(1): 71.

[12] Wang Y, Yang M, Wei S, Qin F, Zhao H, Suo B. Identification of circular RNAsand their targets in leaves of Triticum aestivumL. under dehydration stress.Frontiers in Plant Science, 2017, 7: 2024.

[13] Zhao W, Cheng Y, Zhang C, You Q, Shen X, Guo W, Jiao Y.Genome-wideidentification and characterization of circular RNAs by high throughputsequencing in soybean.Scientific Reports, 2017, 7(1): 5636.

[14] Wei T, Jie Y, Yan H, Li F, Zhou Q, Wei C, Bennetzen J L. Circular RNAarchitecture and differentiation during leaf bud to young leaf developmentin tea (Camellia sinensis). Planta, 2018, 248(10): 1-13.

[15] Darbani B, Noeparvar S, Borg S. Identification of circular RNAs from the parental genes involved in multiple aspects of cellular metabolism in barley. Frontiers in Plant Science, 2016, 7: 776.

[16] Chen L, Zhang P, Fan Y, Lu Q, Li Q, Yan J, Muehlbauer G J, Schnable P S, DaiM, Li L. Circular RNAs mediated by transposons are associated withtranscriptomic and phenotypic variation in maize. New Phytologist, 2018, 217(3): 3.

[17] Tan J, Zhou Z, Niu Y, Sun X, Deng Z. Identification and functional characterization of tomato circrnas derived from genes involved in fruit pigment accumulation. Scientific Reports, 2017, 7: 8594.

[18] Cheng J, Zhang Y, Li Z, Wang T, Zhang X, Zheng B.A lariat-derived circular RNA is required for plant development in. Science China Life Sciences, 2018, 61(2): 204-213.

[19] Pan T, Sun X, Liu Y, Li H, Deng G, Lin H, Wang S.Heat stress alters genome wide profiles of circular RNAs in. Plant Molecular Biology, 2018, 96(3): 217-229.

[20] LI N, LIU T T, GUO F, YANG J W, SHI Y G, WANG S G, SUN D Z. Identification of long non-coding RNA-microRNA-mRNA regulatory modules and their potential roles in drought stress response in wheat (L.). Frontiers in Plant Science, 2022, 10: 1011064.

[21] Quan X, Zeng J, Ye L, Chen G, Han Z, Shah J, Zhang G.Transcriptome profiling analysis for two Tibetan wild barley genotypes in responses to low nitrogen. BMC Plant Biology, 2016, 16(1): 30-45.

[22] Sun Y, Song K, Sun L, Qin Q, Jiang T, Jiang Q, Xue Y.Morpho-Physiological and transcriptome analysis provide insights into the effects of zinc application on nitrogen accumulation and metabolism in wheat (L.). Plant Physiology and Biochemistry, 2020, 149: 111-120.

[23] Kim D, Langmead B, Salzberg S L. HISAT: A fast spliced aligner with low memory requirements. Nature methods, 2015, 12(4): 357-360.

[24] Memczak S, Jens M, Elefsinioti A,Torti F, Krueger J, Rybak A, Maier L, Mackowiak S D, Gregersen L H, Munschauer M, Loewer A, Ziebold U, Landthaler M, Kocks C, Noble F, Rajewsky N. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature, 2013, 495(7441): 333-338.

[25] Zhu Y X, Jia J H, Yang L, Xia Y C, Zhang H L, Jia J B, Zhou R, Nie P Y, Yin J L, Ma D F, Liu L C. Identification of cucumber circular RNAs responsive to salt stress. BMC Plant Biology, 2019, 19(1): 164.

[26] Love M I, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome biology, 2014, 15(12): 550.

[27] Bo X, Wang S. Target Finder: a software for antisense oligonucleotide target site selection based on MAST and secondary structures of target mRNA. Bioinformatics, 2005, 21(8):1401.

[28] Ma S W, Wang M, Wu J H, Guo W L, Chen Y M, Li G W, Wang Y P, Shi W M, Xia G M, Fu D L, Kang Z S, Ni F. WheatOmics: A platform combining multiple omics data to accelerate functional genomics studies in wheat. Molecular Plant, 2021, 14(12): 1965-1968.

[29] Conesa A, G?tz S.Blast2GO: a comprehensive suite for functional analysis in plant genomics. International Journal of Plant Genomics,2008, 2008: 619832.

[30] YANG C, LU X, MA B, CHEN S Y, ZHANG J S. Ethylene signaling in rice and: conserved and diverged aspects. Molecular Plant, 2015, 8(4): 495-505.

[31] Hu B, Jiang Z, Wang W, Qiu Y, Zhang Z, Liu Y, Li A, Gao X, Liu L, Qian Y, Huang X, Yu F, Kang S, Wang Y, Xie J, Cao S, Zhang L, Wang Y, Xie Q, Kopriva S, Chu C. Nitrate-NRT1.1B-SPX4 cascade integrates nitrogen and phosphorus signalling networks in plants. Nature Plants, 2019, 5(4): 401-413.

[32] Park J J, Yi J, Yoon J, Cho L H, Ping J, Jeong H J, Cho S K, Kim W T, An G. OsPUB15, an E3 ubiquitin ligase, functions to reduce cellular oxidative stress during seedling establishment. The Plant Journal, 2011, 65(2): 194-205.

[33] Li X M, Chao D Y, Wu Y, Huang X H, Chen K, Cui L G, Su L, Ye W W, Chen H, Chen H C, Dong N Q, Guo T, Shi M, Feng Q, Zhang P, Han B, Shan J X, Gao J P, Lin H X. Natural alleles of a proteasome α2 subunit gene contribute to thermotolerance and adaptation of African rice. Nature Genetics, 2015, 47(7): 827-833.

[34] Ai P H, Sun S B, Zhao J N, Fan X R, Xin W J, Guo Q, Yu L, Shen Q R, Wu P, Miller A J, Xu G H. Two rice phosphate transporters, OsPht1;2 and OsPht1;6, have different functions and kinetic properties in uptake and translocation. The Plant Journal, 2009, 57(5): 798-809.

[35] Louren?o T, Sapeta H, Figueiredo D D, Rodrigues M, Cordeiro A, Abreu I A, Saibo N J, Oliveira M M. Isolation and characterization of rice (L.) E3-ubiquitin ligase OsHOS1 gene in the modulation of cold stress response. Plant Molecular Biology, 2013, 83(4/5): 351-363.

[36] Sun S K, Xu X, Tang Z, Tang Z, Huang X Y, Wirtz M, Hell R, Zhao F J. A molecular switch in sulfur metabolism to reduce arsenic and enrich selenium in rice grain. Nature Communications, 2021, 12: 1392.

[37] Liao Y D, Lin K H, Chen C C, Chiang C M.protein phosphatase 1a (OsPP1a) involved in salt stress tolerance in transgenic rice. Molecular Breeding, 2016, 36: 22.

[38] Giri J, Vij S, Dansana P K, Tyagi A K. Rice A20/AN1 zinc- finger containing stress-associated proteins (SAP1/11) and a receptor- like cytoplasmic kinase (OsRLCK253) interact via A20 zinc-finger and confer abiotic stress tolerance in transgenicplants.New Phytologist, 2011, 191(3): 721-732.

[39] Mukhopadhyay A, Vij S, Tyagi A K. Overexpression of a zinc- finger protein gene from rice confers tolerance to cold, dehydration, and salt stress in transgenic tobacco. Proceedings of the National Academy of Sciences of the USA, 2004, 101(16): 6309-6314.

[40] Errichelli L, Dini M S, Laneve P, Colantoni A, Legnini I, Capauto D,Rosa A, De Santis R, Scarfo R, Peruzzi G. FUS affects circular RNAexpression in murine embryonic stem cell-derived motor neurons. Nature Communications, 2017, 8: 14741.

[41] Yin J L, Ma D F, Liu L C, Xia Y C, Zhu Y X.Biology features of circular RNAs and their research progress in plants. Acta Botanica Boreali-Occidentalia Sinica, 2017, 37(12): 2510-2518.

[42] Lu T, Cui L, Zhou Y, Zhu C, Fan D, Gong H, Zhao Q, Zhou C, Zhao Y, Lu D. Transcriptome-wide investigation of circular RNAs in rice. RNA, 2015, 21(12): 2076-2087

[43] Hansen T B, Jensen T I, Clausen B H, Bramsen J B, Finsen B, Damgaard C K. Natural RNA circles function as efficient microRNA sponges. Nature,2013, 495: 384-388.

[44] Meloni D, Oliva M, Martinez C, Cambraia J. Photosynthesis and activity of superoxide dismutase, peroxidase and glutathione reductase in cotton under salt stress. Environmental and Experimental Botany, 2003, 49: 69-76.

[45] Gutterson N, Reuber T L.Regulation of disease resistance pathways by AP2/ERF transcription factors.Current Opinion in Plant Biology, 2004, 7(4): 465-471.

[46] NELSON D R. Plant cytochrome P450s from moss to poplar.Phytochemistry Reviews, 2006, 5(2/3): 193-204.

[47] SCHULZ P, HERDE M, ROMEIS T. Calcium-dependent protein kinases: hubs in plant stress signaling and development. Plant physiology, 2013, 163(2): 523-530.

[48] DAS R, PANDEY G. Expressional analysis and role of calcium regulated kinases in abiotic stress signaling. Current genomics, 2010, 11(1): 2-13.

[49] Li M Z, Li M F, Li D D, Wang S M, Yin H J. Overexpression of theaquaporin,, promotes plant growth and stress tolerance. International Journal of Molecular Sciences, 2021, 22(4): 2112.

[50] CaoY F, Wu Y F, ZhengZ, SongF G. Overexpression of the rice EREBP-like geneenhances disease resistance and salt tolerance in transgenic tobacco. Physiological and Molecular Plant Pathology, 2006, 67(3/5): 202-211.

[51] LiuC W, FukumotoT, MatsumotoT, GenaP, FrascariaD, KanekoT, KatsuharaM, ZhongS H, SunX L, ZhuY M, IwasakiI, DingX D, CalamitaG, KitagawaY.Aquaporinpromotes rice salt resistance and seed germination. Plant Physiology and Biochemistry, 2013, 63: 151-158.

Identification of wheat circular RNAs responsive to drought stress

LI Ning, LIU Kun, LIU TongTong, SHI YuGang, WANG ShuGuang, YANG JinWen, SUN DaiZhen

College of Agriculture, Shanxi Agricultural University, Taigu 030801, Shanxi

【Objective】Drought is one of the foremost abiotic stress limiting global wheat production. Exploring the molecular mechanism of wheat response to drought stress have great significance in wheat molecular breeding. Circular RNAs (circRNAs) have been proved to play an important role in the process of plants tolerance to environmental stresses. Therefore, identifying circRNAs involved in drought stress response will help to construct a regulatory network of wheat drought tolerance, and lay a foundation for analyzing the drought resistance mechanism in wheat. 【Method】In this study, two wheat varieties (Zhoumai13 and Jimai38) with significant differences in drought resistance were used and circRNA-seq was performed on their root samples under well-watered and drought conditions. Differentially expressed circRNAs related to drought stress response were screened based on the identified circRNAs and their microRNAs (miRNAs) targets were predicted. Further, potential circRNA-miRNA-mRNA regulatory modules related to wheat drought stress response were constructed according to the expression patterns of miRNAs and their target genes under drought stress. 【Result】A total of 1 409 wheat circRNAs were identified, most of which (68.91%) were exonic circRNAs. Only 133 circRNAs were simultaneously identified in both varieties. A total of 239 differentially expressed circRNAs were identified under drought stress, of which 138 circRNAs were specifically differentially expressed in the drought-resistant variety Zhoumai 13 (ZM13), and 19 circRNAs were differentially expressed simultaneously in both varieties. Besides, 34 targeted miRNAs and 1 408 miRNA target genes were predicted. Based on the expression patterns of these differentially expressed circRNAs, targeted miRNAs, and miRNA target genes, five potential circRNA-miRNA-mRNA regulatory modules centered on tae-miR9664-3p, tae-miR1122b-3p, tae-miR9662a-3p, tae-miR6197-5p and tae-miR1120c-5p in response to drought stress were screened. 【Conclusion】Wheat circRNAs have obvious specificity in different cultivars and different expression patterns among different drought-tolerant wheat cultivars. A total of 239 wheat circRNAs and 5 potential circRNA-miRNA-mRNA regulatory modules in response to drought stress were identified in the present study.

wheat; circular RNAs; drought stress; microRNAs

10.3864/j.issn.0578-1752.2022.23.002

2022-07-25；

2022-09-05

山西省基礎研究計劃（20210302124148）、山西省高等學校科技創(chuàng)新項目（2021L124）、山西農(nóng)業(yè)大學科技創(chuàng)新基金（2020BQ30）

李寧，E-mail：13159862006@163.cm。通信作者楊進文，E-mail：yang_jin_wen@126.com。通信作者孫黛珍，E-mail：sdz64@126.com

（責任編輯 李莉）