意大利蜜蜂幼蟲腸道在球囊菌侵染前期的差異表達microRNA及其調控網絡

郭睿，杜宇，童新宇，熊翠玲，鄭燕珍，徐國鈞，王海朋，耿四海，周丁丁，郭意龍，吳素珍，陳大福

?

意大利蜜蜂幼蟲腸道在球囊菌侵染前期的差異表達microRNA及其調控網絡

郭睿，杜宇，童新宇，熊翠玲，鄭燕珍，徐國鈞，王海朋，耿四海，周丁丁，郭意龍，吳素珍，陳大福

（福建農林大學蜂學學院，福州 350002）

【目的】微小RNA（microRNA，miRNA）是一類重要的基因表達調控因子，可影響宿主與病原間的互作過程。蜜蜂球囊菌（）是一種特異性侵染蜜蜂幼蟲的致死性真菌病原。本研究旨在對意大利蜜蜂（，簡稱意蜂）幼蟲腸道在球囊菌脅迫前期的差異表達miRNA（differentially expressed miRNA，DEmiRNA）及其靶基因進行深入分析，在miRNA組學水平探究意蜂幼蟲在球囊菌侵染前期的脅迫應答，并通過構建顯著DEmiRNA的調控網絡篩選出與宿主應答相關的關鍵miRNA。【方法】利用small RNA-seq（sRNA-seq）技術對正常及球囊菌侵染的意蜂4日齡幼蟲腸道（AmCK和AmT）進行高通量測序，首先對原始數據進行質控和評估，隨后將過濾后的數據與西方蜜蜂（）參考基因組進行比對；將比對上的序列標簽（tags）注釋到miRBase數據庫，得出已知miRNA的表達量；通過TPM（tags per million）算法對各樣本中miRNA的表達量進行歸一化處理，以|log2fold change|≥1且≤0.05作為標準篩選得到顯著DEmiRNA；利用TargetFinder 軟件預測顯著DEmiRNA的靶基因，并對其進行GO和KEGG代謝通路（pathway）富集分析。利用Cytoscape軟件對miRNA-mRNA調控網絡進行可視化。最后，利用莖環反轉錄PCR（Stem-loop RT-PCR）和熒光定量PCR（qPCR）驗證測序數據的可靠性。【結果】AmCK和AmT樣品的測序分別得到13 553 302和10 777 534條原始讀段（raw reads），經嚴格過濾后得到的有效讀段（clean reads）數分別為13 186 921和10 480 913條。各樣品的生物學重復間的Pearson相關性系數分別在0.9822和0.9508以上。共有10個顯著DEmiRNA，包括4個上調miRNA和6個下調miRNA。顯著DEmiRNA在AmT的整體表達水平低于AmCK。10個顯著DEmiRNA可靶向結合3 788個靶基因，其中上調miRNA的1 240個靶基因可注釋到GO數據庫中的39個GO條目，主要富集在結合、細胞進程、代謝進程和應激反應等；下調miRNA的749個靶基因可注釋到34個GO條目，主要富集在細胞進程、結合、代謝進程和應激反應等。KEGG數據庫注釋結果顯示，上調miRNA和下調miRNA的靶基因分別注釋到95和66條代謝通路，富集基因數最多的分別是Wnt信號通路、Hippo信號通路、光傳導以及內吞作用、磷脂酰肌醇信號系統、嘌呤代謝。對于上調和下調miRNA，分別有31和52個靶基因注釋到內吞作用，15和7個靶基因注釋到泛素介導的蛋白水解，11和5個靶基因注釋到Jak-STAT信號通路，1和3個靶基因注釋到MAPK信號通路。顯著DEmiRNA與靶mRNA之間形成復雜的調控網絡，7個顯著DEmiRNA靶向結合96個與Wnt信號通路相關的mRNA，8個顯著DEmiRNA靶向結合55個與內吞作用相關的mRNA。Stem-loop RT-PCR和qPCR結果驗證了測序數據的可靠性。【結論】對意蜂幼蟲腸道在球囊菌侵染前期的DEmiRNA及其靶基因進行預測和分析，并構建和分析了DEmiRNA-mRNA調控網絡，研究結果提供了宿主miRNA的表達譜和差異表達信息，揭示了DEmiRNA通過調控細胞生命活動、新陳代謝以及部分細胞和體液免疫等生物學過程參與宿主的脅迫應答。miR-4331-y、miR-4968-y、miR-8440-y、novel-m0023-5p和novel-m0025-3p共同參與了宿主的Wnt信號通路和內吞作用的調控，可作為白堊病治療的潛在分子靶標。

意大利蜜蜂；幼蟲腸道；發育；差異表達微小RNA；調控網絡

0 引言

