參與調控意大利蜜蜂工蜂中腸基因表達的東方蜜蜂微孢子蟲miRNA的組學解析及其調控網絡

范小雪, 杜 宇, 張文德, 王 杰, 蔣海賓, 范元嬋, 馮睿蓉,萬潔琦, 周紫彧, 熊翠玲,2, 鄭燕珍,2, 陳大福,2,3, 郭 睿,2,3,*

(1. 福建農林大學動物科學學院(蜂學學院), 福州 350002; 2. 福建農林大學蜂療研究所, 福州 350002;3. 福建農林大學蜂產品加工與應用教育部工程研究中心, 福州 350002)

東方蜜蜂微孢子蟲Nosemaceranae是一種專性侵染成年蜜蜂中腸上皮細胞的單細胞真菌病原，可對蜜蜂宿主造成腸道破壞、能量脅迫、免疫抑制、消化紊亂及壽命縮短等不良影響(Higesetal., 2007; Mayack and Naug, 2009; Alauxetal., 2010; Botíasetal., 2013; Chenetal., 2013)，但到目前為止二者互作的分子機制尚未闡明。

微小RNA(microRNA, miRNA)是一類長度約為18～25 nt且高度保守的內源性單鏈非編碼RNA(non-coding RNA, ncRNA)(Bartel, 2004)，可通過靶向結合mRNA使其降解或抑制其翻譯，從而在轉錄后水平調控基因表達，進而發揮廣泛而重要的生物學功能(Pillai, 2005)。此外，miRNA已被證實能夠作為媒介介導真菌病原及其宿主的互作(Linetal., 2016; Crostonetal., 2018)。2012年，張辰宇團隊研究發現植物來源的miR-168a可在動物體內穩定積累表達，并通過靶向抑制小鼠低密度脂蛋白LDLRAP1的表達水平，降低小鼠血漿LDL的去除率，首次證實miRNA的跨界調控作用(Zhangetal., 2012)。較多的研究證實結果表明小RNA(small RNA, sRNA)可在真菌病原及其宿主之間雙向傳播，真菌病原可將sRNA作為效應因子轉移到宿主細胞內，參與宿主的基因表達調控，從而影響自身毒力、宿主免疫力以及疾病進程等。例如，灰霉菌Botrytiscinerea和球孢白僵菌Beauveriabassiana來源的sRNA在細胞囊泡的包裹和保護作用下轉移到宿主體內，通過與宿主Argonaute 1蛋白結合并劫持宿主RNAi機制，進而選擇性沉默宿主免疫基因，以達到有效感染的目的(Weibergetal., 2013; Cuietal., 2019)。對于miRNA介導的西方蜜蜂Apismellifera及東方蜜蜂微孢子蟲互作，前人開展了少量研究(Evans and Huang, 2018; Huangetal., 2019)，但總體進展緩慢。前期研究中，筆者所在團隊聯用鏈特異性建庫的RNA-seq和small RNA-seq(sRNA-seq)技術對東方蜜蜂微孢子蟲感染7 d和10 d的意大利蜜蜂Apismelliferaligustica工蜂中腸及未感染中腸進行測序，基于高質量的全轉錄組數據系統解析了lncRNA, miRNA及ceRNA調控網絡介導的宿主免疫應答(Chenetal., 2019; 付中民等, 2019)；此外，還在lncRNA, miRNA和mRNA組學水平深入分析和探討了東方蜜蜂微孢子蟲侵染意大利蜜蜂工蜂的分子機制(耿四海等, 2020a, 2020b)。近期，筆者所在團隊基于已獲得的東方蜜蜂微孢子蟲感染及未感染的意大利蜜蜂工蜂中腸轉錄組數據(付中民等, 2019)以及純凈的病原孢子全轉錄組數據(耿四海等, 2020a)篩選出宿主的差異表達miRNA(differentially expressed miRNA, DEmiRNA)和病原的差異表達mRNA(differentially expressed mRNA, DEmRNA)，通過生物信息學手段預測宿主DEmiRNA和病原DEmRNA的靶向關系，并對二者間的調控網絡進行了構建和分析，進一步結合本團隊的前期研究結果和前人的相關研究結果探討了意大利蜜蜂工蜂的DEmiRNA對東方蜜蜂微孢子蟲基因和通路的調控作用(未發表數據)。在長期的協同進化中，東方蜜蜂微孢子蟲的基因組大大簡化，多數物質和能量代謝的途徑丟失，其增殖高度依賴宿主細胞提供物質和能量(Burrietal., 2006)。同時東方蜜蜂微孢子蟲也進化出一些特殊的能力，能夠對宿主的生理生化功能進行操縱，如東方蜜蜂微孢子蟲侵染可對西方蜜蜂工蜂造成較強的能量脅迫，致使宿主的能量代謝速率加快且糖水的攝入量提高(Antúnezetal., 2009; Chenetal., 2013)。然而，對于東方蜜蜂微孢子蟲通過何種方式和機制對西方蜜蜂的基因表達和生理活動進行調控，相關研究仍然匱乏，迄今僅有一例報道(Evans and Huang, 2018)。此前，Evans和Huang(2018)通過高通量測序和生物信息學分析方法在東方蜜蜂微孢子蟲感染的西方蜜蜂工蜂中腸中預測出病原的6條miRNA，進一步分析發現其中5條miRNA不僅能靶向病原本身的1 545條mRNA，還可能被分泌到宿主細胞的胞質內，通過靶向918條宿主mRNA調控宿主的新陳代謝、細胞凋亡和免疫防御等過程。

本研究擬通過比較分析篩選出侵染意大利蜜蜂工蜂的東方蜜蜂微孢子蟲的DEmiRNA以及意大利蜜蜂工蜂響應東方蜜蜂微孢子蟲侵染的意大利蜜蜂DEmRNA，通過生物信息學方法預測病原DEmiRNA靶向的宿主mRNA和DEmRNA，并通過注釋數據庫獲得相應的功能和通路信息，進一步結合東方蜜蜂微孢子蟲侵染西方蜜蜂的生物學背景，筆者所在課題組的前期研究結果以及前人的相關研究結果深入探討東方蜜蜂微孢子蟲的DEmiRNA對意大利蜜蜂工蜂中腸基因表達的調控作用。研究結果可為探明東方蜜蜂微孢子蟲的侵染機制、東方蜜蜂微孢子蟲-意大利蜜蜂工蜂互作機制提供重要的參考信息和基礎。

