紫花苜蓿細胞質雄性不育系相關microRNA的鑒定

郭 強，王英哲，陳晶晶，周 仂，徐 博*

(1. 吉林農業大學動物科學技術學院，吉林 長春 130118； 2. 吉林省農業科學院，吉林 長春 130033)

紫花苜蓿(MedicagosativaL.)是一種多年生豆科苜蓿屬植物。1958年，加拿大學者發現紫花苜蓿的雄性不育現象，細胞質雄性不育(Cytoplasmic male sterility，CMS)指開花植物不能產生功能性花藥、花粉或雄配子的一種遺傳學現象[1]，這種現象主要由細胞質因子和細胞器基因控制[2]，其不育特性表現為雄器官的退化、畸形或功能缺失，導致花藥、花粉母細胞和雄配子的生長發育產生異常[3-4]。常用于作物的雜交育種。

隨著雄性不育系的發現，國外的眾多學者開始采用分子技術、組織培養法研究遺傳機理，并利用不育系培育雜交組合和轉基因材料[5-7]；1978年，吳永敷等[8]在我國內蒙地區從‘草原1號’中發現了6株紫花苜蓿雄性不育系；2008年，于洪柱等[9]發現了紫花苜蓿細胞質雄性不育材料—MS-GN。隨著紫花苜蓿細胞質雄性不育系及保持系的發現，對其機理的探索也在不斷深入。同時，國內學者[10-14]對不育系展開了同工酶分析、細胞學觀察、生殖生物學、生理生化機制等內容的研究；利用紫花苜蓿細胞質雄性不育系選育雜交組合并對其F1代進行生物學性狀及抗逆性研究[15-18]；有學者[19-21]利用SSR標記以及BSR技術對紫花苜蓿細胞質雄性不育恢復基因進行了定位分析，并且對紫花苜蓿細胞質雄性不育系及其保持系的花藥發育過程進行切片觀察，測定其在不同發育階段生理生化指標的變化。

植物microRNA是一類長度為20～26 nt的內源性非編碼小RNA[22]。成熟的microRNA長度為20～24 nt，多為21 nt[23]，其中長度為21～24 nt的microRNA占主體地位。植物microRNA參與植物體內許多重要的生命過程，擬南芥(ArabidopsisthalianaL.)體內的miR160可通過負調控生長素響應因子進一步調控根冠細胞的發育[24]；Zea-miR390則可以調控玉米(ZeamayL.)莖尖細胞的發育及葉片的形態建成[25]；miR172，miR159，miR156/157則可以調控開花時期[26-27]。microRNA也參與調控植物的環境適應性，如miR398與植物的耐逆性密切相關[28]。另外，microRNA參與植物激素合成代謝與信號傳導，如miR160，miR164，miR167，miR528等microRNA均與生長素的信號傳遞有關，其中miR167還可以參與脫落酸等其他的信號通路的傳導[29-30]。

基于全基因組信息可知，在許多雄性不育植物的花粉發育過程中均有大量的microRNA參與[31-33]。在擬南芥中，miR165/166參與小孢子的發育，miR156和miR159負責調控植物絨氈層細胞的發育過程，miR159的靶基因屬于GAMYB(Gibberellins acid v-myb avian myeloblastosis viral oncogene homolog)家族，主要調控植物花藥細胞程序性死亡(Programmed Cell Death,PCD)過程中的赤霉素的傳導，miR159的異常表達會使絨氈層無法降解，最終使小孢子不能正常發育而喪失育性[34-36]。擬南芥中若缺乏miR165/166會使花藥只產生兩個花粉囊，最終會導致擬南芥的育性降低[37]；水稻(OryzasativaL.)中的OsmiR528，OsmiR159，OsmiR172d；玉米中Zea-miR397；油菜(BrassicanapusL.)中的bra-miR167，bra-miR172等miRNA都可能與植物的小孢子發育有關[31,38]。另外，有研究發現，某些miRNA可通過參與植物激素合成和應答過程來調控花粉的發育，例如miR167缺失會引起茉莉酸的合成異常，從而導致花藥不開裂以及散粉受阻；miR167的過表達同樣會引起植物不育現象的產生[30,39-41]。

植物microRNA是揭示植物雄性不育機制的有效方式之一，基于之前的研究發現：紫花苜蓿細胞質雄性不育的敗育可能與線粒體基因表達異常、淀粉蔗糖含量虧缺、絨氈層發育異常等現象有關[17,42]，因此，為深入探索不育機制，本試驗利用small RNA-seq技術，找到雄性不育相關的miRNA，為更好地利用紫花苜蓿不育系進行雜交育種工作提供理論基礎，以期選育出更具優良性狀的種質資源。