【研究意義】意大利蜜蜂（，簡稱意蜂）是我國養蜂生產中的主要蜂種，在農業生產和生態保護等方面具有重要價值[1]。白堊病是蜜蜂球囊菌（）特異性侵染蜜蜂幼蟲而引起的致死性真菌病害，對養蜂生產造成巨大損失。蜜蜂幼蟲食入被球囊菌孢子污染的食物，孢子進入中腸后開始低水平萌發，到了預蛹期隨著中腸和后腸的連通，孢子進入后腸后在氧氣的刺激下劇烈萌發，同時菌絲大量生長并陸續穿透腸壁和體壁，從而導致幼蟲死亡[2]。目前，有關蜜蜂幼蟲響應球囊菌脅迫應答及分子調控機制的研究報道較少。微小RNA（microRNA，miRNA）是一類長度約為18—25個核苷酸的非編碼RNA（non-coding RNA，ncRNA），通過與其靶基因3′ UTR（untranslated region）特異性結合，降解或抑制mRNA的蛋白質翻譯過程[3]。作為基因表達的關鍵調控因子，miRNA可廣泛參與細胞增殖與分化[4]、生長發育[5]、細胞凋亡[6]和抗病毒[7]等生物學過程。近年來，較多的研究結果表明miRNA參與對昆蟲及其病原互作的調控[8-9]。因此，解析蜜蜂幼蟲在球囊菌侵染前期的miRNA表達譜，深入研究差異表達miRNA（differentially expressed miRNA，DEmiRNA）及其調控網絡，可為宿主響應球囊菌脅迫的分子調控及宿主-病原互作機制提供必要基礎，也能為宿主的關鍵調控miRNA的篩選和功能研究提供重要的信息和線索。【前人研究進展】ZHANG等[10]研究發現，攜帶登革病毒（DENV）的埃及伊蚊（）被沃爾巴克氏體（）侵染后，aae-miR-2940可特異性表達，并顯著抑制了與登革病毒復制相關的DNA甲基化轉移酶基因（）的表達；李盛杰[11]研究發現，果蠅（）被革蘭氏陽性菌藤黃微球菌（）侵染后，其miR-310家族和miR-958可通過抑制Toll信號通路產生的抗菌肽分子Drosomycin的表達而影響宿主的免疫應答。但相比于果蠅和蚊類等模式昆蟲，蜜蜂與病原互作相關的miRNA研究較為滯后。LOUREN?O等[12]利用革蘭氏陰性菌粘質沙雷氏菌（）和侵染意蜂工蜂，發現miRNA可通過參與體液免疫中Toll、Imd、JNK以及Jak-STAT信號通路的調控，促進抗菌肽與黑化反應的激活來應對細菌的侵染；HUANG等[13]研究發現，東方蜜蜂微孢子蟲（）可能通過影響西方蜜蜂（）代謝調控相關的miRNA表達促進自身的快速繁殖。筆者所在課題組前期已通過二代測序技術對球囊菌侵染的中華蜜蜂（，簡稱中蜂）和意蜂4日齡幼蟲腸道進行轉錄組分析，在mRNA組學水平解析了病原在侵染前期與宿主間的互作[14-15]。然而，蜜蜂幼蟲在球囊菌侵染過程的miRNA表達譜仍然缺失，miRNA在蜜蜂幼蟲的脅迫應答中的作用仍未可知。【本研究切入點】此前的相關研究主要集中在球囊菌侵染后期，即6日齡幼蟲到預蛹的過渡期，有關球囊菌侵染前期的研究還非常滯后。球囊菌在侵染前期已經開始出現低水平的孢子萌發和菌絲生長，必然伴隨著復雜的基因表達和ncRNA調控。miRNA作為關鍵的轉錄后調控因子在昆蟲的免疫防御過程中發揮關鍵的作用[16]。本研究利用small RNA-seq（sRNA-seq）技術對正常及球囊菌脅迫的意蜂4日齡幼蟲腸道進行高通量測序，對宿主miRNA的表達譜、顯著DEmiRNA及其調控網絡進行探究，并篩選免疫防御相關的關鍵miRNA。【擬解決的關鍵問題】通過深入分析顯著DEmiRNA及其調控網絡，在miRNA組學水平解析意蜂幼蟲在球囊菌侵染前期的脅迫應答，揭示顯著DEmiRNA在宿主-病原互作中的作用，為解析意蜂幼蟲-球囊菌的互作機制打下基礎。

1 材料與方法

試驗于2017年9月至2018年7月在福建農林大學蜂學學院蜜蜂保護實驗室完成。

1.1 供試生物材料

供試意蜂幼蟲取自福建農林大學蜂學學院教學蜂場。

1.2 意蜂幼蟲人工飼養及樣品測序

參照前期研究建立的方法[14]活化球囊菌、純化孢子以及人工飼養意蜂幼蟲，配制含球囊菌孢子的飼料（終濃度為1×107個孢子/mL），飼喂處理組3日齡幼蟲，24 h后待含孢子飼料被食盡，飼喂不含球囊菌的正常飼料；對照組幼蟲飼喂不含孢子的正常飼料。腸道是球囊菌與意蜂幼蟲互作的主要部位，因此選擇幼蟲腸道作為測序材料。分別解剖上述對照組和處理組4日齡幼蟲腸道組織，剖取的幼蟲腸道每9只放于1個RNA-Free的EP管中，液氮速凍后保存在-80℃的超低溫冰箱。進行3次生物學重復。對照組的3個生物學重復分別為AmCK-1、AmCK-2和AmCK-3；處理組的3個生物學重復分別為AmT-1、AmT-2和AmT-3。委托廣州基迪奧生物科技有限公司對上述6個腸道樣品進行單端測序，測序平臺為Illumina MiSeq。首先，利用Trizol法從意蜂幼蟲腸道樣本中提取total RNA，瓊脂糖凝膠電泳切膠選擇18—30 nt的片段，連接3′端接頭，連接產物以15%變性PAGE膠電泳分離，切膠選擇36—44 nt的目的條帶；將上述回收產物連接5′端接頭，進而對兩側帶有接頭的small RNA樣本進行RT-PCR，反轉錄產物以3.5%瓊脂糖凝膠電泳分離，切膠回收140—160 bp區域的目的條帶，回收產物即為終文庫，建好的文庫上機測序。原始數據已上傳NCBI SRA數據庫，BioProject號：PRJNA408312。

1.3 測序數據的質控及DEmiRNA的預測

參照前期研究建立的方法[17]對測序數據進行質控和定位，并將得到miRNA的非注釋標簽序列（unannotated tags）。通過Bowit軟件將非注釋標簽序列與西方蜜蜂的參考基因組（assembly Amel 4.5）序列進行比對，得到相關tags在參考基因組上的位置信息，即為mapped tags。利用miRDeep2軟件[18]將mapped tags與miRBase數據庫中已知miRNA前體序列進行比對，從而鑒定已知miRNA的表達。并通過TPM（tags per million）算法對各樣本中miRNA的表達量進行歸一化處理。利用R軟件計算各樣品之間的相關性系數。以|log2fold change|≥1且≤0.05為篩選顯著DEmiRNA的標準。

1.4 顯著DEmiRNA靶基因預測及分析

利用TargetFinder軟件[19]對顯著DEmiRNA的靶基因進行預測，并通過BLAST軟件將預測的靶基因序列與Gene Ontology和KEGG數據庫比對。將顯著DEmiRNA的靶基因分別按照分子功能、細胞組分、生物學進程進行分類，并利用基迪奧OmicShare云平臺（http://www.omicshare.com/tools/Home/Index/index. html）對靶基因進行KEGG pathway富集分析。最后，根據DEmiRNA與mRNA的靶向結合關系，利用Cytoscape軟件構建miRNA-mRNA的調控網絡。

1.5 DEmiRNA的莖環實時熒光定量PCR（Stem-loop RT-qPCR）驗證