1 材料與方法

1.1 東方蜜蜂微孢子蟲的miRNA組學數據來源

前期研究中，筆者所在課題組分別提取了3個東方蜜蜂微孢子蟲的純凈孢子樣品NcCK(NcCK1, NcCK2, NcCK3)的總RNA，并按照已建立的方法(耿四海等, 2020a)構建了相應的cDNA文庫，單端測序由廣州基迪奧生物科技有限公司完成，測序平臺為Illumina MiSeq。通過連續比對西方蜜蜂參考基因組(Assembly Amel_4.5)、GenBank數據庫、Rfam數據庫及東方蜜蜂微孢子蟲參考基因組(Assembly ASM98816v1)篩濾得到比對上東方蜜蜂微孢子蟲參考基因組的東方蜜蜂微孢子蟲純凈孢子miRNA組學數據(NcCK)，侵染意大利蜜蜂工蜂7 d(NcT1)和10 d(NcT2)的東方蜜蜂微孢子蟲miRNA組學數據(耿四海等, 2020a)。測序數據已上傳至NCBI SRA數據庫，NcCK的Bioproject登記號為PRJNA395264。通過對下機的原始數據進行質控得到了最終的有效標簽序列(clean tags)(耿四海等, 2020a)。

1.2 東方蜜蜂微孢子蟲的DEmiRNA篩選

按照耿四海等(2020a)的方法進行東方蜜蜂微孢子蟲miRNA的鑒定和表達量計算以及NcCKvsNcT1和NcCKvsNcT2比較組顯著性DEmiRNA的篩選，篩選出的顯著性DEmiRNA(|log2fold change|≥1且P≤0.05)用于本研究中靶向意大利蜜蜂工蜂中腸mRNA和DEmRNA的預測和分析。

1.3 意大利蜜蜂工蜂中腸的mRNA組學數據來源

筆者所在團隊前期已對接種50%(w/v)蔗糖溶液7 d(AmCK1)和10 d(AmCK2)的意大利蜜蜂工蜂中腸樣品以及接種含東方蜜蜂微孢子蟲孢子的50%(w/v)蔗糖溶液7 d (AmT1)和10 d(AmT2)的意大利蜜蜂工蜂中腸樣品進行RNA抽提、cDNA文庫構建、轉錄組測序以及測序數據質控(付中民等, 2019)。測序數據已上傳至NCBI SRA數據庫，Bioproject號為PRJNA395264。質控后的意大利蜜蜂工蜂中腸mRNA組學數據可用于本研究中東方蜜蜂微孢子蟲的顯著性DEmiRNA的靶標預測及分析。

1.4 東方蜜蜂微孢子蟲顯著性DEmiRNA靶向意大利蜜蜂工蜂中腸mRNA的預測

為探究1.2節篩選出的東方蜜蜂微孢子蟲DEmiRNA調控意大利蜜蜂工蜂中腸基因表達的廣泛性，利用TargetFinder軟件(Allenetal., 2005)預測NcCKvsNcT1和NcCKvsNcT2比較組中顯著性DEmiRNA靶向的1.3節的意大利蜜蜂工蜂中腸的全部mRNA，采用默認參數。

1.5 意大利蜜蜂工蜂中腸響應東方蜜蜂微孢子蟲侵染的DEmRNA的篩選及分析

miRNA一般對基因進行負調控(Pillai, 2005)，故差異變化趨勢相反的miRNA-mRNA關系有重要意義。因此，本研究對東方蜜蜂微孢子蟲上調miRNA-意大利蜜蜂工蜂下調mRNA、東方蜜蜂微孢子蟲下調miRNA-意大利蜜蜂工蜂上調mRNA關系進行重點分析和探討。采用FPKM(fragments per kilobase of transcript per million fragments mapped)法計算和歸一化基因表達量。利用edgeR軟件(Robinsonetal., 2010)篩選AmCK1vsAmT1和AmCK2vsAmT2比較組的顯著性DEmRNA，篩選標準為：|log2fold change|≥1且P≤0.05。

利用OmicShare云平臺(www.omicshare.com)的GO分析軟件和KEGG通路分析軟件對1.4節的靶mRNA和上述靶DEmRNA進行GO數據庫(George, 2008)和KEGG(Kanehisa and Goto, 2006)通路注釋，采用默認參數。

1.6 東方蜜蜂微孢子蟲顯著性DEmiRNA與意大利蜜蜂工蜂中腸DEmRNA的調控網絡構建及分析

根據1.4和1.5節預測出的東方蜜蜂微孢子蟲顯著性DEmiRNA與意大利蜜蜂工蜂中腸顯著性DEmRNA的靶向關系，構建NcCKvsNcT1比較組中DEmiRNA與AmCK1vsAmT1比較組中DEmRNA以及NcCKvsNcT2比較組中DEmiRNA與AmCK2vsAmT2比較組中DEmRNA的調控網絡，并利用Cytoscape軟件(Smootetal., 2011)進行可視化，采用默認參數。

2 結果

2.1 東方蜜蜂微孢子蟲DEmiRNA靶向意大利蜜蜂工蜂中腸mRNA和DEmRNA的預測及分析

NcCKvsNcT1比較組中東方蜜蜂微孢子蟲的165條顯著性DEmiRNA整體上可靶向意大利蜜蜂工蜂中腸的11 711條mRNA，其中東方蜜蜂微孢子蟲的91條顯著上調和55條顯著下調miRNA分別靶向意大利蜜蜂工蜂中腸AmCK1vsAmT1的9 920和8 391條mRNA，其余19條DEmiRNA與意大利蜜蜂工蜂中腸mRNA之間無靶向關系。進一步分析發現東方蜜蜂微孢子蟲的77條顯著上調miRNA和52條顯著下調miRNA可分別靶向意大利蜜蜂工蜂中腸的118條顯著下調mRNA和135條顯著上調mRNA(圖1)。