1 材料與方法

1.1 試驗材料

本試驗所選用的材料為紫花苜蓿細胞質雄性不育系MSGN1A及其對應的同型保持系MSGN1B，MSGN1A和MSGN1B是以具有雄性不育性的‘公農3號’為供體母本株系，通過保持系篩選，回交4代后選育出來的。所有植物材料均由吉林省農業科學院提供。MSGN1A和MSGN1B的個體植株于2016年4月在吉林省長春市吉林農業大學草業科學試驗田種植。選取現蕾期24 h的花藥部分，用于提取總RNA，備用。試驗設置3個生物學重復。

1.2 試驗方法

1.2.1總RNA的提取及質量控制 參照Garg等[43]的方法提取MSGN1A和MSGN1B花藥總RNA，用紫外分光光度計分別在230，260和280 nm處測定吸光值(OD230，OD260和OD280)，OD260/OD280和OD260/OD230均須≥1.8。采用Agilent 2100生物分析儀檢測總RNA的完整性，總RNA完整性值在8.0～10.0的樣品可用于后續試驗。

1.2.2small RNA文庫構建 將提取到的總RNA送北京百邁克生物技術公司，完成small RNA文庫的構建及測序工作。提取兩個樣品總RNA，使用聚丙烯酰胺凝膠電泳法根據片段的長度分離RNA并回收長度在18～30 nt之間的條帶，構建2個small RNA文庫。

1.2.3miRNA的鑒定與分析 對已知miRNA的鑒定，采用miRBase(v21)軟件將比對上參考基因組的reads序列與known miRNA數據庫中的成熟miRNA序列進行序列比對。若序列與known miRNA完全一致，則這個reads將被認為是鑒定到的known miRNA。

考慮到miRNA的生物學特征，對于沒有被鑒定到known miRNA的序列，使用miRDeep 2軟件對novel miRNA做預測。利用miRDeep 2軟件包，將reads比對到基因組上的位置信息，得到可能的前體序列，根據reads在前體序列上的分布信息及前體結構能量信息，使用貝葉斯模型經打分最終實現novel miRNA的預測。

1.2.4差異表達miRNA篩選及其靶基因功能注釋 以|log2(FC)|≥1；FDR≤0.01作為標準篩選差異表達的microRNA。差異倍數(Fold Change，FC)表示兩樣品間表達量的比值。

依據known miRNA和novel miRNA與對應物種的基因序列信息，使用Target Finder軟件[40]對miRNA靶基因進行預測。