通過Stem-loop RT-qPCR檢測4個隨機選取的miRNA（novel-m0031-3p、novel-m0034-5p、miR-3793-x和ame-miR-6000a-5p）在AmCK和AmT中的表達情況，以驗證sRNA-seq數據的可靠性。根據所選miRNA的序列，參照CHEN等[20]的方法，利用DNAMAN軟件（Lynnon Biosoft公司，美國）設計特異性的Stem-loop引物、上游引物和通用的反向引物，委托上海生工生物工程有限公司進行引物合成。選擇snRNA U6作為內參。利用RNA抽提試劑盒（Axygen公司，美國）分別提取AmCK和AmT的總RNA，利用Stem-loop引物進行反轉錄得到相應的cDNA作為模板進行PCR和qPCR。常規PCR體系為50 μL，包括PCR mix（TaKaRa公司，日本）25 μL、無菌水17.5 μL、正、反向引物及cDNA模板各2.5 μL。反應程序如下：95℃ 5 min，95℃ 30 s，50℃ 30 s，72℃ 1 min，34個循環，72℃ 5 min。隨后利用1.8%的瓊脂糖凝膠電泳檢測PCR的產物。qPCR反應體系（20 μL）：SYBR Green Dye 10 μL，正、反向引物各1 μL，cDNA模板1 μL，Rox 0.44 μL，DEPC水補至20 μL。在ABI 7500熒光定量PCR儀（ABI公司，美國）中進行反應，反應條件：95℃ 1 min，95℃ 15 s，49—60℃ 60 s，共40個循環，熔解曲線默認系統程序。所選miRNA的相對表達量采用2-ΔΔCt法計算。每個反應進行3次生物學重復和3次技術重復。最后利用GraphPad Prism 5軟件進行qPCR結果的檢驗及繪圖。本研究使用的引物序列信息詳見表1。

表1 本研究使用的引物

2 結果

2.1 意蜂幼蟲腸道的測序數據的質控與評估

兩組幼蟲腸道樣品的sRNA-seq分別產生13 553 302和10 777 534條raw reads，經嚴格過濾后得到的clean reads數分別為13 186 921和10 480 913條（表2）；此外，AmCK、AmT組內各生物學重復之間的Pearson相關性系數均在0.9822和0.9508以上，說明本研究的高通量測序數據質量良好，可用于進一步分析。

2.2 意蜂幼蟲腸道的DEmiRNA分析

AmCK vs AmT比較組中顯著差異和非顯著差異miRNA數為740個。其中，10個miRNA為顯著性差異表達，包含4個上調miRNA和6個下調miRNA；顯著DEmiRNA中包括3個novel miRNA和7個已知miRNA（表3）。表達譜分析結果顯示，DEmiRNA的差異變化幅度差別明顯（圖1-A），DEmiRNA在AmT中的整體表達水平低于AmCK（圖1-B）。

表2 sRNA-seq數據總覽

A：DEmiRNA的表達量聚類Expression clustering of DEmiRNAs；B：不同樣品中DEmiRNA的表達量比較Comparison of DEmiRNAs expression in different samples

表3 AmCK vs AmT 中DEmiRNA的信息統計

2.3 意蜂幼蟲腸道DEmiRNA的靶基因預測及功能注釋

對于AmCK vs AmT中的顯著DEmiRNA，共預測出3 788個靶基因。GO數據庫注釋結果顯示，上調miRNA的1 240個靶基因共涉及39個GO 條目（term），富集靶基因數最多的分別是結合（733個）、細胞進程（698個）、代謝進程（580個）、單組織進程（576個）、催化活性（439個）、細胞膜（387個）、細胞膜組件（380個）、細胞（286個）、細胞組件（286個）、定位（245個）等（圖2-A）；下調miRNA的749個靶基因可注釋到34個GO條目，富集靶基因數最多的分別是細胞進程（533個）、結合（517個）、代謝進程（394個）、單組織進程（373個）、催化活性（240個）、細胞膜（227個）、細胞膜組件（226個）、細胞（216個）、細胞組件（216個）、細胞器（192個）等（圖2-B）。進一步分析發現，對于上調和下調miRNA，分別有204和113個靶基因注釋到應激反應，說明顯著DEmiRNA參與宿主對球囊菌的脅迫應答的調控。

對顯著DEmiRNA的靶基因進行KEGG數據庫注釋，結果顯示上調miRNA的靶基因可注釋到95條代謝通路（pathway），其中富集數最多的分別是Wnt信號通路（112個）、Hippo信號通路（62個）、光傳導（55個）、嘌呤代謝（48個）、Hedgehog信號通路（43個）、神經活性的配體-受體相互作用（19個）、晝夜節律（41個）、壽命調節途徑（35個）、氨基糖與核苷酸糖代謝（32個）和內吞作用（31個）（圖3-A）；下調miRNA的靶基因注釋到66條代謝通路，其中富集數最多的分別是內吞作用（52個）、磷脂酰肌醇信號系統（30個）、嘌呤代謝（28個）、磷酸肌醇代謝（25個）、Notch信號通路（20個）、神經活性的配體-受體相互作用（19個）、葉酸合成（18個）、背腹軸的形成（14個）、Wnt信號通路（13個）和mRNA監視通路（12個）（圖3-B）。進一步分析結果顯示，在細胞免疫相關通路中，對于上調和下調的miRNA，分別有31和52個靶基因注釋到內吞作用，15和7個靶基因注釋到泛素介導的蛋白水解；在體液免疫相關通路中，對于上調和下調的miRNA，分別有11和5個靶基因注釋到Jak-STAT信號通路，有1和3個靶基因注釋到MAPK信號通路。

2.4 意蜂幼蟲腸道DEmiRNA的調控網絡構建及分析

篩選具有KEGG數據庫注釋的靶基因，并分別構建上調miRNA和下調miRNA與靶基因的調控網絡，結果顯示顯著DEmiRNA位于調控網絡的中心，其中4個上調miRNA可結合499個靶基因（圖4-A），277個靶基因可被6個下調miRNA調控（圖4-B），DEmiRNA與靶基因之間形成復雜的調控網絡。其中，上調miRNA中的miR-8440-y結合靶基因最多，達到419個（圖4-A）；下調miRNA中的novel-m0034-3p可結合149個靶基因（圖4-B）。