NcCKvsNcT2比較組中東方蜜蜂微孢子蟲的124條顯著性DEmiRNA整體上可靶向意大利蜜蜂工蜂中腸的11 295條mRNA，其中東方蜜蜂微孢子蟲的56條顯著上調和49條顯著下調miRNA可分別靶向AmCK2vsAmT2比較組中意大利蜜蜂工蜂中腸的8 138和9 055條mRNA，其余19條東方蜜蜂微孢子蟲DEmiRNA與意大利蜜蜂工蜂中腸mRNA之間無靶向關系。進一步分析結果顯示，東方蜜蜂微孢子蟲的52條顯著上調miRNA可靶向AmCK2vsAmT2中意大利蜜蜂工蜂中腸的97條顯著下調mRNA，東方蜜蜂微孢子蟲的49條顯著下調miRNA可靶向意大利蜜蜂工蜂中腸的210條顯著上調mRNA(圖2)。

圖1 NcCK vs NcT1比較組中東方蜜蜂微孢子蟲的顯著性DEmiRNA與AmCK1 vs AmT1比較組中意大利蜜蜂工蜂中腸的顯著性DEmRNA之間的調控網絡

圖2 NcCK vs NcT2比較組中東方蜜蜂微孢子蟲的顯著性DEmiRNA與AmCK2 vs AmT2比較組中意大利蜜蜂工蜂中腸的顯著性DEmRNA之間的調控網絡

進一步分析發現，NcCKvsNcT1和NcCKvsNcT2中的11條共有的顯著上調miRNA(miR-3968-y, miR-5119-y和miR-317-y等)可靶向AmCK1vsAmT1和AmCK2vsAmT2中的6條皆顯著下調的mRNA (GenBank登錄號: XM_006564187.2, XM_006567501.2和XM_006558063.2等)；19條共有的顯著下調miRNA(novel-m0016-3p, let-7-x和miR-122-x等)靶向AmCK1vsAmT1和AmCK2vsAmT2中14條皆顯著上調的mRNA (GenBank登錄號: XM_016912123.1, XM_006558673.2和XM_016916036.1等)。

NcCKvsNcT1中共63條特有的DEmiRNA，其中36條顯著上調miRNA和18條顯著下調miRNA可分別靶向AmCK1vsAmT1中的97條顯著下調mRNA和91條顯著上調mRNA，其余9條東方蜜蜂微孢子蟲DEmiRNA與意大利蜜蜂工蜂中腸mRNA之間不存在靶向關系。NcCKvsNcT2中共22條特有的DEmiRNA，其中8條顯著上調miRNA和14條顯著下調miRNA可分別靶向AmCK2vsAmT2中51條顯著下調mRNA和139條顯著上調mRNA。

2.2 東方蜜蜂微孢子蟲顯著性DEmiRNA靶向意大利蜜蜂工蜂中腸顯著性DEmRNA的GO數據庫注釋

NcCKvsNcT1中的顯著上調miRNA靶向AmCK1vsAmT1中意大利蜜蜂工蜂中腸的顯著下調mRNA涉及細胞(25條下調mRNA)、細胞組分(25條下調mRNA)和細胞膜(24條下調mRNA)等10個細胞組分相關條目，細胞進程(75條下調mRNA)、代謝進程(60條下調mRNA)和單一組織進程(49條下調mRNA)等15個生物學進程相關條目，以及結合(59條下調mRNA)、催化活性(55條下調mRNA)和轉運活性(11條下調mRNA)等6個分子功能相關條目(圖3: A)。東方蜜蜂微孢子蟲的顯著下調miRNA靶向的意大利蜜蜂工蜂中腸顯著上調mRNA涉及細胞膜(15條上調mRNA)、細胞膜組分(15條上調mRNA)和細胞(12條上調mRNA)等7個細胞組分相關條目，細胞進程(34條上調mRNA)、代謝進程(30條上調mRNA)和單一組織進程(21條上調mRNA)等10個生物學進程相關條目，以及結合(34條上調mRNA)、催化活性(24條上調mRNA)和核酸結合轉錄因子活性(5條上調mRNA)等8個分子功能相關條目(圖3: B)。

圖3 NcCK vs NcT1比較組中東方蜜蜂微孢子蟲顯著性DEmiRNA靶向AmCK1 vs AmT1比較組中意大利蜜蜂工蜂中腸顯著性DEmRNA的GO數據庫注釋

NcCKvsNcT2中東方蜜蜂微孢子蟲的顯著上調miRNA靶向AmCK2vsAmT2中意大利蜜蜂工蜂中腸的顯著下調mRNA涉及細胞膜組分(20條下調mRNA)、細胞膜(20條下調mRNA)和細胞(14條下調mRNA)等11個細胞組分相關條目，細胞進程(52條下調mRNA)、代謝進程(42條下調mRNA)和單一組織進程(35條下調mRNA)等10個生物學進程相關條目，以及結合(49條下調mRNA)、催化活性(41條下調mRNA)和轉運活性(9條下調mRNA)等6個分子功能相關條目(圖4: A)。東方蜜蜂微孢子蟲的顯著下調miRNA靶向的意大利蜜蜂工蜂中腸顯著上調mRNA涉及細胞(24條上調mRNA)、細胞組分(24條上調mRNA)和細胞膜組分(22條上調mRNA)等9個細胞組分相關條目，細胞進程(59條上調mRNA)、代謝進程(47條上調mRNA)和單一組織進程(42條上調mRNA)等12個生物學進程相關條目，以及結合(57條上調mRNA)、催化活性(35條上調mRNA)和轉運活性(14條上調mRNA)等9個分子功能相關條目(圖4: B)。

此外，NcCKvsNcT1和NcCKvsNcT2中共同的顯著上調miRNA靶向AmCK1vsAmT1和AmCK2vsAmT2中共同的顯著下調mRNA可注釋到單一組織進程(2條下調mRNA)、細胞進程(2條下調mRNA)和結合(2條下調mRNA)等7個條目；NcCKvsNcT1和NcCKvsNcT2中共同的顯著下調miRNA靶向AmCK1vsAmT1和AmCK2vsAmT2中共同的顯著上調mRNA可注釋到結合(3條上調mRNA)、代謝進程(2條上調mRNA)和細胞(2條上調mRNA)等10個條目。

