杜宇,周丁丁,萬潔琦,盧家軒,范小雪,范元嬋,陳恒,熊翠玲,鄭燕珍,付中民,徐國鈞,陳大福,郭睿
意大利蜜蜂工蜂中腸發(fā)育過程中的差異基因表達譜及調(diào)控網(wǎng)絡(luò)
杜宇,周丁丁,萬潔琦,盧家軒,范小雪,范元嬋,陳恒,熊翠玲,鄭燕珍,付中民,徐國鈞,陳大福,郭睿
(福建農(nóng)林大學(xué)動物科學(xué)學(xué)院(蜂學(xué)學(xué)院),福州 350002)
【】前期已對意大利蜜蜂(,簡稱意蜂)7日齡工蜂中腸(Am7)、10日齡工蜂中腸(Am10)進行全轉(zhuǎn)錄組測序,本研究基于高質(zhì)量的組學(xué)數(shù)據(jù)探究中腸發(fā)育過程中的差異基因表達譜及其調(diào)控網(wǎng)絡(luò),以期解析意蜂工蜂中腸發(fā)育的分子機理。根據(jù)FPKM(fragments per kilobase of transcript per million mapped reads)算法計算基因表達量,并以|log2fold change|≥1且≤0.05作為標準篩選得到差異表達基因(differentially expressed gene,DEG)。利用TargetFinder軟件預(yù)測ame-miR-6001-3p的靶mRNA。利用相關(guān)生物信息學(xué)軟件,對全部DEG進行GO和KEGG數(shù)據(jù)庫注釋。篩選出與AMPK、P13K-Akt、Wnt、cAMP、FoxO、Hippo、mTOR、Jak-STAT、Toll-like受體、TGF-beta、Notch、MAPK和NF-κB 13條信號通路存在富集關(guān)系的DEG,以及與ame-miR-6001-3p存在靶向結(jié)合關(guān)系的DEG,通過Cytoscape軟件構(gòu)建調(diào)控網(wǎng)絡(luò)將其富集關(guān)系與調(diào)控關(guān)系可視化。利用莖環(huán)反轉(zhuǎn)錄PCR(Stem loop RT-PCR)和實時熒光定量PCR(RT-qPCR)驗證ame-miR-6001-3p以及DEG在Am7和Am10中的差異表達情況。Am7 vs Am10比較組中共有1 038個DEG,包括515個上調(diào)基因和523個下調(diào)基因。這些DEG涉及細胞進程、代謝進程和催化活性等功能條目,并顯著富集在氧化磷酸化、氨基糖與核苷酸糖代謝、脂肪酸代謝和嘌呤代謝等能量和物質(zhì)代謝通路,表明工蜂中腸內(nèi)存在旺盛細胞生命與新陳代謝活動。表達量聚類分析發(fā)現(xiàn)分別有20、18、15和14個DEG富集在AMPK信號通路、P13K-Akt信號通路、內(nèi)吞作用和Hippo信號通路。57個DEG與P13K-Akt、Wnt、Jak-STAT等上述13條與生長發(fā)育和免疫防御相關(guān)的信號通路存在富集關(guān)系,且1個DEG可與多條信號通路存在富集關(guān)系。調(diào)控網(wǎng)絡(luò)分析結(jié)果顯示,分別有54個上調(diào)基因和44個下調(diào)基因可被ame-miR-6001-3p靶向結(jié)合;上調(diào)基因富集在磷酸肌醇代謝、胰島素信號通路、Hippo信號通路和谷胱甘肽代謝等43條代謝通路,而下調(diào)基因富集在Hippo信號通路、新陳代謝途徑、谷胱甘肽代謝和花生四烯酸代謝等20條代謝通路。RT-qPCR結(jié)果顯示隨機挑選的6個DEG的表達量變化趨勢與測序數(shù)據(jù)一致,證實了本研究中基因差異表達真實可靠。此外ame-miR-6001-3p在Am7和Am10內(nèi)均真實表達,并在Am10中的相對表達量顯著較低。對意蜂工蜂中腸發(fā)育過程的DEG表達譜和DEG與ame-miR-6001-3p之間的調(diào)控關(guān)系,以及DEG的潛在作用進行深入分析和探討,發(fā)現(xiàn)DEG可參與TGF-beta、Wnt、Hippo、Notch、PI3K-Akt、mTOR、AMPK和NF-κB等各類信號通路進而影響中腸的生長發(fā)育和免疫防御,DEG可通過與顯著下調(diào)表達的ame-miR-6001-3p形成復(fù)雜的調(diào)控網(wǎng)絡(luò)參與意蜂工蜂中腸發(fā)育過程中胰島素信號通路等多條代謝途徑。
意大利蜜蜂;中腸;發(fā)育;差異表達基因;競爭性內(nèi)源RNA;調(diào)控網(wǎng)絡(luò)
【研究意義】意大利蜜蜂(,簡稱意蜂)是經(jīng)濟價值和生態(tài)價值優(yōu)越的社會性模式昆蟲,廣泛用于國內(nèi)外的養(yǎng)蜂生產(chǎn)。蜜蜂中腸作為營養(yǎng)消化吸收、能量物質(zhì)代謝以及免疫應(yīng)答防御的主要場所,其發(fā)育機理至今仍不明了。因此,通過分析差異表達基因(differentially expressed gene,DEG)在腸道發(fā)育過程中的表達譜、DEG的潛在作用及調(diào)控網(wǎng)絡(luò),可為明確意蜂中腸發(fā)育的分子機理提供理論依據(jù)。【前人研究進展】昆蟲腸道的內(nèi)環(huán)境穩(wěn)態(tài)可影響個體的發(fā)育進程。COON等研究發(fā)現(xiàn),埃及伊蚊()和岡比亞按蚊()腸道內(nèi)的部分菌群是維持個體正常發(fā)育的關(guān)鍵因子[1]。抗生素處理可導(dǎo)致斯氏按蚊()幼蟲發(fā)育緩慢[2],而植物乳桿菌()可以通過影響激素信號受體促進黑腹果蠅()的個體增長[3]。目前有關(guān)蜜蜂腸道的研究主要涉及腸道菌群的結(jié)構(gòu)和功能預(yù)測、外源營養(yǎng)素與腸道發(fā)育的關(guān)系等方面。ZHENG等[4]通過比較無菌西方蜜蜂()與正常西方蜜蜂的生長速率和腸道質(zhì)量,發(fā)現(xiàn)腸道菌群會促進西方蜜蜂腸道的發(fā)育并增加個體質(zhì)量;馮倩倩[5]利用添加維生素A的代用花粉飼喂意蜂,通過檢測中腸蛋白質(zhì)濃度、堿性磷酸酶活性和圍食膜厚度等指標發(fā)現(xiàn)適量的維生素A有益于意蜂腸道的發(fā)育。近年來以RNA-seq為代表的第二代轉(zhuǎn)錄組測序技術(shù)迅猛發(fā)展,為蜜蜂中腸發(fā)育機理的深入研究提供了有力工具。筆者團隊前期利用二代測序技術(shù)和趨勢分析方法解析了意蜂幼蟲腸道發(fā)育過程中的差異基因表達譜[6],并探究了微小RNA(microRNA,miRNA)在幼蟲腸道發(fā)育過程中的調(diào)控作用[7]。前人研究發(fā)現(xiàn)ame-miR-6001-3p作為關(guān)鍵的基因表達調(diào)控因子,可參與調(diào)控蜜蜂蛻皮激素的分泌[8]和級型分化過程[9]。結(jié)合鏈特異性建庫的長鏈非編碼RNA(long non-coding RNA,lncRNA)測序技術(shù)、small RNA-seq(sRNA-seq)技術(shù),筆者團隊前期已對意蜂7和10日齡工蜂中腸(Am7和Am10)進行了深度測序,獲得了高質(zhì)量的全轉(zhuǎn)錄組數(shù)據(jù),并通過深入的生物信息學(xué)分析發(fā)現(xiàn)意蜂工蜂中腸內(nèi)的lncRNA和環(huán)狀RNA(circular RNA,circRNA)可作為競爭性內(nèi)源RNA(ceRNA)吸附結(jié)合ame-miR-6001-3p,從而間接調(diào)控mRNA的表達,影響中腸的細胞分化、新陳代謝、信號傳導(dǎo)和免疫調(diào)控[10-11]。【本研究切入點】然而,前期研究均是在單一ncRNA組學(xué)層面探究意蜂工蜂中腸的發(fā)育機理,缺乏mRNA組學(xué)層面的分析和闡釋,也沒有綜合多個層面組學(xué)數(shù)據(jù)對編碼RNA與ncRNA的相互調(diào)控關(guān)系進行深入探討。【擬解決的關(guān)鍵問題】基于已獲得的高質(zhì)量全轉(zhuǎn)錄組數(shù)據(jù),通過生物信息學(xué)方法對意蜂工蜂中腸發(fā)育過程中的DEG及其潛在作用進行深入分析,并結(jié)合前期在lncRNA組學(xué)和circRNA組學(xué)層面的分析篩選出的ame-miR-6001-3p,構(gòu)建和分析mRNA與miRNA之間的調(diào)控網(wǎng)絡(luò),以期更深入地解析意蜂工蜂中腸的發(fā)育機理。
試驗于2017—2019年在福建農(nóng)林大學(xué)動物科學(xué)學(xué)院(蜂學(xué)學(xué)院)蜜蜂保護實驗室進行。
選取福建農(nóng)林大學(xué)教學(xué)蜂場內(nèi)的意蜂蜂群3群,提取9張老熟封蓋子脾至實驗室的培養(yǎng)箱內(nèi),箱內(nèi)溫度(34±0.5)℃,相對濕度70%,收集剛羽化出房工蜂100只。
中腸樣品的制備步驟簡述如下:每10只剛羽化(記為0日齡)的工蜂放入1個四周打孔通風(fēng)的干凈塑料盒內(nèi),盒子上方裝有50%(w/v)無菌糖水的飼喂器,置于培養(yǎng)箱中培養(yǎng)(條件同1.1)。在超凈工作臺內(nèi)分別剖取飼養(yǎng)至7日齡工蜂中腸(Am7)和10日齡工蜂中腸(Am10),每3只中腸置于1個RNA-Free的EP管中,液氮速凍后迅速放入-80℃超低溫冰箱保存。Am7和Am10均設(shè)置3個生物學(xué)重復(fù)。鑒于東方蜜蜂微孢子蟲()完成一個生活史需要6 d[12]以及10 d仍處于孢子不斷增殖時期[13],同時為滿足探究意蜂工蜂中腸發(fā)育機理及意蜂對主要真菌病原的脅迫應(yīng)答,本研究選擇上述兩個取樣時間點。
上述6個中腸樣品的cDNA文庫構(gòu)建、lncRNA- seq實驗技術(shù)流程參數(shù)參照郭睿等[10]方法,委托廣州基迪奧生物技術(shù)有限公司進行,測序平臺分別為Illumina HiseqTM4000。相關(guān)原始數(shù)據(jù)已上傳至NCBI SRA數(shù)據(jù)庫,BioProject號:PRJNA406998。
利用Perl腳本剔除測序數(shù)據(jù)原始讀段中未知核苷酸含量>5%以及含量>50%的質(zhì)量低的讀段(reads),獲得有效讀段(clean reads)。過濾可比對上核糖體RNA數(shù)據(jù)庫的reads,進而通過TopHat2將unmapped reads比對到西方蜜蜂參考基因組(Amel_4.5)(https://www.ncbi.nlm.nih.gov/assembly/GCF_000002195.4)。校準參數(shù)如下:(1)最大讀取不匹配為2;(2)配對讀取目標區(qū)域為50 bp;(3)配對讀取目標區(qū)域誤差為±80 bp。利用Cufflinks將對比到轉(zhuǎn)錄本上的reads進行組裝。根據(jù)FPKM(fragments per kilobase of transcript per million mapped reads)算法計算基因表達量,并利用R軟件計算Pearson相關(guān)性系數(shù)來評估同組樣品的生物學(xué)重復(fù)性。
將|log2fold change|≥1且值≤0.05的基因定義為DEG。利用Blast軟件將DEG注釋到GO數(shù)據(jù)庫和KEGG數(shù)據(jù)庫,繪制生物學(xué)進程、細胞組分和分子功能三大類的GO富集分類柱狀統(tǒng)計圖,并通過在線分析軟件OmicsShare(http://www.omicshare.com/ tools/Home/Index/index.html)繪制表達量聚類熱圖。ame-miR-6001-3p是基于意蜂工蜂中腸的全轉(zhuǎn)錄組數(shù)據(jù)及前期分析結(jié)果篩選得到的候選靶標,該靶標可靶向結(jié)合多個意蜂lncRNA和circRNA,且位于調(diào)控網(wǎng)絡(luò)的中心[10-11]。利用TargetFinder軟件對ame-miR- 6001-3p的靶mRNA進行預(yù)測,并與本研究的DEG進行比較。最后,根據(jù)DEG與ame-miR-6001-3p以及KEGG代謝通路(pathway)的關(guān)系構(gòu)建調(diào)控網(wǎng)絡(luò),并利用Cytoscape軟件進行網(wǎng)絡(luò)可視化[14]。
為驗證測序數(shù)據(jù)的可靠性,隨機選取3個上調(diào)基因(XM_016912284.1、TCONS_00023605和XM_ 016912389.1)和3個下調(diào)基因(TCONS_00001615、XM_001123053.4和XM_016915022.1)進行RT-qPCR驗證。首先利用RNA抽提試劑盒(TaKaRa公司,中國)提取Am7和Am10樣品的總RNA,反轉(zhuǎn)錄得到cDNA作為模板進行RT-qPCR。按照SYBR Green Dye試劑盒(Vazyme公司,中國)操作說明書進行。反應(yīng)體系20 μL包含:SYBR Green Dye 10 μL,正、反向引物(10.0 μmol·L-1)各1 μL,cDNA模板1 μL,DEPC水7 μL。RT-qPCR反應(yīng)在ABI 7500熒光定量PCR儀(ABI公司,美國)上進行,反應(yīng)條件如下:95℃預(yù)變性1 min;95℃變性15 s,60℃延伸30 s,共40個循環(huán);最后72℃延伸45 s。本試驗進行3次技術(shù)性重復(fù)。以作為內(nèi)參基因,所選基因的相對表達量采用2-??Ct法計算。
利用總RNA和ame-miR-6001-3p的特異性Stem-loop引物進行反轉(zhuǎn)錄得到cDNA模板,結(jié)合上下游引物,進行常規(guī)PCR反應(yīng),PCR產(chǎn)物經(jīng)1.5%瓊脂糖凝膠電泳檢測。以snRNA U6為內(nèi)參,參照筆者團隊前期已建立的RT-qPCR流程[7,15]對ame-miR- 6001-3p在Am7和Am10中的相對表達量進行檢測。
通過GraphPad Prism軟件繪制上述DEG及ame- miR-6001-3p的相對變化量的柱狀圖,利用SPSS軟件計算組間顯著性。本研究所用相關(guān)引物信息詳見表1。
前期研究中已對Am7和Am10的lncRNA-seq數(shù)據(jù)進行了嚴格質(zhì)控,結(jié)果表明意蜂工蜂中腸樣品的生物學(xué)重復(fù)性且全轉(zhuǎn)錄組數(shù)據(jù)質(zhì)量良好[10],可滿足本研究的分析需要。差異分析結(jié)果顯示Am7 vs Am10比較組包含1 038個DEG,其中上調(diào)和下調(diào)基因的數(shù)量分別為515和523個(圖1)。

表1 本研究所用引物

圖1 意蜂工蜂中腸發(fā)育過程中DEG的火山圖
GO分類結(jié)果顯示,DEG共注釋到46個功能條目,其中上調(diào)和下調(diào)基因可分別注釋到36和44個功能條目。進一步分析發(fā)現(xiàn),細胞組分包含基因數(shù)比較多的分別是細胞組件(158個DEG)、細胞(158個DEG)、細胞膜(118個DEG)、細胞膜組件(114個DEG)和細胞器(111個DEG);分子功能包含基因數(shù)比較多的分別是結(jié)合(350個DEG)、催化活性(258個DEG)、轉(zhuǎn)運子活性(59個DEG)、分子傳感器活性(29個DEG)和分子功能調(diào)節(jié)器(28個DEG);生物學(xué)進程包含基因數(shù)比較多的分別是細胞進程(285個DEG)、代謝進程(273個DEG)、單組織進程(240個DEG)、生物學(xué)調(diào)控(112個DEG)和定位(95個DEG)(圖2)。

圖2 意蜂工蜂中腸發(fā)育過程中DEG的GO分類
KEGG代謝通路富集分析結(jié)果顯示,DEG共富集在260條代謝通路。其中,AMPK信號通路(20個DEG)、PI3K-Akt信號通路(18個DEG)、內(nèi)吞作用(15個DEG)、Hippo信號通路(14個DEG)、泛素介導(dǎo)的蛋白質(zhì)水解(14個DEG)、胰島素信號通路(12個DEG)、磷脂酰肌醇信號轉(zhuǎn)導(dǎo)系統(tǒng)(12個DEG)、Wnt信號通路(12個DEG)、磷脂酶D信號通路(11個DEG)、鞘脂信號通路(11個DEG)等通路富集的DEG數(shù)量較多;此外,部分DEG富集于生長發(fā)育、新陳代謝、免疫防御等相關(guān)通路,包括溶酶體(11個DEG)、FoxO信號通路(9個DEG)、吞噬體(9個DEG)、cAMP信號通路(9個DEG)、磷酸肌醇代謝(8個DEG)、Hedgehog信號通路(8個DEG)、碳代謝(8個DEG)、mTOR信號通路(7個DEG)、氨基糖與核苷糖代謝(6個DEG)、脂肪酸代謝(6個DEG)、MAPK信號通路(6個DEG)、嘌呤代謝(5個DEG)、Toll-like受體信號通路(4個DEG)。
進一步對富集在AMPK信號通路、PI3K-Akt信號通路、Hippo信號通路和內(nèi)吞作用上的DEG進行表達量聚類分析,結(jié)果顯示各通路富集的上調(diào)基因數(shù)分別為12、8、13和9個,均高于下調(diào)基因數(shù),表明上述4條代謝通路被不同程度的激活(圖3)。
57個DEG在13條信號通路(AMPK、P13K-Akt、Wnt、cAMP、FoxO、Hippo、mTOR、Jak-STAT、Toll-like受體、TGF-beta、Notch、MAPK和NF-B)的富集網(wǎng)絡(luò)關(guān)系,利用Cytoscape軟件進行可視化,結(jié)果顯示一個基因可同時富集在多條信號通路(圖4),例如磷脂酰肌醇3-激酶調(diào)節(jié)亞基-類編碼基因(XM_016917733.1和XM_392119.6)富集在7條信號通路(AMPK、PI3-Akt、cAMP、FoxO、mTOR、Jak-STAT和Toll-like受體);Ras相關(guān)蛋白rac1亞型編碼基因(TCONS_00021690和TCONS_00038549)富集在5條信號通路(P13K-Akt、cAMP(NF-B、AMPK)、Toll-like受體、Wnt、MAPK)(圖4);Smad(mothers against decapentaplegic homolog)3亞型編碼基因(TCONS_00001946)與4條信號通路(Wnt、FoxO、Hippo和TGF-beta)存在富集關(guān)系;類胰島素樣受體編碼基因(XM_016918001.1)、mTOR的調(diào)節(jié)相關(guān)蛋白編碼基因(TCONS_00038733)、磷酸烯醇丙酮酸羧激酶編碼基因(XM_006560464.2)等8個DEG均與3條信號通路存在富集關(guān)系。
利用TargetFinder軟件對ame-miR-6001-3p與本研究中的DEG進行靶向關(guān)系預(yù)測,發(fā)現(xiàn)ame-miR-6001-3p可靶向結(jié)合98個DEG,其中上調(diào)基因54個,下調(diào)基因44個(圖5)。進一步分析發(fā)現(xiàn),ame-miR-6001-3p靶向結(jié)合的上調(diào)基因富集在新陳代謝途徑(8個上調(diào)基因)、磷酸肌醇代謝(2個上調(diào)基因)、胰島素信號通路(2個上調(diào)基因)、Hippo信號通路(2個上調(diào)基因)、谷胱甘肽代謝(1個上調(diào)基因)、氧化磷酸化(1個上調(diào)基因)、鈣離子信號通路(1個上調(diào)基因)、cAMP信號通路(1個上調(diào)基因)、FoxO信號通路(1個上調(diào)基因)、AMPK信號通路(1個上調(diào)基因)和PI3K-Akt信號通路(1個上調(diào)基因)等43條代謝通路;靶向結(jié)合的下調(diào)基因富集在Hippo信號通路(3個下調(diào)基因)、新陳代謝途徑(3個下調(diào)基因)、谷胱甘肽代謝(1個下調(diào)基因)和花生四烯酸代謝(1個下調(diào)基因)等20條代謝通路。
利用RT-qPCR對隨機選擇的3個上調(diào)基因和3個下調(diào)基因進行表達量檢測,結(jié)果顯示上述6個DEG與轉(zhuǎn)錄組測序中相應(yīng)的基因表達量變化趨勢一致(圖6-A—F),證明了本研究中基因差異表達情況真實可靠。
此外,瓊脂糖凝膠電泳結(jié)果顯示ame-miR-6001-3p在Am7和Am10內(nèi)均有表達(圖6-G);RT-qPCR檢測結(jié)果顯示ame-miR-6001-3p在Am10中的相對表達量較其在Am7中的相對表達量顯著下調(diào)(= 0.0003)(圖6-H)。
腸道是蜜蜂吸收消化以及免疫防御的主要場所[16-17],腸道的發(fā)育情況及內(nèi)環(huán)境穩(wěn)態(tài)共同影響著不同蟲態(tài)蜜蜂的正常生命活動。前期研究中,筆者團隊對意蜂幼蟲腸道發(fā)育及響應(yīng)蜜蜂球囊菌()脅迫的mRNA和miRNA應(yīng)答進行了全面細致的組學(xué)研究[7,15,18-19];也分別在lncRNA和circRNA組學(xué)層面探究了意蜂工蜂中腸的發(fā)育機理[10-11]。基于前期已獲得的意蜂工蜂中腸的全轉(zhuǎn)錄組數(shù)據(jù),本研究對意蜂工蜂中腸發(fā)育過程的DEG及其潛在作用,以及融合DEG與ame-miR-6001-3p的調(diào)控網(wǎng)絡(luò)進行深入分析,發(fā)現(xiàn)Am7 vs Am10比較組的上調(diào)基因和下調(diào)基因數(shù)分別為515和523個,其中分別有285、273和258個DEG涉及細胞進程、代謝進程和催化活性等功能條目,分別有13、8、8、6、6和5個DEG富集在氧化磷酸化、磷酸肌醇代謝、碳代謝、氨基糖與核苷酸糖代謝、脂肪酸代謝和嘌呤代謝等能量和物質(zhì)代謝通路,說明意蜂工蜂中腸發(fā)育過程伴隨著較為旺盛細胞生命與新陳代謝活動。
昆蟲體內(nèi)具有高度復(fù)雜的信號通路,同一基因可能涉及多條信號通路(圖4)。有研究表明,TGF-beta信號通路主要參與調(diào)節(jié)細胞增殖、組織分化和免疫防御等生物學(xué)過程,并能與Wnt、FoxO、Hippo、MAPK、Notch等信號通路共同作用,影響昆蟲色素沉淀和附肢發(fā)育[20],蜜蜂幼蟲的腸道發(fā)育[7],以及調(diào)控蜂王的卵巢激活和產(chǎn)卵活動[21]等。本研究中,分別有14、12、9、6、3、2個DEG富集在Hippo、Wnt、FoxO、MAPK、TGF-beta、Notch信號通路,其中SMAD3亞型編碼基因(TCONS_00001946)顯著下調(diào)表達(=0.0489,log2fold change=-7.1)并同時富集在TGF-beta、Wnt、Hippo和FoxO 4條信號通路。寧越等[22]研究發(fā)現(xiàn),SMAD1可正向調(diào)控秦川牛()原代成肌細胞的分化。SMAD轉(zhuǎn)導(dǎo)因子位于TGF-beta通路的下游,在果蠅幼蟲發(fā)育過程中Activin/TGF-beta通路可通過SMAD2促進細胞的生長[23]。因此,推測SMAD3亞型編碼基因通過TGF-beta信號通路影響細胞的增殖和分化,進而調(diào)控意蜂工蜂中腸的發(fā)育過程。AMPK信號通路作為一種高度保守的細胞能量調(diào)節(jié)器,在維持細胞內(nèi)能量平衡方面發(fā)揮重要的調(diào)節(jié)作用,并可通過介導(dǎo)NF-B信號通路進而抑制具有調(diào)節(jié)昆蟲蛻皮激素分泌、海藻糖代謝及組織發(fā)育等功能的胰島素信號通路的活躍程度[24-25]。本研究中,絲裂原活化蛋白激酶亞型編碼基因(TCONS_00038549)顯著下調(diào)表達(=0.0002,log2fold change=-7.3)且同時富集在AMPK和NF-B信號通路;類胰島素受體樣編碼基因(XM_016918001.1和XM_016918002.1)顯著上調(diào)(=0.0004,log2fold change=9.2;=0.0031,log2fold change=4.2),并富集在胰島素耐受性信號通路;此外,還發(fā)現(xiàn)有3個與海藻糖轉(zhuǎn)運蛋白相關(guān)的編碼基因(XM_016916206.1、NM_001113739.1和XM_006565793.2)的表達水平顯著上調(diào)(=0.0292,log2fold change=10.0;=3.46E-39,log2fold change= 9.5;=0.0233,log2fold change=1.5)。推測絲裂原活化蛋白激酶亞型編碼基因通過下調(diào)其表達水平,減輕對胰島素信號通路的抑制作用間接增強類胰島素受體樣編碼基因的表達,從而影響促進意蜂工蜂中腸對海藻糖等能量物質(zhì)的吸收和代謝作用,進而影響中腸的發(fā)育過程。

A:AMPK信號通路AMPK signaling pathway;B:內(nèi)吞作用Endocytosis;C:P13K-Akt信號通路PI3K-Akt signaling pathway;D:Hippo信號通路Hippo signaling pathway

六邊形代表信號通路hexagons indicate signaling pathways;紅色圓形代表上調(diào)基因red circles indicate up-regulated genes;綠色圓形代表下調(diào)基因green circles indicate down-regulated genes;六邊形和圓形的大小代表連接DEG或信號通路的數(shù)量the size of hexagons and circles indicates the number of DEGs or signaling pathways connected

圖5 意蜂工蜂中腸的DEG與ame-miR-6001-3p的調(diào)控網(wǎng)絡(luò)

圖6 DEG與ame-miR-6001-3p的驗證
PI3K-Akt信號通路作為細胞內(nèi)重要信號轉(zhuǎn)導(dǎo)通路之一,通過影響下游多種效應(yīng)分子的活化狀態(tài),調(diào)節(jié)昆蟲細胞增殖、組織分化以及個體發(fā)育等多個生物學(xué)進程[26]。此外,PI3K-Akt和mTOR信號通路可共同作用,通過加強干擾素刺激基因的翻譯過程激活抗病原免疫應(yīng)激反應(yīng)[27]。本研究發(fā)現(xiàn),磷脂酰肌醇3-激酶調(diào)節(jié)亞基-類編碼基因(XM_016917733.1和XM_392119.6)、mTOR的調(diào)節(jié)相關(guān)蛋白編碼基因(TCONS_00038733)和結(jié)核蛋白亞型編碼基因(XM_016917423.1)4個DEG同時富集在PI3K-Akt和mTOR信號通路,此外,干擾素相關(guān)發(fā)育調(diào)節(jié)因子編碼基因(TCONS_00019888)呈顯著上調(diào)表達(=0.0033,log2fold change=12.4)。推測在意蜂工蜂腸道發(fā)育過程中PI3K-Akt和mTOR通路共同作用刺激干擾素相關(guān)基因的表達,從而激活中腸的抗病原免疫應(yīng)激反應(yīng),這仍有待于進一步的實驗驗證。
昆蟲腸道發(fā)育過程復(fù)雜,伴隨著非編碼RNA和編碼RNA的相互作用[28]。為探究意蜂幼蟲腸道的發(fā)育機理,筆者團隊利用二代測序和趨勢分析初步解析了意蜂幼蟲腸道發(fā)育過程中差異基因表達譜,發(fā)現(xiàn)呈顯著上調(diào)趨勢的DEG廣泛富集在核苷酸代謝、脂代謝、氨基酸代謝等各類新陳代謝通路,以及FoxO、Hippo、mTOR、MAPK、Notch、Wnt、TGF-beta、Hedgehog、Jak-STAT等信號通路上,表明隨著幼蟲腸道發(fā)育時間的延長,與發(fā)育和免疫相關(guān)通路上富集基因的表達水平逐漸增強[6]。此外,筆者團隊還發(fā)現(xiàn)在意蜂幼蟲腸道發(fā)育過程中,共有560個miRNA與16 479個靶mRNA存在靶向結(jié)合關(guān)系,宿主miRNA通過調(diào)控相關(guān)通路及富集基因參與對腸道發(fā)育和免疫等生命活動的調(diào)控[29]。本研究發(fā)現(xiàn)隨著意蜂工蜂中腸發(fā)育時間的延長,富集在AMPK、PI3K-Akt、Hippo等免疫和發(fā)育相關(guān)信號通路上的上調(diào)基因數(shù)均高于下調(diào)基因數(shù),表明上述通路存在不同程度的激活。綜上,推測意蜂幼蟲腸道和工蜂中腸隨著發(fā)育歷期的延長,其內(nèi)部與發(fā)育和免疫相關(guān)的信號通路及富集基因激活表達,以適應(yīng)腸道生長和發(fā)育所需的物質(zhì)能量代謝,并維持腸道內(nèi)環(huán)境穩(wěn)態(tài)以應(yīng)答外界復(fù)雜環(huán)境。
LncRNA和circRNA可作為ceRNA,通過miRNA應(yīng)答元件競爭性結(jié)合miRNA,以減輕miRNA對mRNA的負調(diào)控作用,間接影響基因的表達,最終調(diào)節(jié)各類生物進程[30]。郭睿等[10-11]曾對意蜂7和10日齡工蜂中腸發(fā)育過程的lncRNA和circRNA的差異表達譜及ceRNA調(diào)控網(wǎng)絡(luò)進行解析,研究結(jié)果顯示這兩種ncRNA均能通過靶向結(jié)合ame-miR-6001-3p、let-7-z和miR-136-x等對靶mRNA的表達水平進行間接調(diào)控,經(jīng)比較分析發(fā)現(xiàn),本研究中的283個DEG能被112個DElncRNA靶向的111個miRNA所調(diào)控,120個DEG可被141個DEcircRNA靶向的107個miRNA所調(diào)控。功能注釋結(jié)果顯示上述DEG參與意蜂工蜂中腸發(fā)育過程中的各類信號轉(zhuǎn)導(dǎo)及細胞生命活動。前期分析結(jié)果顯示lncRNA和circRNA也可分別通過調(diào)節(jié)上下游基因和來源基因的表達水平,參與意蜂工蜂中腸發(fā)育過程中的TGF-beta、Wnt、Hippo和Notch等信號通路,以及各類細胞生命和新陳代謝活動。上述結(jié)果表明意蜂工蜂中腸發(fā)育過程伴隨著極其復(fù)雜的分子事件,mRNA、lncRNA和circRNA可通過競爭性結(jié)合miRNA形成復(fù)雜的ceRNA調(diào)控網(wǎng)絡(luò)。前人研究發(fā)現(xiàn)ame-miR-6001-3p在蜜蜂蛻皮激素分泌[8]和級型分化[9]過程中具有調(diào)控作用。意蜂工蜂中腸發(fā)育過程中分別有29個顯著性上調(diào)lncRNA,以及14個顯著性上調(diào)circRNA參與競爭性結(jié)合ame-miR-6001- 3p[10-11]。本研究以顯著下調(diào)表達的ame-miR-6001-3p為介導(dǎo),構(gòu)建DEG與ame-miR-6001-3p的調(diào)控網(wǎng)絡(luò),發(fā)現(xiàn)ame-miR-6001-3p靶向結(jié)合54個上調(diào)基因和44個下調(diào)基因。根據(jù)ceRNA假說,推測lncRNA和circRNA通過改變表達水平影響對ame-miR-6001-3p的吸附結(jié)合,從而調(diào)節(jié)ame-miR-6001-3p的表達水平,進而部分解除ame-miR-6001-3p對靶mRNA的調(diào)控。進一步對ame-miR-6001-3p靶向結(jié)合的上調(diào)基因進行功能及代謝通路注釋,發(fā)現(xiàn)它們主要涉及胰島素信號通路(2個上調(diào)基因)、Hippo信號通路(2個上調(diào)基因)、PI3K-Akt信號通路(1個上調(diào)基因)等多條信號通路,以及磷酸肌醇代謝(2個上調(diào)基因)和谷胱甘肽代謝(1個上調(diào)基因)等各類代謝途徑。綜上,推測ceRNA網(wǎng)絡(luò)在意蜂工蜂中腸發(fā)育過程中發(fā)揮廣泛的調(diào)控作用。下一步的工作重點是結(jié)合意蜂工蜂中腸的miRNA組學(xué)數(shù)據(jù)構(gòu)建mRNA-miRNA-lncRNA- circRNA調(diào)控網(wǎng)絡(luò),進一步篩選出關(guān)鍵代謝通路、關(guān)鍵基因和關(guān)鍵ncRNA作為候選靶標,并通過過表達和敲減進行功能研究。
意蜂工蜂中腸發(fā)育過程中DEG通過參與TGF- beta、Wnt、Hippo、Notch、PI3K-Akt、mTOR、AMPK和NF-B等信號通路影響中腸的生長發(fā)育和免疫防御;ceRNA網(wǎng)絡(luò)在中腸發(fā)育中發(fā)揮廣泛的調(diào)控作用;顯著下調(diào)表達的ame-miR-6001-3p靶向結(jié)合的54個上調(diào)基因可參與胰島素信號通路等多條代謝途徑。
[1] COON K L, VOGEL K J, BROWN M R, STRAND M R. Mosquitoes rely on their gut microbiota for development., 2014, 23(11): 2727-2739.
[2] CHOUAIA B, ROSSI P, MONTAGNA M, RICCI I, CROTTI E, DAMIANI C, EPIS S, FAYE I, SAGNON N, ALMA A, FAVIA G, DAFFONCHIO D, BANDI C. Molecular evidence for multiple infections as revealed by typing ofbacterial symbionts of four mosquito species., 2010, 76(22): 7444-7450.
[3] STORELLI G, DEFAYE A, ERKOSAR B, HOLS P, ROYET J, LEULIER F.promotessystemic growth by modulating hormonal signals through TOR-dependent nutrient sensing., 2011, 14(3): 403-414.
[4] ZHENG H, POWELL J E, STEELE M I, DIETRICH C, MORAN N A. Honeybee gut microbiota promotes host weight gain via bacterial metabolism and hormonal signaling., 2017, 114(18): 4775-4780.
[5] 馮倩倩. 意大利蜜蜂春繁階段代用花粉中VA、VE適宜水平的研究[D]. 泰安: 山東農(nóng)業(yè)大學(xué), 2012.
FENG Q Q. Optimal level of VA and VE in pollen substitutes ofduring the period of spring multiplication[D]. Taian: Shandong Agricultural University, 2012. (in Chinese)
[6] 郭睿, 解彥玲, 熊翠玲, 尹偉軒, 鄭燕珍, 付中民, 陳大福. 意大利蜜蜂4、5和6日齡幼蟲腸道發(fā)育過程中差異表達基因的趨勢分析. 上海交通大學(xué)學(xué)報(農(nóng)業(yè)科學(xué)版), 2018, 36(4): 14-21, 29.
GUO R, XIE Y L, XIONG C L, YIN W X, ZHENG Y Z, FU Z M, CHEN D F. Trend analysis for differentially expressed genes in the developmental process of 4-, 5- and 6-day-old larval guts of., 2018, 36(4): 14-21, 29. (in Chinese)
[7] 郭睿, 杜宇, 熊翠玲, 鄭燕珍, 付中民, 徐國鈞, 王海朋, 陳華枝, 耿四海, 周丁丁, 石彩云, 趙紅霞, 陳大福. 意大利蜜蜂幼蟲腸道發(fā)育過程中的差異表達microRNA及其調(diào)控網(wǎng)絡(luò). 中國農(nóng)業(yè)科學(xué), 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)
[8] MELLO T R P, ALEIXO A C, PINHEIRO D G, NUNES F M F, BITONDI M M G, HARTFELDER K, BARCHUK A R, SIM?ES Z L P. Developmental regulation of ecdysone receptor (EcR) and EcR-controlled gene expression during pharate-adult development of honeybees ()., 2014, 5: 445.
[9] COLLINS D H, MOHORIANU I, BECKERS M, MOULTON V, DALMAY T, BOURKE AF. MicroRNAs associated with caste determination and differentiation in a primitively eusocial insect., 2017, 7: 45674.
[10] 郭睿, 耿四海, 熊翠玲, 鄭燕珍, 付中民, 王海朋, 杜宇, 童新宇, 趙紅霞, 陳大福. 意大利蜜蜂工蜂中腸發(fā)育過程中長鏈非編碼RNA的差異表達分析. 中國農(nóng)業(yè)科學(xué), 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)
[11] 郭睿, 陳華枝, 熊翠玲, 鄭燕珍, 付中民, 徐國鈞, 杜宇, 王海朋, 耿四海, 周丁丁, 劉思亞, 陳大福. 意大利蜜蜂工蜂中腸發(fā)育過程中的差異表達環(huán)狀RNA及其調(diào)控網(wǎng)絡(luò)分析. 中國農(nóng)業(yè)科學(xué), 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)
[12] HUANG Q, CHEN Y P, WANG R W, CHENG S, EVANS J D. Host-parasite interactions and purifying selection in a microsporidian parasite of honey bees., 2016, 11(2): e0147549.
[13] HUANG W F, SOLTER L F. Comparative development and tissue tropism ofand., 2013, 113(1): 35-41.
[14] SMOOT M E, ONO K, RUSCHEINSKI J, WANG P L, IDEKER T. Cytoscape 2.8: new features for data integration and network visualization.(, 2011, 27(3): 431-432.
[15] 郭睿, 杜宇, 童新宇, 熊翠玲, 鄭燕珍, 徐國鈞, 王海朋, 耿四海, 周丁丁, 郭意龍, 吳素珍, 陳大福. 意大利蜜蜂幼蟲腸道在球囊菌侵染前期的差異表達microRNA及其調(diào)控網(wǎng)絡(luò). 中國農(nóng)業(yè)科學(xué), 2019, 52(1): 166-180.
GUO R, DU Y, TONG X Y, XIONG C L, ZHENG Y Z, XU G J, WANG H P, GENG S H, ZHOU D D, GUO Y L, WU S Z, CHEN D F. Differentially expressed microRNAs and their regulation networks inlarval gut during the early stage ofinfection., 2019, 52(1): 166-180.(in Chinese)
[16] MARTINSON V G, MOY J, MORAN N A. Establishment of characteristic gut bacteria during development of the honeybee worker., 2012, 78(8): 2830-2840.
[17] KWONG W K, MANCENIDO A L, MORAN N A. Immune system stimulation by the native gut microbiota of honey bees., 2017, 4(2): 170003.
[18] 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.
[19] 郭睿, 杜宇, 周倪紅, 劉思亞, 熊翠玲, 鄭燕珍, 付中民, 徐國鈞, 王海朋, 耿四海, 周丁丁, 陳大福. 意大利蜜蜂幼蟲腸道在球囊菌脅迫后期的差異表達微小RNA及其靶基因分析. 昆蟲學(xué)報, 2019, 62(1): 49-60.
GUO R, DU Y, ZHOU N H, LIU S Y, XIONG C L, ZHENG Y Z, FU Z M, XU G J, WANG H P, GENG S H, ZHOU D D, CHEN D F. Comprehensive analysis of differentially expressed microRNAs and their target genes in the larval gut ofduring the late stage ofstress, 2019, 62(1): 49-60.(in Chinese)
[20] 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.
[21] 陳曉. 蜜蜂卵巢激活和產(chǎn)卵過程差異表達的編碼RNA與非編碼RNA的篩選和鑒定[D]. 北京: 中國農(nóng)業(yè)科學(xué)院, 2017.
CHEN X. Identification of differentially expressed coding and noncoding RNAs during ovary activation and oviposition in honey bees[D]. Beijing: Chinese Academy of Agricultural Sciences, 2017. (in Chinese)
[22] 寧越, 米雪, 陳星伊, 邵建航, 昝林森. SMAD1基因的沉默和過表達及對秦川牛原代成肌細胞生肌的影響. 中國農(nóng)業(yè)科學(xué), 2019, 52(10): 1818-1829.
NING Y, MI X, CHEN X Y, SHAO J H, ZAN L S. Silencing and overexpressing SMAD family member 1 () gene and its effect on myogenesis in primary myoblast of Qinchuan cattle ()., 2019, 52(10): 1818-1829.(in Chinese)
[23] BRUMMEL T, ABDOLLAH S, HAERRY T E, SHIMELL M J, MERRIAM J, RAFTERY L, WRANA J L, O’CONNOR M B. Theactivin receptor baboon signals through dSmad2 and controls cell proliferation but not patterning during larval development., 1999, 13(1): 98-111.
[24] 張振, 盧忠燕. 家蠶類胰島素信號通路及其功能的研究進展. 蠶學(xué)通訊, 2014, 34(2): 22-26.
ZHANG Z, LU Z Y. Research progress on insulin signaling pathway and its function in silkworm., 2014, 34(2): 22-26. (in Chinese)
[25] GAUTAM S, ISHRAT N, YADAV P, SINGH R, NARENDER T, SRIVASTAVA A K. 4-Hydroxyisoleucine attenuates the inflammation- mediated insulin resistance by the activation of AMPK and suppression of SOCS-3 coimmunoprecipitation with both the IR-subunit as well as IRS-1., 2016, 414(1/2): 95-104.
[26] NEUFELD T P. Body building: regulation of shape and size by PI3K/TOR signaling during development., 2003, 120(11): 1283-1296.
[27] KAUR S, SASSANO A, DOLNIAK B, JOSHI S, MAJCHRZAK- KITA B, BAKER D P, HAY N, FISH E N, PLATANIAS L C. Role of the Akt pathway in mRNA translation of interferon-stimulated genes., 2008, 105(12): 4808-4813.
[28] LIU S, LUCAS K J, ROY S, HA J, RAIKHEL A S.Mosquito-specific microRNA-1174 targets serine hydroxymethyltransferase to control key functions in the gut., 2014, 111(40): 14460-14465.
[29] 熊翠玲, 杜宇, 陳大福, 鄭燕珍, 付中民, 王海朋, 耿四海, 陳華枝, 周丁丁, 吳素珍, 石彩云, 郭睿. 意大利蜜蜂幼蟲腸道的miRNAs的生物信息學(xué)預(yù)測及分析. 應(yīng)用昆蟲學(xué)報, 2018, 55(6): 1023-1033.
XIONG C L, DU Y, CHEN D F, ZHENG Y Z, FU Z M, WANG H P, GENG S H, CHEN H Z, ZHOU D D, WU S Z, SHI C Y, GUO R. Bioinformatic prediction and analysis of miRNAs in thelarval gut., 2018, 55(6): 1023-1033.(in Chinese)
[30] KARRENTH F A, TAY Y, PERNA D, ALA U, TAN S M, RUST A G, DENICOLA G, WEBSTER K A, WEISS D, PEREZ-MANCERA P A, KRAUTHAMMER M, HALABAN R, PROVERO P, ADAMS D J, TUVESON D A, PANDOLFI P P.identification of tumor-suppressive PTEN ceRNAs in an oncogenic BRAF-induced mouse model of melanoma., 2011, 147(2): 382-395.
Profiling and Regulation Network of Differentially Expressed Genes During the Development Process ofWorker’s Midgut
DU Yu, ZHOU DingDing, WAN JieQi, LU JiaXuan, FAN XiaoXue, FAN YuanChan, CHEN Heng, XIONG CuiLing, ZHENG YanZhen, FU ZhongMin, XU GuoJun, CHEN DaFu, GUO Rui
(College of Animal Sciences (College of Bee Science), Fujian Agriculture and Forestry University, Fuzhou 350002)
【】The whole transcriptome sequencing of7- and 10-day-old workers’ midguts (Am7 and Am10) was previously conducted. In this study, the differential expression profile and regulation network of genes were investigated to reveal the molecular mechanism underlying the midgut development.【】Gene expressions were calculated based on FPKM (fragments per kilobase of transcript per million mapped reads) algorithm, and differentially expressed genes (DEGs) were gained following the standard |log2fold change|≥1 and≤0.05. Target mRNAs of ame-miR-6001-3p were predicted utilizing TargetFinder. Annotations of all DEGs in GO and KEGG databases were performed using related software. In addition, DEGs enriched in 13 signaling pathways including AMPK, P13K-Akt, Wnt, cAMP, FoxO, Hippo, mTOR, Jak-STAT, Toll-like receptor, TGF-beta, Notch, MAPK and NF-B, as well as DEGs targeted by ame-miR-6001-3p were screened out, followed by visualization of enrichment networks and regulation networks with Cytoscape. Finally, Stem loop RT-PCR and RT-qPCR were used to verify the expression and differential expression pattern of ame-miR-6001-3p and DEGs in Am7 and Am10.【】A total of 1 038 DEGs were identified in Am7 vs Am10 comparison group, including 515 up- and 523 down-regulated genes, respectively. These DEGs were associated with cellular process, metabolic process and catalytic activity, and significantly enriched in some material and energy metabolisms such as oxidative phosphorylation, amino sugar and nucleotide sugar metabolisms, fatty acid metabolism and purine metabolism, indicative of the active cellular and metabolic activities. Expression cluster analysis suggested that 20, 18, 15 and 14 DEGs were respectively enriched in AMPK signaling pathway, PI3K-Akt signaling pathway, endocytosis and Hippo signaling pathway. In addition, 57 DEGs were enriched in the aforementioned 13 signaling pathways associated with growth and development as well as immune defense, among them one DEG was enriched in several signaling pathways. Moreover, regulation network analysis showed that 54 up-regulated genes and 44 down-regulated genes were targets of ame-miR-6001-3p, respectively; these up-regulated genes were enriched in 43 pathways including inositol phosphate metabolism, Hippo signaling pathway, glutathione metabolism and insulin signaling pathway, while these down-regulated genes were enriched in 20 pathways including Hippo signaling pathway, metabolic pathways, glutathione metabolism and arachidonic acid metabolism. Moreover, RT-qPCR result showed that the variation trend of six randomly selected DEGs were consistent with that in sequencing data, confirming the reliability of DEGs. Finally, ame-miR-6001-3p was definitely expressed and significantly down-regulated in Am10.【】In this work, the expression pattern of DEGs and the regulation network between DEGs and ame-miR-6001-3p as well as the potential role of DEGs during the developmental process ofworker’s midgut were deeply analyzed. The results revealed that the DEGs may participate in various signaling pathways including TGF-beta, Wnt, Hippo, Notch, PI3K-Akt, mTOR, AMPK, NF-B signaling pathways, thus affecting the growth and development as well as immune defense of the midgut; DEGs were likely to regulate several signaling pathways such as insulin signaling pathway during the midgut development via formation regulation networks with significantly down-regulated ame-miR-6001-3p.
; midgut; development; differentially expressed gene (DEG); competitive endogenous RNA; regulation network
10.3864/j.issn.0578-1752.2020.01.019

2019-07-09;
2019-09-02
國家現(xiàn)代農(nóng)業(yè)產(chǎn)業(yè)技術(shù)體系建設(shè)專項資金(CARS-44-KXJ7)、福建省自然科學(xué)基金(2018J05042)、福建省教育廳中青年教師教育科研項目(JAT170158)、福建農(nóng)林大學(xué)杰出青年科研人才計劃(xjq201814)、福建農(nóng)林大學(xué)科技創(chuàng)新專項基金(CXZX2017342,CXZX2017343)、福建省大學(xué)生創(chuàng)新創(chuàng)業(yè)訓(xùn)練計劃(3165602032,3155006018)
杜宇,E-mail:m18505700830@163.com。周丁丁,E-mail:ZDD03569981@163.com。杜宇和周丁丁為同等貢獻作者。通信作者郭睿,E-mail:ruiguo@fafu.edu.cn
(責(zé)任編輯 岳梅)