提取注釋到Wnt信號通路、內吞作用的靶基因，構建它們與相關顯著DEmiRNA的調控網絡，分析結果顯示7個顯著DEmiRNA結合96個注釋到Wnt信號通路的靶基因，二者形成1個較大和3個較小的調控網絡，其中結合靶基因數由多到少分別為miR-8440-y（78個）、miR-3793-x（27個）、miR-4968-y（9個）、miR-7085-x（5個）、novel-m0025-3p（3個）、novel-m0023-5p（2個）和miR-4331-y（1個）（圖5-A）；8個顯著DEmiRNA結合55個注釋到內吞作用的靶基因，二者形成1個較大的和2個較小的調控網絡，其中miR-8440-y、novel-m0034-3p、ame-miR-6000a- 5p、novel-m0025-3p、novel-m0023-5p、miR-4968-y、miR-4331-y和miR-971-y結合的靶基因數分別為26、26、12、10、5、2、1和1個（圖5-B）。進一步分析發現，miR-4331-y、miR-4968-y、miR-8440-y、novel-m0023-5p和novel-m0025-3p共同了參與上述2條代謝通路調控，而miR-3793-x和miR-7085-x僅參與對Wnt信號通路的調控，ame-miR-6000a-5p、miR-971-y和novel-m0034-3p僅涉及對內吞作用的調控（圖5-A、5-B）。

2.5 意蜂幼蟲腸道DEmiRNA的Stem-loop RT-PCR和qPCR驗證

為整體驗證測序數據的準確性，從所有顯著差異表達和非顯著差異表達的740個miRNA中隨機挑選2個非顯著差異表達miRNA（novel-m0031-3p和novel-m0034-5p）和2個顯著差異miRNA（miR-3793-x和ame-miR-6000a-5p）進行Stem-loopRT-PCR驗證，電泳結果顯示上述miRNA在AmCK和AmT中均可擴增出符合預期的目的條帶（圖6-A）。進一步對上述4個DEmiRNA進行qPCR檢測，結果顯示它們的表達水平的變化趨勢與sRNA-seq數據中的變化趨勢一致（圖6-B、6-C、6-D、6-E），證實了測序結果的可靠性。

A：上調miRNA的靶基因Target genes of up-regulated miRNAs；B：下調miRNA的靶基因Target genes of down-regulated miRNAs

A：上調miRNA的靶基因Target genes of up-regulated miRNAs；B：下調miRNA的靶基因Target genes of down-regulated miRNAs

3 討論

球囊菌是一種特異性侵染蜜蜂幼蟲腸道的真菌病原，能嚴重影響成年蜜蜂的數量和蜂群的生產力，造成較大的經濟損失[21]。miRNA作為一種轉錄后調控的關鍵因子，通過調節宿主相關免疫信號通路[22]，改變宿主基因的表達而影響病原復制[23]，而在宿主與病原間的互作中扮演著特殊角色。筆者所在課題組前期已在mRNA組學水平全面解析意蜂及中蜂幼蟲與球囊菌間的互作[24-26]。然而在ncRNA組學層面，至今還沒有二者互作的相關研究報道。本研究利用sRNA-seq技術對正常及球囊菌侵染的意蜂4日齡幼蟲腸道進行測序，預測出4個上調miRNA和6個下調miRNA，DEmiRNA在AmCK中的表達水平總體高于AmT（圖1-B），表明宿主在球囊菌侵染前期，通過下調相關miRNA的表達水平降低其對相應靶基因的抑制作用，從而響應對病原的應答。

3.1 球囊菌脅迫可影響與蜜蜂幼蟲腸道發育、代謝和免疫相關的miRNA的表達水平

球囊菌在侵染意蜂幼蟲的前期，孢子處于低水平的萌發狀態，幼蟲并不表現出明顯的白堊病癥狀，但宿主與病原之間存在復雜的互作[15]。Hippo信號通路參與腸道干細胞的增殖調控[27]，并通過與Wnt、Notch等信號通路共同作用影響的器官大小[28]。郭睿等[29-30]研究發現，意蜂工蜂中腸內的長鏈非編碼RNA和環狀RNA均可通過間接調控Hippo信號通路上的富集基因，影響蜜蜂腸道的發育過程。此外，該信號通路也被證明與蜜蜂的級型分化和卵巢發育狀態密切相關[31]。本研究中，對于意蜂幼蟲腸道的上調miRNA的靶基因，富集在Hippo信號通路的數量（62個），遠高于下調miRNA的靶基因富集在該通路上的數量（6個），表明球囊菌在侵染前期通過調控Hippo信號通路對宿主的腸道發育產生影響。郝向偉[32]研究發現，家蠶（）被大腸桿菌（）和蘇云金芽孢桿菌（）侵染后，其腸道的表達量顯著上調，作者推測Hedgehog信號通路可能參與腸干細胞的增殖和分化以及腸道損傷后的修復。本研究發現，對于富集在Hedgehog信號通路的靶基因，上調miRNA的靶基因數（43個）遠多于下調miRNA的靶基因數（1個），說明多數該通路上的基因受到抑制，表明Hedgehog信號通路對不同的病原微生物存在應答差異。本研究還發現，對于上調miRNA的靶基因，富集在晝夜節律（41個）、嘌呤代謝（48個）、氨基糖與核苷酸糖代謝（32個）和嘧啶代謝（5個）等代謝通路上的數量均高于下調miRNA的靶基因富集在上述通路上的數量。HUANG等[13]研究發現，東方蜜蜂微孢子蟲感染西方蜜蜂的前6 d內，宿主的17個DEmiRNA主要參與調控嘌呤代謝、嘧啶代謝和氧化磷酸化等9條代謝通路，作者推測這與通過改變宿主的代謝以促進其更快速的繁殖機制有關。本研究表明在球囊菌侵染前期，意蜂幼蟲腸道的部分生命活動、新陳代謝受到病原的抑制，推測球囊菌通過調節宿主的上述生物學過程為自身的孢子萌發創造有利條件。

A：上調miRNA的調控網絡Regulation networks of up-regulated miRNAs； B：下調miRNA的調控網絡Regulation networks of down-regulated miRNAs。綠色圓形代表靶基因，紅色三角代表miRNA Green circles indicate target genes, red triangles indicate miRNAs

A： Wnt信號通路相關DEmiRNA的調控網絡Regulation networks of DEmiRNAs related to Wnt signaling pathway； B：內吞作用相關DEmiRNA的調控網絡Regulation networks of DEmiRNAs related to endocytosis。綠色圓形代表靶基因，紅色三角代表miRNA Green circles indicate target genes, red triangles indicate miRNAs

圖6 DEmiRNA的Stem-loop RT-PCR（A）和qPCR驗證（B—E）