NcCKvsNcT1中特有的顯著上調miRNA靶向AmCK1vsAmT1中顯著下調mRNA涉及26個條目，包括細胞膜組分(9條下調mRNA)、細胞膜(9條下調mRNA)和細胞(8條下調mRNA)等7個細胞組分相關條目，細胞進程(27條下調mRNA)、單一組織進程(21條下調mRNA)和代謝進程(21條下調mRNA)等13個生物學進程相關條目，結合(24條下調mRNA)、催化活性(19條下調mRNA)、分子轉導活性(5條下調mRNA)等6個分子功能相關的條目。NcCKvsNcT1特有的東方蜜蜂微孢子蟲顯著下調miRNA靶向AmCK1vsAmT1中顯著上調mRNA涉及22個功能條目，包括9個生物學進程相關條目，7個細胞組分相關條目，以及6個分子功能相關條目。NcCKvsNcT2中特有的顯著上調miRNA靶向AmCK2vsAmT2中顯著下調mRNA可注釋到22個功能條目，包括細胞(10條下調mRNA)、細胞組分(10條下調mRNA)和細胞器(9條下調mRNA)等7個細胞組分相關條目，細胞進程(19條下調mRNA)、代謝進程(16條下調mRNA)和單一組織進程(10條下調mRNA)等8個生物學進程相關條目，結合(24條下調mRNA)、催化活性(13條下調mRNA)和核酸結合轉錄因子活性(4條下調mRNA)等7個分子功能相關條目；NcCKvsNcT2中特有的東方蜜蜂微孢子蟲顯著下調miRNA靶向AmCK2vsAmT2中意大利蜜蜂工蜂顯著上調mRNA涉及28個功能條目，包括9個細胞組分相關條目，11個生物學進程相關條目，以及8個分子功能相關條目。

2.3 東方蜜蜂微孢子蟲顯著性DEmiRNA靶向意大利蜜蜂工蜂中腸顯著性DEmRNA的KEGG通路注釋

NcCKvsNcT1中東方蜜蜂微孢子蟲的顯著上調miRNA靶向AmCK1vsAmT1中意大利蜜蜂工蜂中腸的顯著下調mRNA可注釋到剪接體(4條下調mRNA)、光傳導-果蠅(3條下調mRNA)和磷酸肌醇代謝(3條下調mRNA)等113條通路，其中包括胞吞作用(2條下調mRNA)等細胞免疫通路，以及泛素介導的蛋白水解(2條下調mRNA)、黑色素生成(3條下調mRNA)等體液免疫通路(圖5: A)。NcCKvsNcT1中東方蜜蜂微孢子蟲的顯著下調miRNA靶向意大利蜜蜂工蜂中腸的顯著上調mRNA可注釋到Hedgehog信號通路-果蠅(3條上調mRNA)、粘著斑(3條上調mRNA)和生理節律-果蠅(2條上調mRNA)等107條通路，其中包括胞吞作用(4條上調mRNA)、泛素介導的蛋白水解(1條上調mRNA)等細胞免疫通路，黑色素生成(1條上調mRNA)等體液免疫通路，以及色氨酸代謝(1條上調mRNA)等3條氨基酸代謝通路，檸檬酸循環(TCA循環)(2條上調mRNA)等4條碳水化合物代謝通路和氧化磷酸化(1條上調mRNA)等2條能量代謝通路(圖5: B)。注釋數量前20位的通路信息詳見圖5。

圖5 NcCK vs NcT1比較組中東方蜜蜂微孢子蟲顯著性DEmiRNA靶向AmCK1 vs AmT1比較組中意大利蜜蜂工蜂中腸顯著性DEmRNA注釋的前20位KEGG通路

NcCKvsNcT2中東方蜜蜂微孢子蟲的顯著上調miRNA靶向AmCK2vsAmT2中意大利蜜蜂工蜂中腸的顯著下調mRNA涉及代謝通路(11條下調mRNA)、鞘脂信號通路(5條下調mRNA)及硫代謝(1條下調mRNA)等97條通路，其中包括溶酶體(1條下調mRNA)、黑色素生成(1條下調mRNA)、胞吞作用(1條下調mRNA)等細胞免疫通路。NcCKvsNcT2中東方蜜蜂微孢子蟲的顯著下調miRNA靶向意大利蜜蜂工蜂中腸的顯著上調mRNA涉及代謝通路(9條上調mRNA)、Hippo信號通路(5條上調mRNA)和其他聚糖降解(4條上調mRNA)等127條通路，其中包括溶酶體(4條上調mRNA)、胞吞作用(2條上調mRNA)等細胞免疫通路，以及Toll/Imd信號通路(1條上調mRNA)、Toll樣受體信號通路(1條上調mRNA)等體液免疫通路，嘧啶代謝(1條上調mRNA)等2條核苷酸代謝通路，角質、亞黃素和蠟的生物合成(1條上調mRNA)等3條脂質代謝通路，谷胱甘肽代謝(1條上調mRNA)等2條其他氨基酸的代謝通路，1條萜類和聚酮類化合物的代謝通路。注釋數量前20位的通路信息詳見圖6。

圖6 NcCK vs NcT2比較組中東方蜜蜂微孢子蟲顯著性DEmiRNA靶向AmCK2 vs AmT2比較組中意大利蜜蜂工蜂中腸顯著性DEmRNA注釋的前20位KEGG通路

此外，NcCKvsNcT1和NcCKvsNcT2中共同的顯著上調miRNA靶向AmCK1vsAmT1和AmCK2vsAmT2中共同的顯著下調mRNA未注釋到任何通路；共同的東方蜜蜂微孢子蟲顯著下調miRNA靶向共同的意大利蜜蜂工蜂中腸顯著上調mRNA可注釋到9條通路，包括生理節律(1條上調mRNA)、生理節律-果蠅(1條上調mRNA)、Hippo信號通路-多物種(1條上調mRNA)、Hedgehog信號通路(1條上調mRNA)、Hedgehog信號通路-果蠅(1條上調mRNA)、FoxO信號通路(1條上調mRNA)、Hippo信號通路(1條上調mRNA)、Hippo信號通路-果蠅(1條上調mRNA)和Wnt信號通路(1條上調mRNA)等。