基因功能注釋由北京百邁客生物技術公司進行，差異表達miRNA靶基因的功能注釋所使用的數據庫為GO (Gene Ontology，http://www.geneontology.org/)和KEGG (Kyoto Encyclopedia of Genes and Genomes，http://www.genome.jp/kegg/)。

2 結果與分析

2.1 測序數據分析

對測序得到的數據進行質量控制，每個樣品的Clean Data均大于14.91 M。細胞質雄性不育系MSGN1A及其同型保持系MSGN1B的Raw reads分別為27.83 M和28.28 M，Clean reads分別為20.24 M和14.91 M，Q30(質量值大于30)分別為98.87%和99.03%。

表1 測序數據產出統計表

2.2 miRNA的鑒定

為了進一步鑒定已知的miRNA，將比對上參考基因組的reads序列與已知miRNA數據庫miRBase(v21)中的成熟miRNA序列進行比對。序列與已知miRNA完全相同的reads被認為是鑒定到的已知miRNA。針對miRNA的生物特征，對于未鑒定到known miRNA的序列，利用miRDeep 2軟件進行novel miRNA的預測。通過鑒定，本試驗共檢測到1 478個miRNA，其中已知miRNA 1 351個，novel miRNA 127個(表2)。

表2 miRNA統計結果

試驗檢測到1 351個known miRNA，分屬于136個miRNA家族，部分結果如圖1所示。每個家族中所含miRNA的數量和表達水平不盡相同，其中包括屬于miR166家族的247個miRNA，miR167家族的159個miRNA，miR160的家族124個的miRNA，miR172家族的105個miRNA，mi156家族的98個miRNA，miR319家族的88個miRNA等。成熟miRNA長度主要集中在20 nt到24 nt的范圍內，其中植物的miRNA以21 nt或24 nt為主，試驗鑒定出的novel miRNA的長度主要為21 nt，其次是24 nt和22 nt。

圖1 部分已知miRNA在miRNA家族中所含數量

2.3 差異表達miRNA的篩選與鑒定

在差異表達miRNA檢測過程中，使用|log2(FC)|≥1；FDR≤0.01作為篩選標準。試驗共計篩選出差異表達miRNA 130個，包括90個上調表達的miRNA和40個下調表達的miRNA。試驗共鑒定到40個差異表達的已知miRNA，分屬于15個miRNA家族，與保持系MSGN1B相比，紫花苜蓿細胞質雄性不育系MSGN1A中有24個上調表達miRNA，有16個下調表達的差異miRNA。

表3 差異表達的已知miRNA

2.4 miRNA靶基因預測

根據known miRNA和novel miRNA與其對應物種的基因序列信息，使用Target Finder軟件進行靶基因的預測，1 351個known miRNA中，有1 347個miRNA預測到靶基因，靶基因數目為9 848個；127個novel miRNA中，有126個miRNA預測到8 252個靶基因。注釋到的靶基因數目為16 458個。

表4 miRNA靶基因數目預測統計

2.5 差異表達miRNA靶基因注釋

試驗共得到16 458個靶基因，其中10 637個靶基因獲得功能注釋信息，紫花苜蓿細胞質雄性不育系MSGN1A與其同型保持系MSGN1B間注釋的差異表達miRNA靶基因在各注釋庫中的數目統計見,表5，共有10 637個miRNA靶基因被各個注釋庫注釋到，其中COG注釋庫中有3 104個miRNA靶基因被注釋到，包括815個差異表達miRNA靶基因；GO注釋庫中有5 147個miRNA靶基因被注釋到，包括1 264個差異表達miRNA靶基因；KEGG注釋庫中有3 479個miRNA靶基因被注釋到，包括861個差異表達miRNA靶基因；KOG注釋庫中有5 620個miRNA靶基因被注釋到，包括1 286個差異表達miRNA靶基因；Pfam注釋庫中有8 309個miRNA靶基因被注釋到，包括1 871個差異表達miRNA靶基因；Swissprot注釋庫中有6 535個miRNA靶基因被注釋到，包括1 662個差異表達miRNA靶基因；eggNOG注釋庫中有9 133個miRNA靶基因被注釋到，包括2 261個差異表達miRNA靶基因；Nr注釋庫中有6 535個miRNA靶基因被注釋到，包括2 612個差異表達miRNA靶基因。

表5 miRNA靶基因注釋統計表

2.5.1GO注釋分類 對差異表達miRNA靶基因進行GO分類，共計1 264個miRNA靶基因得到注釋。按照GO的三個大類進行下一次級分類，共分成41個功能群，其中包括18個生物學過程，12個細胞成分，11個分子功能(圖2)。

圖2 差異表達miRNA靶基因GO注釋分類

在生物學過程中，“代謝過程”所占比例最大，其次是“細胞過程”、“信號生物過程”、“生物調控”、“刺激應答”、“定位”、“細胞成分組織或生物合成”等生物學過程，它們的占比均大于10%；“發育過程”、“多細胞生物過程”、“繁殖過程”、“信號”、“多生物過程”、“生長”、“繁殖”、“免疫系統過程”等生物學過程的占比均在1%～10%之間；“生物附著”、“生物相”、“節律過程”的占比在0.1%～1.0%之間。

在細胞組分類別中，“細胞部分”和“細胞”占比最大，緊接著是“細胞器”、“膜”、“細胞器部分”、“高分子配合物”，其占比均大于10%；“胞外區”和“細胞連接”的占比在1%～10%之間；“膜封閉腔”、“病毒”以及“病毒部分”的占比在0.1%～1.0%之間。

在分子功能類別中，“合成”和“催化活性”是最大的功能類別；其次是“轉運蛋白活性”、“核酸合成轉錄因子活性”、“電子傳遞體活性”、“結構分子活性”、“酶調節活性”以及“分子傳感器活性”，占比均在占比在1%～10%之間；最后是“抗氧化劑活性”、“受體活性”、“蛋白合成轉錄因子活性”和“鳥苷酸交換因子活性”，占比均不到1%。

圖3 差異表達miRNA靶基因KEGG通路注釋

2.5.2KEGG通路分析 對差異表達miRNA靶基因進行KEGG通路分析，差異表達miRNA靶基因所涉及的通路的類別主要包括：“細胞進程”、“環境信息過程”、“基因信息過程”、“代謝”、“生物系統”等五個大類，在此基礎上進一步分為82個類別，其中有24個靶基因注釋到“代謝”中的“氨基酸生物合成”通路中，占比為5.61%；18個靶基因注釋為“內吞作用”通路中，占比4.21%；16個靶基因注釋到“核糖體”通路中，占比為3.74%；“RNA轉運”和“碳代謝”通路中分別有15個靶基因被注釋到，占比均為3.50%；“內質網中的蛋白質過程”、“泛素介導的蛋白水解作用”、“植物激素信號傳導”通路中分別有14個靶基因被注釋到，占比均為3.27%；有13個靶基因被注釋到“剪接體”通路中，占比為3.04%；“氨基糖和核糖代謝”和“氧化磷酸化”通路中有12個靶基因被注釋到，占比均為2.8%；“淀粉蔗糖代謝”通路中有8個靶基因被注釋到，占比為1.87%。

2.6與花藥發育相關的miRNA

在差異表達的miRNA中，mtr-miR5256的靶基因的功能注釋為“glucan ando-1,3-beta-glucosidase 4”,該基因參與淀粉蔗糖代謝。miR159靶基因的功能注釋為“MYB transcription factor”，“MYB轉錄因子”是花分生組織識別基因的啟動子。miR172靶基因功能注釋為“AP2 Domain transcription factor”，AP2類似轉錄因子主要調控花器分化和花期。miR396靶基因功能注釋為“生長調控因子HSP70蛋白”，miR169靶基因功能注釋為“CBF HAP2 like protein”，其主要作用是調控花期，這些差異表達的miRNA均有可能導致花粉敗育。

3 討論

MiR159靶基因的功能注釋為“MYB transcription factor”，“MYB轉錄因子”是花分生組織識別基因的啟動子。MYB是谷物糊粉蛋白細胞中赤霉素(Gibberellic acid，GA)信號的組成部分，在花的發育過程中起著重要作用[44]。陶園[45]認為，擬南芥轉錄因子ICE1(Inducer of CBF Expression 1)通過調控GAMYB基因來調節赤霉素信號傳導，從而影響花粉育性并調控植物生長。miR159的過量表達會引起與花發育過程相關的轉錄因子MYB的表達異常，從而進一步影響到花粉絨氈層的正常代謝，包括絨氈層的形態和結構異常等[46]。在本研究中miR159b和miRNA159c均下調表達，miRNA159d上調表達，這些miRNA的表達異常讓我們有理由認為miR159家族與紫花苜蓿細胞質雄性不育有著密切聯系。

試驗所得結果中miR172，miR172a，miR172d相較于保持系來說，在不育系中均上調表達，miR172靶基因功能注釋為“AP2 Domain transcription factor”，AP2-like轉錄因子主要調控花器分化和花期。AP2-like轉錄因子對花器官的生成是至關重要的，如果花器官發育出現異常情況(如花絲和花粉囊發育異常)就很可能會使花粉粒無法正常發育，從而導致花粉敗育[47]。Chen等[48]研究認為，miR172過表達會導致植株的花期提前以及花器官發育畸形，其表型與ap2突變體表型類似，最終導致配子異常發育。因此，miR172在不育系中的表達量降低也可能導致花粉敗育，具體的調控模式還需進一步研究。

miR396靶基因功能注釋為“生長調控因子HSP70蛋白”。試驗結果顯示miR396家族的miR396，miR396a，miR396b-5p，miR396c-5p，miR396d，miR396e等miRNA在不育系中上調表達，而miR396在初級和成熟的miRNA中都具有高度保守性，是“擬南芥生長調節因子”家族成員[49-50]。Yang等[51]認為miRNA396介導的NtGRF-like基因表達調控通路與煙草葉片和花的發育過程密不可分，這些miRNA對應的靶基因調控花的發育過程，也就是說，miR396的異常表達可能影響花的發育。

由結果可知，與保持系相比，miR169c在不育系中上調表達，miR169靶基因功能注釋為“CBF HAP2 like protein”，其主要作用是調控花期，擬南芥miR169家族成員是發育過程和刺激通路中普遍存在的調控因子[52]。因此，miR169的表達量降低會可能導致花期的提前或延遲，影響花藥的發育等過程。

4 結論

差異表達microRNA靶基因功能注釋結果顯示，與紫花苜蓿細胞質雄性不育相關的差異表達microRNA家族主要有4個，包括miR159，miR172，miR396，miR169；其靶基因的功能主要是調控花藥的發育，其次是調控花期以及信號傳導；micro RNA靶基因主要涉及的通路有：“代謝過程”、“信號生物過程”、“植物激素信號傳導”、“氧化磷酸化”和“淀粉蔗糖代謝”等，而差異表達microRNA靶向的轉錄因子中多為控制花期、花藥發育和信號傳導等過程的因子，對不育系的敗育起到一定的影響。本試驗結果可初步揭示紫花苜蓿不育系的分子機理，為利用紫花苜蓿不育系進行雜交育種提供理論基礎，從而選育出具有優良性狀的紫花苜蓿種質資源。