Wnt信號通路廣泛參與昆蟲的細胞增殖、組織分化、器官形成等過程[33-34]。該通路還可通過與TGF-beta、Hippo、Notch、MAPK、FoxO等信號通路的共同作用而參與蜂王卵巢激活和產卵過程[35]。本研究中，共有7個顯著DEmiRNA（miR-8440-y、miR-3793-x和miR-4968-y等）的96個靶基因注釋到Wnt信號通路，表明球囊菌侵染會影響宿主的Wnt信號通路基因的表達。此外，有研究表明Wnt信號通路參與了人類的先天免疫反應[36-37]，例如SMITH等[38]發現miR-34家族對Hela細胞內的Wnt/-catenin信號通路具有抑制作用，并表現出強烈的抗黃病毒作用；喬瑩[39]發現被隱核蟲（）感染的大黃魚（）體內miR-466-x、miR-1895-y和miR-8485-y等miRNA廣泛參與了對Wnt和Jak-STAT信號通路以及內吞作用和溶酶體等免疫通路的調控。本研究中，miR-4331-y、miR-4968-y、miR-8440-y、novel-m0023-5p和novel-m0025-3p均參與了對Wnt信號通路和內吞作用的調控，表明它們通過調控上述免疫通路對意蜂幼蟲腸道的免疫應答進行調節。

3.2 miR-4331-y作為潛在的分子靶標參與宿主與病原復雜的互作過程

miRNA在物種間具有高度的保守性、表達時序性和組織特異性[40]。郭睿等[41]對意蜂幼蟲腸道發育過程中miRNA的表達譜進行分析，發現miR-4331-y僅在4和5日齡幼蟲腸道內差異表達，并間接調控富集在Wnt、FoxO、TGF-beta、Notch和Jak-STAT等信號通路上的靶基因。ZHAO等[42]研究發現，豬腎細胞（PK-15）被傳染性腸胃炎病毒（，TGEV）侵染后，其miR-4331的表達量顯著上調，并通過靶向結合作用于p38 MAPK信號通路，從而調控TGEV誘導的PK-15線粒體損傷。此外，在患有H1N1流感病毒的病人以及感染低致病性H5N3流感病毒的雞的體內，均發現miR-204的表達顯著提高[43-44]。由于病原種類、侵染機制以及檢測時間點的差異，不同宿主的同源miRNA在響應病原侵染過程中，其表達量變化趨勢可能并不完全一致。例如ZHANG等[45]發現甲型H1N1病毒（SIV-H1N1/2009）可通過抑制新生豬氣管細胞（NPTr）中ssc-miR-204和ssc-miR-4331的表達促進自身在細胞內的復制；ssc-miR-204和ssc-miR-4331可抑制病毒的HA表面糖蛋白和NS1干擾素拮抗劑蛋白的表達，從而對病毒的復制進行負調控。本研究中，miR-4331-y的表達量在球囊菌脅迫后呈顯著下調，并通過對靶基因的調控同時參與對內吞作用、吞噬體、Jak-STAT和Wnt信號通路等免疫通路的間接調控，表明宿主在球囊菌侵染前期通過下調miR-4331-y的表達水平降低對靶基因的抑制作用，從而激活相關免疫通路，以響應球囊菌的脅迫；球囊菌也可能抑制宿主腸道內miR-4331-y的表達量，以利于孢子的萌發。未來可通過合成相應的miRNA mimic和miRNA inhibitor對miR-4331-y進行過表達和敲減，從而深入研究其在宿主的脅迫應答中的作用。

3.3 意蜂幼蟲腸道響應球囊菌脅迫前期過程中mRNA與miRNA組學比較

為探究不同蜂種的蜜蜂幼蟲在球囊菌侵染前期的應答及與病原間的互作，筆者課題組曾對正常及球囊菌脅迫的意蜂和中蜂4日齡幼蟲腸道進行轉錄組測序和分析[14-15]，發現意蜂幼蟲的差異表達基因涉及代謝進程、免疫系統進程和應激反應等19個功能分類，并有5個角質層蛋白基因和2個抗菌肽編碼基因表達水平下調，進而在mRNA組學水平解析了宿主的脅迫應答[15]。昆蟲腸道圍食膜是抵御經口攝入病原的重要物理屏障，角質層蛋白作為圍食膜的主要成分之一，在保護腸道健康過程中發揮重要作用[46]。本研究發現，上調miRNA中的miR-8440-y靶向結合2個角質層蛋白編碼基因，表明球囊菌可通過提高宿主的miR-8440-y的表達水平降低角質層蛋白編碼基因的水平，以利于球囊菌的萌發和侵染過程，這與此前的研究結果一致。

4 結論

在miRNA組學水平對意蜂幼蟲腸道在球囊菌侵染前期的脅迫應答及宿主-病原互作進行研究，結果表明顯著DEmiRNA可通過調控細胞生命活動、新陳代謝和部分細胞和體液免疫等生物學過程參與宿主的脅迫應答。其中miR-4331-y、miR-8440-y、miR-4968-y、novel-m0023-5p和novel-m0025-3p共同參與了宿主的Wnt信號通路和內吞作用的調控，且miR-8440-y結合靶基因的數量最多。篩選出的關鍵miRNA可用于后續的功能研究，有望為白堊病治療提供潛在的分子靶標。

[1] 羅術東, 王彪, 褚忠橋, 柳萌, 吳杰. 不同蜂為設施辣椒授粉的授粉效果比較. 環境昆蟲學報, 2015, 37(2): 381-386.

LUO S D, WANG B, CHU Z Q, LIU M, WU J. Comparison of the pollination effects for pepper between different bees in greenhouse., 2015, 37(2): 381-386. (in Chinese)

[2] 李江紅, 鄭志陽, 陳大福, 梁勤. 影響蜜蜂球囊菌侵染蜜蜂幼蟲的因素及侵染過程觀察. 昆蟲學報, 2012, 55(7): 790-797.

LI J H, ZHENG Z Y, CHEN D F, LIANG Q. Factors influencinginfection on honeybee larvae and observation on the infection process., 2012, 55(7): 790-797. (in Chinese)

[3] ASGARI S. MicroRNA functions in insects., 2013, 43(4): 388-397.

[4] CHEN J F, MANDEL E M, THOMSON J M, WU Q, CALLIS T E, HAMMOND S M, CONLON F L, WANG D Z. The role of microRNA-1 and microRNA-133 in skeletal muscle proliferation and differentiation., 2006, 38(2): 228-233.