NcCKvsNcT1中特有的顯著上調miRNA靶向AmCK1vsAmT1中特有的顯著下調mRNA涉及代謝途徑(5條下調mRNA)、剪接體(4條下調mRNA)和光傳導-果蠅(3條下調mRNA)等106條通路。NcCKvsNcT1東方蜜蜂微孢子蟲特有的顯著下調miRNA靶向意大利蜜蜂工蜂中腸特有的顯著上調mRNA涉及代謝途徑(4條上調mRNA)、Hippo信號通路(3條上調mRNA)和生理節律-果蠅(2條上調mRNA)等92條通路。

NcCKvsNcT2中特有的顯著上調miRNA靶向AmCK2vsAmT2中顯著下調mRNA涉及代謝通路(5條下調mRNA)、糖胺聚糖降解(1條下調mRNA)及Notch信號通路(1條下調mRNA)等36條通路。NcCKvsNcT2中東方蜜蜂微孢子蟲特有的顯著下調miRNA靶向意大利蜜蜂工蜂中腸特有的顯著上調mRNA涉及各種類型的N-聚糖的生物合成(4條上調mRNA)、溶酶體(4條上調mRNA)及cAMP信號通路(4條上調mRNA)等112條通路。

2.4 東方蜜蜂微孢子蟲顯著性DEmiRNA及其靶向意大利蜜蜂工蜂中腸的免疫防御相關DEmRNA的調控網絡

根據2.3節通路注釋的結果，東方蜜蜂微孢子蟲顯著性DEmiRNA靶向意大利蜜蜂工蜂中腸DEmRNA涉及的免疫防御相關通路包括：胞吞作用、黑色素生成、溶酶體、自噬、Toll樣受體信號通路、細胞凋亡、Ras信號通路、泛素介導的蛋白水解和MAPK信號通路。

圖7 東方蜜蜂微孢子蟲的顯著性DEmiRNA靶向意大利蜜蜂工蜂中腸免疫防御和能量代謝相關通路注釋的顯著性DEmRNA的調控網絡

NcCKvsNcT1中顯著上調的miR-14-y [log2(fold change)=16.41,P=3.76E-05]靶向意大利蜜蜂工蜂中腸的mRNA (GenBank登錄號: XM_395448.6) [log2(fold change)=-10.02,P=0.0076]； miR-283-x [log2(fold change)=5.50,P=0.049]和miR-216-x [log2(fold change)=16.33,P=1.43E-07]共同靶向意大利蜜蜂工蜂中腸的E3泛素蛋白連接酶Nedd-4亞型X8編碼基因的mRNA (GenBank登錄號: XM_016917135.1) [log2(fold change)=-10.57,P=0.0020]。NcCKvsNcT1中顯著下調的miR-181-x[log2(fold change)=-2.22,P=0.0051]和miR-128-y[log2(fold change)=-1.87,P=5.55E-04]共同靶向的意大利蜜蜂工蜂中腸mRNA (GenBank登錄號: XM_016915020.1) [log2(fold change)=1.04,P=0.019], miR-340-x[log2(fold change)=-1.12,P=0.024]靶向的mRNA(GenBank登錄號: XM_006570195.2)[log2(fold change)=9.52,P=0.043], miR-484-z[log2(fold change)=-1.64,P=0.029]靶向的mRNA (GenBank登錄號: XM_006571076.2)[log2(fold change)=9.40,P=0.033]和miR-320-y[log2(fold change)=-4.73,P=1.18E-04]靶向的mRNA (GenBank登錄號: XM_394836.6) [log2(fold change)=4.01,P=0.0086]均可注釋到胞吞作用(圖7: A)。顯著上調的miR-3202-x[log2(fold change)=12.47,P=1.75E-07], miR-317-y[log2(fold change)=19.06,P=1.72E-06]和miR-252-x[log2(fold change)=14.99,P=9.90E-05]等6條東方蜜蜂微孢子蟲DEmiRNA共同靶向的意大利蜜蜂工蜂中腸mRNA (GenBank登錄號: XM_016911563.1) [log2(fold change)=-8.28,P=0.0093]， miR-283-x, miR-8-y [log2(fold change)=7.03,P=0.032]和miR-216-x等7條DEmiRNA共同靶向的1-磷脂酰肌醇-4,5-二磷酸磷酸二酯酶I類和II類亞型X1編碼基因mRNA (GenBank登錄號: XM_003250942.3) [log2(fold change)=-3.41,P=0.0084]，miR-339-y [log2(fold change)=12.47,P=1.75E-07], miR-9226-y [log2(fold change)=21.32,P=1.98E-10]和miR-317-y 4條DEmiRNA共同靶向mRNA (GenBank登錄號: XM_006559228.1) [log2(fold change)=-2.45,P=0.045]以及顯著下調的novel-m0016-3p [log2(fold change)=-7.07,P=1.82E-05]和novel-m0028-5p [log2(fold change)=-14.19,P=3.86E-07]共同靶向mRNA (GenBank登錄號: XM_006567354.2) [log2(fold change)=2.08,P=9.21E-04]均可注釋到黑色素生成(圖7: A)。顯著上調的miR-14-y靶向的意大利蜜蜂工蜂中腸mRNA (GenBank登錄號: XM_395448.6) [log2(fold change)=-10.02,P=0.0076], miR-283-x和miR-216-x共同靶向的mRNA (GenBank登錄號: XM_016917135.1)以及顯著下調的miR-106-x [log2(fold change)=-2.30,P=0.026], miR-17-x [log2(fold change)=-1.26,P=0.036], let-7-y [log2(fold change)=-3.76,P=0.0070], miR-20-x [log2(fold change)=-3.70,P=9.29E-04)和miR-146-x [log2(fold change)=-2.86,P=0.037)共同靶向的mRNA (GenBank登錄號: XM_006564220.2) [log2(fold change)=10.41,P=0.010]均可注釋到泛素介導的蛋白水解(圖7: A)。顯著上調的miR-283-x, miR-216-x和miR-4796-y [log2(fold change)=15.51,P=1.21E-06]共同靶向的mRNA (GenBank登錄號: XM_016914371.1)，顯著下調的miR-126-y [log2(fold change)=-3.77,P=0.0098]和miR-128-y共同靶向的mRNA (GenBank登錄號: XM_394432.6) [log2(fold change)=7.97,P=0.035]均可注釋到MAPK信號通路(圖7: A)。

在NcCKvsNcT2中，東方蜜蜂微孢子蟲顯著下調的novel-m0016-3p [log2(fold change)=-4.37,P=3.24E-04]靶向意大利蜜蜂工蜂中腸的mRNA (GenBank登錄號: XM_016917619.1) [log2(fold change)=12.12,P=4.14E-05]和miR-320-y [log2(fold change)=-4.87,P=0.0016]靶向的mRNA (GenBank登錄號: XM_006572368.2) [log2(fold change)=8.59,P=0.0041]以及顯著上調的miR-5106-y [log2(fold change)=19.69,P=9.27E-09]靶向的mRNA (GenBank登錄號: XM_006561133.2)[log2(fold change)=-1.13,P=0.045]皆可注釋到胞吞作用(圖7: B)。顯著下調的novel-m0021-5p [log2(fold change)=-1.52,P=0.036]靶向意大利蜜蜂工蜂中腸的4條mRNAs (GenBank登錄號: XM_006567428.2) [log2(fold change)=1.52,P=0.0042]，(GenBank登錄號: XM_006567431.2) [log2(fold change)=1.34,P=0.039]，(GenBank登錄號: XM_016914010.1) [log2(fold change)=1.22,P=1.47E-04]和(GenBank登錄號: XM_016914011.1) [log2(fold change)=3.09,P=0.042]以及顯著上調的novel-m0022-5p [log2(fold change)=14.13,P=7.97E-06]和miR-216-x [log2(fold change)=16.33,P=1.43E-07]共同靶向的類α-乙酰氨基葡萄糖苷酶編碼基因mRNA (GenBank登錄號: XM_016913157.1) [log2(fold change)=-9.33,P=0.010]皆可注釋到溶酶體(圖7: B)。顯著下調的novel-m0002-3p [log2(fold change)=-5.85,P=6.42E-05], novel-m0021-5p和miR-374-x[log2(fold change)=-15.48,P=1.23E-07]共同靶向的mRNA (GenBank登錄號: XM_016911751.1) [log2(fold change)=1.65,P=0.0040]以及miR-26-x [log2(fold change)=-1.81,P=0.012]和miR-486-x [log2(fold change)=-8.86,P=0.028]共同靶向的mRNA (GenBank登錄號: XM_006557238.2) [log2(fold change)=3.73,P=0.0047]皆可注釋到自噬-動物(圖7: B)。顯著上調的bantam-y [log2(fold change)=16.57,P=1.84E-06], miR-8-y [log2(fold change)=5.53,P=0.044], miR-216-x和miR-5119-y [log2(fold change)=9.19,P=0.022]共同靶向的1-磷脂酰肌醇-4,5-二磷酸磷酸二酯酶I類和II類亞型X3編碼基因mRNA (GenBank登錄號: XM_006567225.2) [log2(fold change)=-8.75,P=0.016]可注釋到黑色素生成(圖7: B)。miR-424-x [log2(fold change)=13.86,P=9.14E-06], miR-7465-x [log2(fold change)=15.36,P=4.75E-07]和miR-547-x [log2(fold change)=13.98,P=5.96E-08]共同靶向的mRNA (GenBank登錄號: XM_006564132.2) [log2(fold change)=-8.92,P=0.038]可注釋到MAPK信號通路-果蠅。顯著下調的miR-26-x和miR-486-x共同靶向的mRNA (GenBank登錄號: XM_006557238.2)可同時注釋到細胞凋亡、MAPK信號通路和Toll樣受體信號通路(圖7: B)。顯著上調的miR-5106-y靶向的mRNA (GenBank登錄號: XM_006561133.2),顯著下調的miR-26-x和miR-486-x共同靶向的mRNA (GenBank登錄號: XM_006557238.2), novel-m0021-5p, miR-221-z [log2(fold change)=-2.89,P=0.0039], miR-374-x, miR-423-y [log2(fold change)=-2.04,P=0.012], miR-17-x [log2(fold change)=-1.92,P=0.039]和miR-16-y [log2(fold change)=-16.85,P=5.38E-08]共同靶向的mRNA (GenBank登錄號: XM_006569471.2) [log2(fold change)=5.64,P=0.027]以及顯著上調的miR-8-y, bantam-y和novel-m0022-5p [log2(fold change)=14.13,P=7.97E-06]等10條miRNA共同靶向的1-磷脂酰肌醇-4,5-二磷酸磷酸二酯酶ε-1編碼基因mRNA (GenBank登錄號: XM_016911813.1) [log2(fold change)=-7.89，P=0.023]皆可注釋到Ras信號通路(圖7: B)。

此外，根據2.3節通路注釋的結果，東方蜜蜂微孢子蟲顯著性DEmiRNA靶向的意大利蜜蜂工蜂中腸DEmRNA可注釋到氧化磷酸化和硫代謝等能量代謝相關通路。NcCKvsNcT1中顯著下調的miR-584-x [log2(fold change)=-1.32,P=0.042]和miR-151-y [log2(fold change)=-2.21,P=6.92E-04]等2條miRNA共同靶向的細胞色素c氧化酶亞基4亞型X1，線粒體亞型X1編碼基因mRNA (GenBank登錄號: XM_006572254.2) [log2(fold change)=-2.21,P=6.92E-04],顯著上調的miR-459-x[log2(fold change)=16.24,P=1.76E-07], miR-315-x[log2(fold change)=12.80,P=4.70E-06]和miR-5119-y分別靶向的mRNAs (GenBank登錄號: XM_001122063.4)[log2(fold change)=-1.07,P=0.022]、(GenBank登錄號: XM_006567193.2)[log2(fold change)=-1.02,P=8.31E-04]和(GenBank登錄號: XM_623229.5) [log2(fold change)=-2.86,P=0.0095]皆可注釋到氧化磷酸化(圖7: C)。NcCKvsNcT2中顯著上調的miR-301-y [log2(fold change)=15.58,P=6.07E-09]靶向的mRNA (GenBank登錄號: XM_396499.6) [log2(fold change)=-2.09,P=0.0048]可注釋到硫代謝(圖7: D)。