[5] ALVAREZ-GARCIA I, MISKA E A. MicroRNA functions in animal development and human disease., 2005, 132(21): 4653-4662.

[6] GUO C J, PAN Q, LI D G, SUN H, LIU B W. MiR-15b and miR-16 are implicated in activation of the rat hepatic stellate cell: an essential role for apoptosis., 2009, 50(4): 766-778.

[7] SCARIA V, HARIHARAN M, MAITI S, PILLAI B, BRAHMACHARI S K. Host-virus interaction: a new role for microRNAs., 2006, 3: 68.

[8] LI S, SHEN L, SUN L, XU J, JIN P, CHEN L, MA F. Small RNA-Seq analysis reveals microRNA-regulation of the Imd pathway duringinfection in., 2017, 70: 80-87.

[9] QIAN P, JIANG T, WANG X, SONG F, CHEN C, SHEN X. Bmo-miR-275 down-regulates expression ofsericin gene 2., 2018, 13(1): e0190464.

[10] ZHANG G, HUSSAIN M, O’NEILL S L, ASGARI S.uses a host microRNA to regulate transcripts of a methyltransferase, contributing to dengue virus inhibition in., 2013, 110(25): 10276-10281.

[11] 李盛杰. microRNA在果蠅Toll信號免疫響應中的調控作用研究[D]. 南京: 南京師范大學, 2017.

Li S J. Regulation of microRNA on Toll signal immune response in[D]. Nanjing: Nanjing Normal University, 2017. (in Chinese)

[12] LOUREN?O A P, GUIDUGLILAZZARINI K R, FREITAS F C, BITONDI M M, Sim?es Z L. Bacterial infection activates the immune system response and dysregulates microRNA expression in honey bees., 2013, 43(5): 474-482.

[13] HUANG Q, CHEN Y P, Rui W W, SCHWARZ R S, EVANS J D. Honeybee microRNAs respond to infection by the microsporidian parasite., 2015, 5: Article number 17494.

[14] 陳大福, 郭睿, 熊翠玲, 梁勤, 鄭燕珍, 徐細建, 張曌楠, 黃枳腱, 張璐, 王鴻權, 解彥玲, 童新宇. 中華蜜蜂幼蟲腸道響應球囊菌早期脅迫的轉錄組學. 中國農業科學, 2017, 50(13): 2614-2623.

CHEN D F, GUO R, XIONG C L, LIANG Q, ZHENG Y Z, XU X J, ZHANG Z N, HUANG Z J, ZHANG L, WANG H Q, XIE Y L, TONG X Y. Transcriptome oflarval gut under the stress of., 2017, 50(13): 2614-2623. (in Chinese)

[15] 郭睿, 熊翠玲, 鄭燕珍, 張璐, 童新宇, 梁勤, 陳大福. 意大利蜜蜂幼蟲腸道響應球囊菌早期脅迫的轉錄組學分析. 應用昆蟲學報, 2017, 54(4): 553-560.

GUO R, XIONG C L, ZHENG Y Z, ZHANG L, TONG X Y, LIANG Q, CHEN D F. Transcriptome analysis oflarval gut during the early stage of stress induced by., 2017, 54(4): 553-560. (in Chinese)

[16] HUSSAIN M, ASGARI S. MicroRNAs as mediators of insect host-pathogen interactions and immunity., 2014, 70: 151-158.

[17] 郭睿, 王海朋, 陳華枝, 熊翠玲, 鄭燕珍, 付中民, 趙紅霞, 陳大福. 蜜蜂球囊菌的microRNA鑒定及其調控網絡分析. 微生物學報, 2018, 58(6): 1077-1089.

GUO R, WANG H P, CHEN H Z, XIONG C L, ZHENG Y Z, FU Z M, ZHAO H X, CHEN D F. Identification ofmicroRNAs and investigation of their regulation networks., 2018, 58(6): 1077-1089. (in Chinese)

[18] FRIEDL?NDER M R, MACKOWIAK S D, LI N, CHEN W, RAJEWSKY N. MiRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades., 2012, 40(1): 37-52.

[19] ALLEN E, XIE Z, GUSTAFSON A M, Carrington J C. MicroRNA-directed phasing during trans-acting siRNA biogenesis in plants., 2005, 121(2): 207-221.

[20] CHEN C, RIDZON D A, BROOMER A J, Lee D H, Nguyen J T, Barbisin M, Xu N L, Mahuvakar V R, Andersen M R, Lao K Q, Livak K J, Guegler K J. Real-time quantification of microRNAs by stem-loop RT-PCR., 2005, 33(20): e179.

[21] 趙紅霞, 梁勤, 羅岳雄, 李江紅, 張學鋒, 曾鑫年. 蜜蜂白堊病的研究進展. 環境昆蟲學報, 2014, 36(2): 233-239.

ZHAO H X, LIANG Q, LUO Y X, LI J H, ZHANG X F, ZENG X N. Chalkbrood disease in honeybee., 2014, 36(2): 233-239. (in Chinese)

[22] ZHAO L, ZHU J, ZHOU H, ZHAO Z, ZOU Z, LIU X, LIN X, ZHANG X, DENG X, WANG R, CHEN H, JIN M. Identification of cellular microRNA-136 as a dual regulator of RIG-I-mediated innate immunity that antagonizes H5N1 IAV replication in A549 cells., 2015, 5: Article number14991.

[23] HU Y, JIANG L, LAI W, QIN Y, ZHANG T, WANG S, YE X. MicroRNA-33a disturbs influenza A virus replication by targeting ARCN1 and inhibiting viral ribonucleoprotein activity., 2016, 97(1): 27-38.

[24] CHEN D, GUO R, XU X, XIONG C, LIANG Q, ZHENG Y, LUO Q, ZHANG Z, HUANG Z, KUMAR D, XI W, ZOU X, LIU M. Uncovering the immune responses oflarval gut toinfection utilizing transcriptome sequencing., 2017, 621: 40-50.

[25] 陳大福, 郭睿, 熊翠玲, 梁勤, 鄭燕珍, 徐細建, 黃枳腱, 張曌楠, 張璐, 李汶東, 童新宇, 席偉軍. 脅迫意大利蜜蜂幼蟲腸道的球囊菌的轉錄組分析. 昆蟲學報, 2017, 60(4): 401-411.