NcCKvsNcT1中特有的miR-186-x [log2(fold change)=1.02,P=0.047767213]靶向意大利蜜蜂工蜂中腸的mRNAs (GenBank登錄號: XM_006559228.1)和(GenBank登錄號: NM_001011611.2) [log2(fold change)=-4.44,P=0.0048]可分別注釋到黑色素生成通路和免疫系統進程條目，NcCKvsNcT2中特有的novel-m0022-5p靶向意大利蜜蜂工蜂中腸的mRNA (GenBank登錄號: XM_006559395.1) [log2(fold change)=-1.23,P=0.0088], mRNA(GenBank登錄號: XM_006560586.2) [log2(fold change)=-1.13,P=1.03E-07]和mRNA(GenBank登錄號: XM_016911813.1)可注釋到代謝途徑。

NcCKvsNcT1中特有的miR-22-x [log2(fold change)=-13.99,P=6.33E-06]靶向意大利蜜蜂工蜂中腸的mRNAs (GenBank登錄號: XM_006563625.2) [log2(fold change)=8.42,P=0.027]和(GenBank登錄號: XM_006559419.2) [log2(fold change)=1.30,P=0.0078]可分別注釋到生物附著和應激反應條目。NcCKvsNcT2中特有的miR-26-x靶向意大利蜜蜂工蜂中腸的mRNA (GenBank登錄號: XM_006557238.2)可同時注釋到Toll樣受體信號通路、Hippo信號通路、MAPK信號通路及生物進程調節功能條目。

3 討論

3.1 東方蜜蜂微孢子蟲DEmiRNA對意大利蜜蜂工蜂中腸的基因表達具有潛在的廣泛影響

本研究中，NcCKvsNcT1中東方蜜蜂微孢子蟲的顯著下調miRNA靶向的意大利蜜蜂工蜂中腸顯著上調mRNA涉及7個細胞組分相關條目，10個生物學進程相關條目，8個分子功能相關條目(圖3: B)；以及3條氨基酸代謝通路，4條碳水化合物代謝通路，2條能量代謝通路(圖5: B)。NcCKvsNcT2中東方蜜蜂微孢子蟲的顯著下調miRNA靶向意大利蜜蜂工蜂中腸的顯著上調mRNA涉及9個細胞組分相關條目，12個生物學進程相關條目，9個分子功能相關條目(圖4: B)；以及2條核苷酸代謝通路，3條脂質代謝通路，2條其他氨基酸的代謝通路，1條萜類和聚酮類化合物的代謝通路(圖6: B)。這些結果暗示東方蜜蜂微孢子蟲在侵染意大利蜜蜂工蜂中腸的不同時間點通過差異表達部分miRNA對相應的意大利蜜蜂工蜂基因表達進行調控，進而影響上述與意大利蜜蜂工蜂物質和能量代謝相關的多個方面。

此外，NcCKvsNcT1和NcCKvsNcT2比較組共同的顯著下調miRNA靶向AmCK1vsAmT1和AmCK2vsAmT2比較組共同的顯著上調mRNA可注釋到10個GO條目，以及9條KEGG通路。上述結果表明東方蜜蜂微孢子蟲在侵染意大利蜜蜂工蜂中腸的過程中可能通過持續差異表達一些miRNA對宿主的細胞活動、生化反應、生理活動和信號轉導產生廣泛的影響。

3.2 東方蜜蜂微孢子蟲可能通過下調表達部分miRNA對意大利蜜蜂工蜂中腸的能量代謝進行調控

在長期的協同進化中，東方蜜蜂微孢子蟲的基因組大大簡化，多數物質和能量代謝的途徑丟失，其增殖高度依賴宿主細胞提供物質和能量(Burrietal., 2006)。Chen等(2013)研究發現東方蜜蜂微孢子蟲侵染可對西方蜜蜂工蜂造成較強的能量脅迫，致使宿主的能量代謝速率加快且糖水的攝入量提高。ATP/ADP轉移酶和ABC轉運蛋白參與微孢子蟲對宿主物質和能量的竊取(Paldietal., 2010)。前期研究中，筆者所在課題組發現東方蜜蜂微孢子蟲的7個ABC轉運蛋白編碼基因在NcCKvsNcT1和NcCKvsNcT2中差異表達(耿四海等, 2020b)。本研究發現，NcCKvsNcT1中miR-151-y和miR-584-x顯著下調表達，且能靶向意大利蜜蜂工蜂中腸的mRNA (GenBank登錄號: XM_006572254.2)，該基因可注釋到能量代謝相關的氧化磷酸化信號通路。上述結果暗示東方蜜蜂微孢子蟲的DEmiRNA參與意大利蜜蜂工蜂中腸的關鍵能量代謝過程，控制宿主能量代謝增強以產生更多的ATP，從而掠奪更多的ATP利于自身增殖。

3.3 東方蜜蜂微孢子蟲可能通過上調表達部分miRNA調控意大利蜜蜂工蜂中腸的免疫防御