CHEN D F, GUO R, XIONG C L, LIANG Q, ZHENG Y Z, XU X J, HUANG Z J, ZHANG Z N, ZHANG L, LI W D, TONG X Y, XI W J. Transcriptomic analysis ofstressing larval gut of(Hyemenoptera: Apidae) ., 2017, 60(4): 401-411. (in Chinese)

[26] 郭睿, 陳大福, 黃枳腱, 梁勤, 熊翠玲, 徐細建, 鄭燕珍, 張曌楠, 解彥玲, 童新宇, 候志賢, 江亮亮, 刀晨. 球囊菌脅迫中華蜜蜂幼蟲腸道過程中病原的轉錄組學研究. 微生物學報, 2017, 57(12): 1865-1878.

GUO R, CHEN D F, HUANG Z J, LIANG Q, XIONG C L, XU X J, ZHENG Y Z, ZHANG Z N, XIE Y L, TONG X Y, HOU Z X, JIANG L L, DAO C. Transcriptome analysis ofstressing larval gut of., 2017, 57(12): 1865-1878. (in Chinese)

[27] AYYAZ A, LI H, JASPER H. Hemocytes control stem cell activity in theintestine., 2015, 17(6): 736-748.

[28] BARRY E R, CAMARGO F D. The Hippo superhighway: signaling crossroads converging on the Hippo/Yap pathway in stem cells and development., 2013, 25(2): 247-253.

[29] 郭睿, 耿四海, 熊翠玲, 鄭燕珍, 付中民, 王海朋, 杜宇, 童新宇, 趙紅霞, 陳大福. 意大利蜜蜂工蜂中腸發育過程中長鏈非編碼RNA的差異表達分析. 中國農業科學, 2018, 51(18): 3600-3613.

GUO R, GENG S H, XIONG C L, ZHENG Y Z, FU Z M, WANG H P, DU Y, TONG X Y, ZHAO H X, CHEN D F. Differential expression analysis of long non-coding RNAs during the developmental process ofworker’s midgut., 2018, 51(18): 3600-3613. (in Chinese)

[30] 郭睿, 陳華枝, 熊翠玲, 鄭燕珍, 付中民, 徐國鈞, 杜宇, 王海朋, 耿四海, 周丁丁, 劉思亞, 陳大福. 意大利蜜蜂工蜂中腸發育過程中的差異表達環狀RNA及其調控網絡分析. 中國農業科學, 2018, 51(23): 4575-4590.

GUO R, CHEN H Z, XIONG C L, ZHENG Y Z, FU Z M, XU G J, DU Y, WANG H P, GENG S H, ZHOU D D, LIU S Y, CHEN D F. Analysis of differentially expressed circular RNAs and their regulation networks during the developmental process ofworker’s midgut., 2018, 51(23): 4575-4590. (in Chinese)

[31] ASHBY R, FORêT S, SEARLE I, MALESZKA R. MicroRNAs in honeybee caste determination., 2016, 6: Article number18794.

[32] 郝向偉. 家蠶Hedgehog信號通路相關基因的克隆、鑒定及其功能分析[D]. 重慶: 西南大學, 2013.

HAO X W. Cloning, characterization and functional analysis of genes in the hedgehog signaling pathway in the intestine of silkworm,[D]. Chongqing: Southwest University, 2013. (in Chinese)

[33] CHA Y H, KIM N H, PARK C, LEE I, KIM H S, YOOK J I. MiRNA-34 intrinsically links p53 tumor suppressor and Wnt signaling., 2012, 11(7): 1273-1281.

[34] KIM N H, KIM H S, KIM N G, LEE I, CHOI H S, LI X Y, KANG S E, CHA S Y, RYU J K, NA J M, PARK C, KIM K, LEE S, GUMBINER B M, YOOK J I, WEISS S J. P53 and microRNA-34 are suppressors of canonical Wnt signaling., 2011, 4(197): ra71.

[35] 陳曉. 蜜蜂卵巢激活和產卵過程差異表達的編碼RNA與非編碼RNA的篩選和鑒定[D]. 北京: 中國農業科學院, 2017.

CHEN X. Identification of differentially expressed coding and noncoding RNAs during ovary activation and oviposition in honeybees[D]. Beijing: Chinese Academy of Agricultural Sciences, 2017. (in Chinese)

[36] BARIL M, ES-SAAD S, CHATEL-CHAIX L, FINK K, PHAM T, RAYMOOD V A, AUDETTE K, GUENIER A S, DUCHAINE J, SERVANT M, BILODEAU M, COHEN E, GRANDVAUX N, LAMARRE D. Genome-wide RNAi screen reveals a new role of a Wnt/CTNNB1 signaling pathway as negative regulator of virus-induced innate immune responses., 2013, 9(6): e1003416.

[37] HACK K, REILLY L, PROBY C, FLEMING C, LEIGH I, FOERSTER J. Wnt5a inhibits the CpG oligodeoxynucleotide- triggered activation of human plasmacytoid dendritic cells., 2012, 37(5): 557-561.

[38] SMITH J L, JENG S, MCWEENEY S K, HIRSCH A J. A microRNA screen identifies the Wnt signaling pathway as a regulator of the interferon response during flavivirus infection., 2017, 91(8): e02388-16.

[39] 喬瑩. 刺激隱核蟲感染大黃魚的miRNA和mRNA組學測序及關聯分析[D]. 廈門: 廈門大學, 2016.

QIAO Y. Sequencing and correlation analysis of miRNAome and transcriptome ofinfected by[D]. Xiamen: Xiamen University, 2016. (in Chinese)

[40] HE L, HANNON G J. MicroRNAs: small RNAs with a big role in gene regulation., 2004, 5(7): 522-531.

[41] 郭睿, 杜宇, 熊翠玲, 鄭燕珍, 付中民, 徐國鈞, 王海朋, 陳華枝, 耿四海, 周丁丁, 石彩云, 趙紅霞, 陳大福. 意大利蜜蜂幼蟲腸道發育過程中的差異表達microRNA及其調控網絡. 中國農業科學, 2018, 51(21): 4197-4209.

GUO R, DU Y, XIONG C L, ZHENG Y Z, FU Z M, XU G J, WANG H P, CHEN H Z, GENG S H, ZHOU D D, SHI C Y, ZHAO H X, CHEN D F. Differentially expressed microRNA and their regulation networks during the developmental process oflarval gut., 2018, 51(21): 4197-4209. (in Chinese)

[42] ZHAO X, BAI X, GUAN L, LI J, SONG X, MA X, GUO J, ZHANG Z, DU Q, HUANG Y, TONG D. MicroRNA-4331 promotes Transmissible Gastroenteritis Virus (TGEV)-induced mitochondrial damage via targeting RB1, upregulating interleukin-1 receptor accessory protein (IL1RAP), and activating p38 MAPK Pathway., 2018, 17(2): 190-204.

[43] WANG Y, BRAHMAKSHATRIYA V, LUPIANI B, REDDY S M, SOIBAM B, BENHAM A L, GUNARATNE P, LIU H C, TRAKOOLJUL N, ING N, OKIMOTO R, ZHOU H. Integrated analysis of microRNA expression and mRNA transcriptome in lungs of avian influenza virus infected broilers., 2012, 13: 278.

[44] TAMBYAH P A, SEPRAMANIAM S, MOHAMED ALI J, CHAI S C, SWAINATHAN P, ARMUGAM A, JEYASEELAN K. MicroRNAs in circulation are altered in response to influenza a virus infection in humans., 2013, 8(10): e76811.

[45] ZHANG S, WANG R, SU H, WANG B, SIZHU S, LEI Z, JIN M, CHEN H, CAO J, ZHOU H.miR-204 and miR-4331 negatively regulate swine H1N1/2009 influenza a virus replication by targeting viral HA and NS, respectively., 2017, 18(4): E749.

[46] WANG P, GRANADOS R R. Molecular structure of the peritrophic membrane (PM): identification of potential PM target sites for insect control., 2001, 47(2): 110-118.

Differentially expressed microRNAs and their regulation networks inlarval gut during the early stage ofinfection

GUO Rui, DU Yu, TONG XinYu, XIONG CuiLing, ZHENG YanZhen, XU GuoJun, WANG HaiPeng, GENG SiHai, ZHOU DingDing, GUO YiLong, WU SuZhen, CHEN DaFu

(College of Bee Science, Fujian Agriculture and Forestry University, Fuzhou 350002)

【Objective】MicroRNA (miRNA) is a kind of key gene expression regulator, which can affect the interactions between host and pathogen.is a lethal fungal pathogen that specifically infects honeybee larvae. The objective of this study is to analyze the differentially expressed miRNAs (DEmiRNAs) and their target genes in thelarval gut during the early infection stage of, reveal DEmiRNA’ roles in the stress responses of host at the miRNA omics level, and to screen the key miRNAs related to host response by constructing regulation networks of significant DEmiRNAs. 【Method】Normal and-infected 4-day-old larval gut of(AmCK and AmT) were deep-sequenced using small RNA-seq (sRNA-seq) technology, followed by quality-control of raw data and then mapping of the filtered data with the reference genome of. The mapped tags were compared to the miRBase database to identify the expression of known miRNAs. The expression of miRNAs in each sample was normalized by TPM (tags per million) algorithm and significant DEmiRNAs were gained according to the standard |log2fold change|≥1 and≤0.05. Target genes of significant DEmiRNAs were predicted utilizing TargetFinder, and then annotated to the GO and KEGG databases. Cytoscape was used to visualize the regulation networks between significant DEmiRNAs and target mRNAs. Finally, Stem-loop RT-PCR and qPCR were conducted to verify the reliability of the sequencing data.【Result】sRNA-seq of AmCK and AmT produced 13 553 302 and 10 777 534 raw reads, and after strict filtration, 13 186 921 and 10 480 913 clean reads were obtained, respectively. The Pearson correlation coefficients among different biological replicates in each sample were above 0.9822 and 0.9508. There were 10 significant DEmiRNAs including 4 up-regulated miRNAs and 6 down-regulated miRNAs, and the overall expression level of DEmiRNAs in AmT was lower than that in AmCK. In total, 10 significant DEmiRNAs could link 3 788 target genes. The 1 240 target genes of up-regulated miRNAs could be annotated to 39 GO terms, and the mostly enriched terms were binding, cellular processes, metabolic processes, and response to stimulus. The 749 target genes of down-regulated miRNAs could be annotated to 34 GO terms, and the mostly enriched terms were cellular processes, binding, metabolic processes, and response to stimulus. The result of KEGG database annotation suggested that the target genes of up- and down-regulated miRNAs were respectively annotated in 95 and 66 pathways, the most abundant pathways were Wnt signaling pathway, Hippo signaling pathway, phototransduction and endocytosis, phosphatidylinositol signaling system, as well as purine metabolism. For up- and down-regulated miRNAs, there were 31 and 52 target genes could be annotated to endocytosis, 15 and 7 target genes could be annotated to ubiquitin-mediated proteolysis, 11 and 5 target genes could be annotated to Jak-STAT signaling pathway, 1 and 3 target genes could be annotated to the MAPK signaling pathway, respectively. Complex regulation networks existed between significant DEmiRNAs and their target mRNAs, among them 7 significant DEmiRNAs targeted 96 mRNAs associated with Wnt signaling pathway, and 8 significant DEmiRNAs targeted 55 mRNAs involved in endocytosis. Finally, the results of Stem-loop RT-PCR and qPCR verified the reliability of the sequencing data.【Conclusion】larval gut’s DEmiRNAs and their target genes during the early infection stage ofwere predicted and analyzed. DEmiRNA-mRNA regulation networks in the host were constructed and investigated. The results provide the expression profile and differential expression information of host miRNAs, and reveal that these DEmiRNAs likely participate in the stress responses of host via regulating biological processes such as cellular activity, metabolism, and immune defense. miR-4331-y, miR-4968-y, miR-8440-y, novel-m0023-5p and novel-m0025-3p jointly regulate Wnt signaling pathway and endocytosis of host and can be used as potential molecular targets for chalkbrood control.

; larval gut; development; differentially expressed microRNA; regulation network

10.3864/j.issn.0578-1752.2019.01.015

2018-08-02；

2018-10-01

國家自然科學基金（31702190）、國家現代農業產業技術體系建設專項資金（CARS-44-KXJ7）、福建省科技計劃（2018J05042）、福建省教育廳中青年教師教育科研項目（JAT170158）、福建農林大學科技創新專項基金（CXZX2017343）、福建省大學生創新創業訓練計劃（201710389058，201810389082）

郭睿，E-mail：ruiguo@fafu.edu.cn。杜宇，E-mail：m18505700830@163.com。郭睿和杜宇為同等貢獻作者。通信作者陳大福，E-mail：dfchen826@fafu.edu.cn

（責任編輯 岳梅）