昆蟲的天然免疫系統包括細胞免疫和體液免疫，前者主要是指東方蜜蜂微孢子蟲體被血淋巴細胞等識別、包被和吞噬，而后者主要包括各種識別受體、胞內信號因子、抗微生物肽、蛋白酶以及酶級聯反應導致血淋巴凝結、黑化等(Leulieretal., 2000)。較多的研究證明，東方蜜蜂微孢子蟲在侵染宿主的過程中可通過多種方式跨界調控宿主的基因表達，并抑制宿主免疫防御(Wangetal., 2017; Shahidetal., 2018; Cuietal., 2019; Huangetal., 2019)。菟絲子將Bc-siR3.2裝載到AGO1蛋白，靶向擬南芥Arabidopsisthaliana的絲裂原活化蛋白激酶2基因(MPK2)和MPK1基因的轉錄本，從而抑制宿主的免疫反應(Shahidetal., 2018)。球孢白僵菌將bba-milR-1載入囊泡并轉運到斯氏按蚊Anophelesstephensi的細胞內，進而跨界調控宿主基因CLIPB9和Spz4的表達，抑制Toll樣受體信號通路和逃避黑化反應(Cuietal., 2019)。本研究發現，NcCKvsNcT1和NcCKvsNcT2比較組中miR-216-x, miR-8-y, miR-4796-y, bantam-y和miR-5119-y等多個共同顯著上調的miRNA可靶向注釋到黑色素生成、MAPK信號通路、Ras信號通路、泛素介導的蛋白水解和溶酶體等通路的多條顯著下調mRNA。其中，NcCKvsNcT1中的miR-216-x同時靶向AmCK1vsAmT1中可注釋到宿主免疫相關的黑色素生成、MAPK信號通路及泛素介導的蛋白水解通路的mRNAs (GenBank登錄號分別為XM_003250942.3, XM_016914371.1和XM_016917135.1)(圖7: A)；NcCKvsNcT2中的miR-216-x還可同時靶向AmCK2vsAmT2中可注釋到意大利蜜蜂工蜂免疫相關的黑色素生成、溶酶體及Ras信號通路通路的mRNAs (GenBank登錄號分別為XM_006567225.2, XM_016913157.1和XM_016911813.1)(圖7: B)。NcCKvsNcT1中miR-8-y靶向AmCK1vsAmT1中的mRNA可注釋到黑色素生成通路(圖7: A)；NcCKvsNcT2中miR-8-y靶向AmCK2vsAmT2中的2個mRNA可分別注釋到黑色素生成和Ras信號通路(圖7: B)。NcCKvsNcT1中miR-4796-y靶向AmCK1vsAmT1中的mRNA可注釋到MAPK信號通路(圖7: A)；NcCKvsNcT2中miR-4796-y靶向AmCK2vsAmT2中的mRNA可注釋到Ras信號通路(圖7: B)。NcCKvsNcT1中bantam-y靶向AmCK1vsAmT1中的mRNA可注釋到黑色素生成(圖7: A)；NcCKvsNcT2中bantam-y同時靶向AmCK2vsAmT2中2條mRNA可分別注釋到Ras信號通路和黑色素生成(圖7: B)。NcCKvsNcT1中miR-5119-y靶向AmCK1vsAmT1中的mRNA可注釋到黑色素生成(圖7: A)；NcCKvsNcT2中miR-5119-y還同時靶向AmCK2vsAmT2中的2條mRNA可分別注釋到黑色素生成和MAPK信號通路(圖7: B)。

值得注意的是，筆者所在課題組前期研究發現，NcCKvsNcT2中bantam-y可靶向東方蜜蜂微孢子蟲的己糖激酶編碼基因mRNA (GenBank登錄號: XM_002995838.1)，該基因可注釋到東方蜜蜂微孢子蟲能量代謝相關的糖酵解/糖異生途徑；并且bantam-y和miR-5119-y可分別靶向東方蜜蜂微孢子蟲自身的ABC轉運蛋白編碼基因mRNA (GenBank登錄號分別為XM_002995835.1和XM_002996675.1)(耿四海等, 2020a)。據此猜測東方蜜蜂微孢子蟲在侵染意大利蜜蜂工蜂中腸的過程中通過不同程度地差異表達miR-216-x和miR-8-y等miRNA對意大利蜜蜂工蜂中腸的黑色素生成、溶酶體和Ras信號通路進行調控，此外東方蜜蜂微孢子蟲還能通過差異表達bantam-y和miR-5119-y等miRNA一方面調控意大利蜜蜂工蜂中腸的黑色素生成、Ras信號通路和MAPK信號通路，另一方面調節自身的物質和能量運輸，從而促進增殖，體現了miRNA作為具有多面性的關鍵調控因子在東方蜜蜂微孢子蟲和意大利蜜蜂工蜂中腸的互作中扮演著特殊角色。

3.4 東方蜜蜂微孢子蟲在不同侵染時間點通過特異性差異表達部分miRNA對意大利蜜蜂工蜂中腸基因進行表達調控

NcCKvsNcT1中特有的miR-186-x靶向意大利蜜蜂工蜂中腸的mRNAs (GenBank登錄號分別為XM_006559228.1和NM_001011611.2)可分別注釋到黑色素生成通路和免疫系統進程條目，NcCKvsNcT2中特有的novel-m0022-5p靶向意大利蜜蜂工蜂可注釋到代謝途徑的mRNAs (GenBank登錄號分別為XM_006559395.1, XM_006560586.2和XM_016911813.1)。我們推測東方蜜蜂微孢子蟲在感染意大利蜜蜂工蜂7 d通過上調miR-186-x抑制意大利蜜蜂工蜂中腸的免疫防御相關基因的表達，在感染10 d時通過上調novel-m0022-5p對意大利蜜蜂工蜂中腸的能量代謝相關基因進行抑制，從而加強對宿主細胞內能量的掠奪，供其增殖所需。此外，NcCKvsNcT1中特有的miR-22-x靶向意大利蜜蜂工蜂可分別注釋到生物附著和應激反應條目的mRNAs (GenBank登錄號分別為XM_006563625.2和XM_006559419.2)，NcCKvsNcT2中特有的miR-26-x靶向意大利蜜蜂工蜂的mRNA (GenBank登錄號: XM_006557238.2)可同時注釋到Toll樣受體信號通路、MAPK信號通路及生物進程調節功能條目。這些結果表明東方蜜蜂微孢子蟲在感染意大利蜜蜂工蜂7 d和10 d可能通過分別下調miR-22-x和miR-26-x對上述重要的功能條目和通路相關基因表達進行抑制，從而影響宿主的正常生命活動。

綜上所述，本研究結合前期獲得的高質量miRNA和mRNA組學數據對東方蜜蜂微孢子蟲DEmiRNA靶向的意大利蜜蜂工蜂中腸mRNA和DEmRNA進行預測、分析和探討，揭示病原的DEmiRNA對宿主的基因表達調控具有廣泛性，病原可能通過上調部分miRNA調控宿主的免疫防御以促進增殖，通過下調部分miRNA跨界調控宿主的能量代謝以加強能量竊取并促進增殖。研究結果為闡明東方蜜蜂微孢子蟲跨界調控意大利蜜蜂工蜂中腸的分子機制提供了理論依據和數據基礎，也為深入理解二者的互作機制提供了有價值的參考信息。