• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    意大利蜜蜂工蜂中腸發(fā)育過(guò)程中長(zhǎng)鏈非編碼RNA的差異表達(dá)分析

    2018-10-11 02:24:08郭睿耿四海熊翠玲鄭燕珍付中民王海朋杜宇童新宇趙紅霞陳大福
    關(guān)鍵詞:意蜂中腸工蜂

    郭睿,耿四海,熊翠玲,鄭燕珍,付中民,王海朋,杜宇,童新宇,趙紅霞,陳大福

    ?

    意大利蜜蜂工蜂中腸發(fā)育過(guò)程中長(zhǎng)鏈非編碼RNA的差異表達(dá)分析

    郭睿1,耿四海1,熊翠玲1,鄭燕珍1,付中民1,王海朋1,杜宇1,童新宇1,趙紅霞2,陳大福1

    (1福建農(nóng)林大學(xué)蜂學(xué)學(xué)院,福州 350002;2廣東省生物資源應(yīng)用研究所,廣州 510260)

    【目的】長(zhǎng)鏈非編碼RNA(lncRNA)在真核生物的基因表達(dá)、表觀遺傳和細(xì)胞周期調(diào)控等方面發(fā)揮重要功能。本研究旨在探究意大利蜜蜂(, 簡(jiǎn)稱意蜂)工蜂中腸發(fā)育過(guò)程中l(wèi)ncRNA的表達(dá)譜及其作用。【方法】利用RNA-seq技術(shù)和鏈特異性建庫(kù)方法對(duì)意蜂7和10日齡工蜂中腸(Am7、Am10)進(jìn)行深度測(cè)序,下機(jī)的原始數(shù)據(jù)經(jīng)過(guò)Perl腳本過(guò)濾得到高質(zhì)量有效讀段。利用bowtie工具將有效讀段比對(duì)核糖體數(shù)據(jù)庫(kù),進(jìn)一步利用TopHat2軟件將未比對(duì)到核糖體數(shù)據(jù)庫(kù)上的數(shù)據(jù)比對(duì)到參考基因組。利用CPC和CNCI軟件對(duì)轉(zhuǎn)錄本的編碼能力進(jìn)行預(yù)測(cè)。通過(guò)RT-PCR對(duì)部分lncRNA進(jìn)行鑒定。利用edgeR軟件進(jìn)行差異表達(dá)lncRNA(DElncRNA)分析,進(jìn)而預(yù)測(cè)lncRNA的上下游基因,并對(duì)上下游基因進(jìn)行GO及KEGG代謝通路富集分析。聯(lián)用RNAhybrid、Miranda和TargetScan軟件預(yù)測(cè)DElncRNA靶向結(jié)合的miRNA及miRNA靶向結(jié)合的靶基因,并通過(guò)Cytoscape軟件構(gòu)建DElncRNAs-miRNAs-mRNAs的調(diào)控網(wǎng)絡(luò)。最后,通過(guò)RT-qPCR驗(yàn)證測(cè)序數(shù)據(jù)的可靠性?!窘Y(jié)果】Am7和Am10的深度測(cè)序分別獲得134 802 058和147 051 470條原始讀段,經(jīng)嚴(yán)格過(guò)濾分別得到134 166 157和146 293 288條有效讀段;共得到3 890個(gè)DElncRNA,包括2 005個(gè)上調(diào)lncRNA與1 885個(gè)下調(diào)lncRNA。RT-PCR驗(yàn)證結(jié)果顯示共有8個(gè)lncRNA能擴(kuò)增出符合預(yù)期的目的片段,表明預(yù)測(cè)出的lncRNA真實(shí)存在。DElncRNA的上下游基因數(shù)為1 793個(gè),它們涉及42個(gè)GO條目,包括代謝進(jìn)程、發(fā)育進(jìn)程、細(xì)胞進(jìn)程、應(yīng)激反應(yīng)和免疫系統(tǒng)進(jìn)程等;這些上下游基因還涉及251條代謝通路,包括碳代謝、嘌呤代謝和脂肪酸的生物合成等物質(zhì)代謝通路,硫代謝、甲烷代謝和氧化磷酸化等能量代謝通路,Hippo信號(hào)通路、Wnt信號(hào)通路和Notch信號(hào)通路等信號(hào)通路,溶酶體、內(nèi)吞作用和泛素介導(dǎo)的蛋白水解等細(xì)胞免疫通路,以及MAPK信號(hào)通路、Jak-STAT信號(hào)通路和NF-kappa B信號(hào)通路等體液免疫通路,上述結(jié)果表明DElncRNA在意蜂中腸發(fā)育過(guò)程中參與物質(zhì)和能量代謝、細(xì)胞生命活動(dòng)和免疫調(diào)控。進(jìn)一步分析發(fā)現(xiàn)TCONS_00020918可通過(guò)調(diào)控西方蜜蜂王漿主蛋白1編碼基因在意蜂工蜂中腸的營(yíng)養(yǎng)吸收、級(jí)型分化中發(fā)揮功能。DElncRNA的調(diào)控網(wǎng)絡(luò)分析結(jié)果顯示DElncRNA與miRNA、mRNA間存在復(fù)雜的調(diào)控關(guān)系,部分DElncRNA處于調(diào)控網(wǎng)絡(luò)的中心位置且能靶向結(jié)合較多的miRNA,也有部分miRNA可被多個(gè)DElncRNA共同靶向,表明這些DElncRNA可能在中腸發(fā)育中發(fā)揮重要作用。隨機(jī)挑取5個(gè)DElncRNA進(jìn)行RT-qPCR驗(yàn)證,結(jié)果顯示它們的表達(dá)量變化趨勢(shì)與測(cè)序結(jié)果一致,證實(shí)了本研究測(cè)序數(shù)據(jù)的可靠性?!窘Y(jié)論】差異表達(dá)長(zhǎng)鏈非編碼RNA(DElncRNA)廣泛參與意蜂工蜂中腸的新陳代謝、細(xì)胞活動(dòng)和免疫調(diào)控并作為競(jìng)爭(zhēng)性內(nèi)源RNA(ceRNA)發(fā)揮作用,研究結(jié)果為關(guān)鍵lncRNA的篩選和功能研究提供了必要的數(shù)據(jù)支持。

    意大利蜜蜂;中腸;發(fā)育;長(zhǎng)鏈非編碼RNA;上下游基因

    0 引言

    【研究意義】蜜蜂作為社會(huì)學(xué)模式昆蟲(chóng),在發(fā)育學(xué)、行為學(xué)和神經(jīng)生物學(xué)等方面具有重要價(jià)值[1]。同時(shí),蜜蜂也是最重要的授粉昆蟲(chóng),具有不可替代的經(jīng)濟(jì)價(jià)值和生態(tài)價(jià)值[2]。意大利蜜蜂(,簡(jiǎn)稱意蜂)具有優(yōu)越的采集能力、造脾能力和分泌蜂王漿能力[3],在世界各地的養(yǎng)蜂生產(chǎn)中廣泛使用。目前,有關(guān)蜜蜂中腸發(fā)育機(jī)理及調(diào)控機(jī)制的研究極為滯后,非編碼RNA(non-coding RNA,ncRNA)在中腸發(fā)育過(guò)程中的作用的相關(guān)信息也十分有限。在轉(zhuǎn)錄組水平對(duì)長(zhǎng)鏈非編碼RNA(long non-coding RNA,lncRNA)在意蜂工蜂中腸發(fā)育中的作用進(jìn)行研究,可為深入解析意蜂工蜂中腸發(fā)育的分子機(jī)理提供新的思路和線索?!厩叭搜芯窟M(jìn)展】隨著測(cè)序技術(shù)的發(fā)展,起初被認(rèn)為基因轉(zhuǎn)錄“噪音”的ncRNA成為近幾年的研究熱點(diǎn)。在人類基因組中,編碼蛋白序列所占比例不到2%,超過(guò)98%的序列是不編碼蛋白質(zhì)的[4]。ncRNA根據(jù)其長(zhǎng)短和形狀分為小RNA、lncRNA和環(huán)狀RNA等。其中,lncRNA已被證明在劑量補(bǔ)償效應(yīng)[5]、表觀遺傳[6]、細(xì)胞周期[7]和細(xì)胞分化調(diào)控[8]等眾多生命活動(dòng)中發(fā)揮重要作用,因此成為生命科學(xué)領(lǐng)域的研究熱點(diǎn)。lncRNA的作用方式多樣,包括作為信號(hào)分子調(diào)控上下游基因轉(zhuǎn)錄[9],作為誘導(dǎo)分子充當(dāng)分子阻斷劑[10],作為引導(dǎo)分子與蛋白結(jié)合[11],作為支架分子發(fā)揮“腳手架”作用[12],作為小RNA的前體[13],以及作為競(jìng)爭(zhēng)內(nèi)源RNA結(jié)合miRNA對(duì)mRNA起剪接調(diào)控作用[14]。近年來(lái),有關(guān)lncRNA在腫瘤發(fā)生機(jī)制方面作用的研究取得了顯著進(jìn)展,在脂肪代謝、肌肉發(fā)育及免疫抗病等機(jī)理方面作用的研究也取得了許多突破性進(jìn)展[15-16]。然而,lncRNA在昆蟲(chóng)領(lǐng)域的相關(guān)研究較為滯后,僅在果蠅和家蠶等極少數(shù)模式昆蟲(chóng)中有較少的lncRNA相關(guān)研究報(bào)道[17]。此前,有研究報(bào)道lncRNA在蜜蜂的級(jí)型分化[18]、卵巢發(fā)育[19-20]和抵抗病毒入侵[21]中起重要的調(diào)控作用,但總體上有關(guān)蜜蜂lncRNA的研究仍非常有限。腸道是昆蟲(chóng)的重要消化器官和免疫器官。蜜蜂腸道的相關(guān)研究主要集中在腸道微生物方面[22-27],如Engel等[26]利用二代測(cè)序技術(shù)對(duì)意蜂工蜂后腸進(jìn)行宏基因組測(cè)序,基于數(shù)據(jù)分析發(fā)現(xiàn)意蜂后腸內(nèi)的少量細(xì)菌種類存在極大的遺傳多樣性,進(jìn)一步的比較分析發(fā)現(xiàn)不同的細(xì)菌種類涉及不同的功能,包括宿主互作、菌膜形成及碳水化合物水解等。【本研究切入點(diǎn)】蜜蜂中腸發(fā)育機(jī)理及調(diào)控機(jī)制仍不明確,關(guān)于lncRNA在蜜蜂中腸發(fā)育過(guò)程中作用的研究未見(jiàn)報(bào)道?!緮M解決的關(guān)鍵問(wèn)題】結(jié)合RNA-seq技術(shù)和鏈特異性cDNA建庫(kù)方法對(duì)意蜂7和10日齡的工蜂中腸進(jìn)行測(cè)序,對(duì)中腸發(fā)育過(guò)程的lncRNA進(jìn)行差異表達(dá)分析,進(jìn)一步通過(guò)上下游基因分析和lncRNA-miRNA-mRNA調(diào)控網(wǎng)絡(luò)分析DElncRNA在意蜂工蜂中腸發(fā)育過(guò)程中的作用。

    1 材料與方法

    試驗(yàn)于2017年在福建農(nóng)林大學(xué)蜂學(xué)學(xué)院蜜蜂保護(hù)實(shí)驗(yàn)室進(jìn)行。

    1.1 生物材料

    供試意大利蜜蜂工蜂取自福建農(nóng)林大學(xué)蜂學(xué)學(xué)院教學(xué)蜂場(chǎng)。

    1.2 樣品制備與Illumina測(cè)序

    從群勢(shì)較強(qiáng)的意蜂蜂群中選取有封蓋子的巢脾,快速提至實(shí)驗(yàn)室,放入(34±0.5)℃培養(yǎng)箱中培養(yǎng)至工蜂出房,將剛出房的工蜂(記為0 d)放入干凈的四周打孔以通風(fēng)的塑料盒中,盒子上方插入一支裝有50%(w/v)無(wú)菌糖水的飼喂器(圖1),(34±0.5)℃條件下恒溫培養(yǎng),每日檢查工蜂存活情況,及時(shí)清理死亡的工蜂。進(jìn)行3次生物學(xué)重復(fù)。

    圖1 意蜂工蜂的人工飼養(yǎng)

    選取意蜂7和10日齡工蜂中腸作為測(cè)序材料。分別快速拉取意蜂7日齡工蜂中腸(Am7)和10日齡工蜂中腸(Am10),液氮速凍后置于-80℃超低溫冰箱保存?zhèn)溆?,每只中腸的取樣時(shí)間嚴(yán)格控制在15 s以內(nèi)。首先用RNA抽提試劑盒(AxyPrepTMMultisource Total RNA Miniprep Kit)(TaKaRa公司,日本)抽提蜜蜂中腸樣品的總RNA,為最大限度地保留所有非編碼RNA(ncRNA),去除核糖體RNA后的mRNA和ncRNA用裂解緩沖液隨機(jī)打斷為小片段,作為模板用六堿基隨機(jī)引物、緩沖液、dNTPs、RNase H和DNA polymerase I合成cDNA第二鏈。經(jīng)過(guò)QiaQuick PCR試劑盒(Qiagen公司,德國(guó))純化并加EB緩沖液洗脫經(jīng)末端修復(fù)、加堿基A,加測(cè)序接頭,然后通過(guò)尿嘧啶-N-糖基化酶(UNG)降解cDNA第二鏈。消化產(chǎn)物經(jīng)瓊脂糖凝膠電泳和PCR擴(kuò)增,PCR產(chǎn)物經(jīng)高碘酸鈉(Sigma 公司,美國(guó))處理。委托廣州基迪奧生物科技有限公司對(duì)上述cDNA文庫(kù)進(jìn)行雙端(paired-end)測(cè)序,測(cè)序平臺(tái)為Illumina HiSeq 4000。測(cè)序數(shù)據(jù)已上傳NCBI SRA數(shù)據(jù)庫(kù),BioProject號(hào):PRJNA406998。

    1.3 高通量數(shù)據(jù)分析

    1.3.1 LncRNA的生物信息學(xué)預(yù)測(cè) 對(duì)于下機(jī)數(shù)據(jù),利用Perl腳本去除含有adaptor、未知核苷酸比例>10%和低質(zhì)量reads,獲得有效讀段(clean reads)。使用短reads比對(duì)工具bowtie[28]分別將樣品Am7-1、Am7-2、Am7-3、Am10-1、Am10-2和Am10-3的clean reads比對(duì)(mapping)到核糖體數(shù)據(jù)庫(kù)(最多允許5個(gè)錯(cuò)配),去除mapping上核糖體的reads,利用TopHat2軟件[29]將保留下來(lái)的數(shù)據(jù)進(jìn)一步比對(duì)西方蜜蜂的參考基因組(Amel_4.5)[30]。利用FPKM(Fragments Per Kilobase of transcript per Million mapped reads)法計(jì)算基因表達(dá)量。利用R軟件(version 2.16.2)計(jì)算各樣品之間的相關(guān)性系數(shù)。

    lncRNA與mRNA相比在序列保守型、ORF長(zhǎng)度等特點(diǎn)上具有一定差異,利用軟件CPC[31]和CNCI[32]對(duì)新轉(zhuǎn)錄本的編碼能力進(jìn)行預(yù)測(cè),選取二者的交集作為可靠的預(yù)測(cè)結(jié)果。

    1.3.2 DElncRNA及其上下游基因分析 利用edgeR軟件[33]進(jìn)行l(wèi)ncRNA的差異分析,得到DElncRNA。LncRNA的功能與其編碼基因位置臨近的蛋白編碼基因關(guān)系密切,位于上游的lncRNA可能與啟動(dòng)子或共表達(dá)基因的其他順式作用元件存在交集,從而在轉(zhuǎn)錄或轉(zhuǎn)錄后水平進(jìn)行基因表達(dá)調(diào)控;位于3′ UTR或基因下游的lncRNA可能參與其他調(diào)控作用。因此,對(duì)lncRNA進(jìn)行注釋,如果其位于一個(gè)基因的上游或下游,這些lncRNA有可能與順式作用元件所在區(qū)域有交集,從而參與轉(zhuǎn)錄調(diào)控。利用Omicshare平臺(tái)(http://www.omicshare.com/tools/index.php/)對(duì)DElncRNA的上下游基因進(jìn)行GO(Gene Ontology)分類和KEGG pathway富集分析。

    1.3.3 DElncRNA的調(diào)控網(wǎng)絡(luò)構(gòu)建 LncRNA可以競(jìng)爭(zhēng)結(jié)合miRNA,從而減少miRNA結(jié)合mRNA,表明lncRNA可以通過(guò)miRNA調(diào)控mRNA的表達(dá)。利用RNAhybrid[34]、Miranda[35]和TargetScan[36]軟件預(yù)測(cè)DElncRNA靶向結(jié)合的miRNA,以及miRNA結(jié)合的mRNA,根據(jù)上述靶向結(jié)合關(guān)系構(gòu)建lncRNA-miRNA- mRNA的調(diào)控網(wǎng)絡(luò)。利用Cytoscape軟件[37]對(duì)上述調(diào)控網(wǎng)絡(luò)進(jìn)行可視化。

    1.4 lncRNA的RT-PCR驗(yàn)證

    隨機(jī)選取9個(gè)lncRNA,根據(jù)它們的序列利用DNAMAN軟件設(shè)計(jì)相應(yīng)的特異性引物。委托上海生工生物工程有限公司合成引物。利用RNA抽提試劑盒(TaKaRa,中國(guó))提取Am7和Am10的總RNA,作為模板進(jìn)行反轉(zhuǎn)錄,得到的cDNA作為模板進(jìn)行PCR。PCR反應(yīng)體系(20 μL)包括cDNA模板1 μL,上游引物1 μL,下游引物1 μL,Mixture 10 μL,無(wú)菌水補(bǔ)至20 μL。PCR程序:94℃預(yù)變性5 min;94℃變性50 s,55℃退火30 s,72℃延伸50 s,共34個(gè)循環(huán);72℃再延伸10 min。PCR產(chǎn)物經(jīng)1.5%瓊脂糖凝膠電泳和凝膠成像儀(上海培清,中國(guó))檢測(cè)。

    1.5 DElncRNA的實(shí)時(shí)熒光定量PCR(RT-qPCR)驗(yàn)證

    為了驗(yàn)證RNA-seq數(shù)據(jù)的可靠性,隨機(jī)選取5個(gè)DElncRNA進(jìn)行RT-qPCR驗(yàn)證。利用DNAMAN軟件設(shè)計(jì)相應(yīng)的特異性引物。委托上海生工生物工程有限公司合成引物。利用RNA抽提試劑盒(TaKaRa,中國(guó))提取Am7和Am10的總RNA,作為模板進(jìn)行反轉(zhuǎn)錄,得到的cDNA作為模板進(jìn)行PCR。RT-qPCR反應(yīng)按照SYBR Green Dye試劑盒(Vazyme公司,中國(guó))操作說(shuō)明書(shū)進(jìn)行,每個(gè)反應(yīng)進(jìn)行3次重復(fù)。反應(yīng)體系(20 μL)包含正、反向引物(10.0 μmol·L-1)各1 μL,cDNA模板DNA 1 μL,SYBR Green Dye 10 μL,DEPC水7 μL。qRT-PCR反應(yīng)在ABI 7500熒光定量PCR儀(ABI公司,美國(guó))上進(jìn)行,反應(yīng)條件:95℃預(yù)變性1 min,95℃變性15 s,60℃延伸30 s,共40個(gè)循環(huán),最后72℃延伸45 s。利用2-ΔΔCt法對(duì)上述基因的相對(duì)表達(dá)量進(jìn)行計(jì)算。

    2 結(jié)果

    2.1 數(shù)據(jù)質(zhì)控與評(píng)估

    在實(shí)驗(yàn)室條件下對(duì)意蜂工蜂進(jìn)行人工飼養(yǎng)(圖1),Am7和Am10中腸樣品的測(cè)序得到raw reads平均分別為134 802 058和147 051 470條,過(guò)濾后得到clean reads平均分別為134 166 157和146 293 288條。各樣品的平均Q20和Q30分別為97.34%和93.86%(表1)。Am7與Am10的組內(nèi)各生物學(xué)重復(fù)之間的Pearson相關(guān)系數(shù)均在0.97以上,說(shuō)明各樣品的重復(fù)性較好(圖2)。上述結(jié)果說(shuō)明本研究的測(cè)序數(shù)據(jù)質(zhì)量良好,可用于進(jìn)一步分析。

    隨機(jī)挑選9個(gè)lncRNA進(jìn)行RT-PCR驗(yàn)證,電泳結(jié)果顯示其中有8個(gè)lncRNA均能擴(kuò)增出符合預(yù)期的目的片段,說(shuō)明本研究預(yù)測(cè)出的絕大多數(shù)的lncRNA真實(shí)存在(圖3)。相關(guān)引物信息詳見(jiàn)表2。進(jìn)一步對(duì)上述8個(gè)驗(yàn)證的lncRNA的上下游基因進(jìn)行預(yù)測(cè),其中TCONS_00020918可調(diào)控西方蜜蜂王漿主蛋白1編碼基因,TCONS_00021005可調(diào)控西方蜜蜂未知蛋白LOC725PPP,TCONS_00025221可調(diào)控西方蜜蜂未知蛋白LOC726321、LOC100577788和LOC102656594。

    橫縱坐標(biāo)均表示基因表達(dá)量(FPKM)The horizontal and vertical coordinates represent gene expression level (FPKM)

    表1 RNA-seq數(shù)據(jù)概覽

    M: DNA marker; 1: TCONS_00020918; 2: TCONS_00021005; 3: TCONS_00019675; 4: TCONS_00019678; 5: TCONS_00025221; 6: TCONS_00025232; 7: TCONS_00025235; 8: TCONS_00025236

    表2 RT-PCR與RT-qPCR引物信息

    2.2 DElncRNA及其上下游基因分析

    LncRNA的差異表達(dá)分析結(jié)果顯示,在Am7 vs Am10比較組中共有3 890個(gè)DElncRNA,包括2 005個(gè)上調(diào)lncRNA和1885個(gè)下調(diào)lncRNA。lncRNA可通過(guò)影響其編碼基因的上下游基因而發(fā)揮作用。共預(yù)測(cè)出1 698個(gè)DElncRNA的上下游基因,其中上游基因、下游基因和重疊基因分別為485、535和593個(gè)。對(duì)DElncRNA編碼基因的上下游基因進(jìn)行GO分類,結(jié)果顯示它們分布在生物學(xué)進(jìn)程、細(xì)胞組分和分子功能3大類,涉及42個(gè)GO term,富集基因數(shù)最多的前10個(gè)GO term分別為結(jié)合(355 genes)、細(xì)胞進(jìn)程(315 genes)、代謝進(jìn)程(299 genes)、單細(xì)胞進(jìn)程(262 genes)、催化活性(259 genes)、細(xì)胞(147 genes)、細(xì)胞組分(147 genes)、細(xì)胞膜(138 genes)、細(xì)胞膜組分(131 genes)和細(xì)胞器(117 genes)(圖4)。

    進(jìn)一步對(duì)DElncRNA編碼基因的上下游基因進(jìn)行KEGG pathway富集分析,結(jié)果顯示它們富集在251個(gè)pathway,富集基因數(shù)最多的前10位pathway是Hippo信號(hào)通路(36 genes)、Wnt信號(hào)通路(19 genes)、癌癥microRNA(16 genes)、癌癥通路(16 genes)、剪接體(15 genes)、PI3K-AKt信號(hào)通路(14 genes)、碳代謝(12 genes)、單純皰疹感染(12 genes)、癌癥蛋白聚糖(12 genes)和神經(jīng)營(yíng)養(yǎng)因子信號(hào)通路(11 genes)(圖5-A)。Hippo信號(hào)通路的概貌如圖5-B所示。

    2.3 DElncRNA的調(diào)控網(wǎng)絡(luò)分析

    LncRNA可以作為一種競(jìng)爭(zhēng)性內(nèi)源RNA,可通過(guò)結(jié)合miRNA減少靶向結(jié)合mRNA的miRNA數(shù),從而減輕miRNA對(duì)mRNA的抑制作用[38]。利用軟件預(yù)測(cè)DElncRNA靶向結(jié)合的miRNA,以及miRNA靶向結(jié)合的mRNA,其中靶向結(jié)合miRNA最多的DElncRNA為TCONS_00004891(36 miRNAs)、XR_001704571.1(28 miRNAs)和TCONS_00024939(27 miRNAs),靶向結(jié)合mRNA最多的miRNA為ame-miR-6001-3p(98 mRNAs)、mir-136-y(61 mRNAs)和mir-410-y(38 mRNAs)。進(jìn)一步利用Cytoscape軟件對(duì)DElncRNA-miRNA-mRNA調(diào)控網(wǎng)絡(luò)進(jìn)行可視化,結(jié)果顯示上調(diào)、下調(diào)lncRNA與miRNA和mRNA間均形成復(fù)雜的調(diào)控網(wǎng)絡(luò)(圖6),共有125個(gè)lncRNA預(yù)測(cè)到靶向結(jié)合的miRNA,其中僅有6個(gè)lncRNA靶向結(jié)合1個(gè)miRNA,其余的119個(gè)lncRNA均靶向結(jié)合2個(gè)或更多個(gè)miRNA。

    圖4 DElncRNA上下游基因的GO分類

    2.4 DElncRNA的qRT-PCR驗(yàn)證

    隨機(jī)挑取6個(gè)DElncRNA進(jìn)行RT-qPCR驗(yàn)證,結(jié)果顯示其中有5個(gè)DElncRNA表達(dá)水平的變化趨勢(shì)和轉(zhuǎn)錄組數(shù)據(jù)中相應(yīng)DElncRNA表達(dá)水平的變化趨勢(shì)一致(圖7),說(shuō)明本研究中的測(cè)序數(shù)據(jù)真實(shí)可靠。相關(guān)引物序列信息詳見(jiàn)表2。對(duì)驗(yàn)證的5個(gè)DElncRNA進(jìn)行上下游基因分析,發(fā)現(xiàn)TCONS_00029692可調(diào)控肥大細(xì)胞脫粒肽編碼基因,TCONS_00038012可調(diào)控未知蛋白LOC100576733編碼基因。

    3 討論

    此前,蜜蜂腸道的相關(guān)研究集中在腸道微生物方面,有關(guān)蜜蜂腸道的發(fā)育機(jī)理、ncRNA在腸道發(fā)育過(guò)程中作用的研究未見(jiàn)報(bào)道。前人研究結(jié)果表明lncRNA與生物的神經(jīng)發(fā)育、細(xì)胞組織發(fā)育、細(xì)胞周期調(diào)控、干細(xì)胞分化、免疫應(yīng)答和癌癥發(fā)生等過(guò)程關(guān)系密切。近年來(lái),隨著二代測(cè)序技術(shù)的迅速發(fā)展和廣泛應(yīng)用,越來(lái)越多的lncRNA在動(dòng)物、植物和微生物中被鑒定出來(lái),但蜜蜂lncRNA的研究較為滯后,相關(guān)信息極為有限。為探究意蜂工蜂中腸的發(fā)育機(jī)理及調(diào)控機(jī)制,本研究結(jié)合RNA-seq技術(shù)和鏈特異性建庫(kù)方法對(duì)意蜂7和10日齡工蜂中腸進(jìn)行測(cè)序,基于高質(zhì)量的測(cè)序數(shù)據(jù)共預(yù)測(cè)出6 353個(gè)lncRNA,RT-PCR驗(yàn)證的擴(kuò)增成功率高達(dá)88.9%,表明本研究預(yù)測(cè)出的絕大多數(shù)lncRNA真實(shí)存在。相對(duì)于非特異性建庫(kù),特異性建庫(kù)的方法在合成cDNA第二鏈時(shí),將dTTP替換為dUTP,加上接頭后再用UNG酶降解第二鏈,因而該方法可以確定轉(zhuǎn)錄本來(lái)自正鏈或負(fù)鏈,故能更準(zhǔn)確地獲得基因的結(jié)構(gòu)和基因表達(dá)信息。lncRNA的表達(dá)具有時(shí)間和組織特異性,本研究的測(cè)序?qū)ο笫且夥涔し渲心c組織,它們?cè)谄渌M織中是否表達(dá)以及表達(dá)水平有待于進(jìn)一步研究。本研究中人工飼養(yǎng)工蜂的過(guò)程中僅飼喂糖水,工蜂的腸道內(nèi)會(huì)有少量流體且多聚集在后腸,為排除雜質(zhì)對(duì)測(cè)序結(jié)果的影響,舍棄了后腸,將中腸組織用于測(cè)序。測(cè)得的原始數(shù)據(jù)包含工蜂本身的數(shù)據(jù)和少量的中腸內(nèi)微生物的數(shù)據(jù),為排除后者的干擾,對(duì)原始數(shù)據(jù)進(jìn)行了嚴(yán)格的過(guò)濾和質(zhì)控,先比對(duì)核糖體數(shù)據(jù)庫(kù)以去除rRNA數(shù)據(jù),再比對(duì)西方蜜蜂基因組去除未比對(duì)上的數(shù)據(jù),保留下來(lái)的能比對(duì)上的數(shù)據(jù)用于后續(xù)的生物信息學(xué)分析,考慮到昆蟲(chóng)和微生物的物種親緣關(guān)系較遠(yuǎn),二者的基因保守性很低,因而比對(duì)上參考基因組的數(shù)據(jù)理論上皆為意蜂工蜂中腸本身的數(shù)據(jù)。

    A:Pathway富集分析enrichment analysis;B:Hippo信號(hào)通路概貌General picture of Hippo signaling pathway

    A:上調(diào)lncRNA的lncRNA-miRNA-mRNA網(wǎng)絡(luò)lncRNA-miRNA-mRNA networks of up-regulated lncRNAs;B:下調(diào)lncRNA的lncRNA-miRNA-mRNA網(wǎng)絡(luò) lncRNA-miRNA-mRNA networks of down-regulated lncRNAs

    A: XR_409562.2; B: TCONS_00011956; C: TCONS_00012589; D: TCONS_00038012; E: TCONS_00029692

    LncRNA的功能與其編碼基因坐標(biāo)臨近的蛋白編碼基因相關(guān),位于上游的lncRNA可能與啟動(dòng)子或共表達(dá)基因的其他順式作用元件有交集,從而在轉(zhuǎn)錄或轉(zhuǎn)錄后水平對(duì)基因的表達(dá)進(jìn)行調(diào)控[39]。本研究中,DElncRNA的上下游基因達(dá)到1 697個(gè),GO富集分析結(jié)果顯示分別有299個(gè)和4個(gè)上下游基因涉及代謝進(jìn)程和發(fā)育進(jìn)程,有259個(gè)上下游基因涉及氧化磷酸化,表明DElncRNA參與意蜂工蜂中腸的物質(zhì)和能量代謝過(guò)程調(diào)控。本研究發(fā)現(xiàn),315個(gè)上下游基因涉及細(xì)胞進(jìn)程,表明DElncRNA廣泛參與意蜂工蜂中腸的細(xì)胞生命活動(dòng)。此外,還發(fā)現(xiàn)分別有92個(gè)和3個(gè)上下游基因涉及應(yīng)激反應(yīng)和免疫系統(tǒng)進(jìn)程,暗示相應(yīng)的DElncRNA在意蜂工蜂中腸發(fā)育中發(fā)揮調(diào)控作用。王漿主蛋白(MRJP)為幼蟲(chóng)發(fā)育提供必要的可利用氮元素,在蜜蜂行為中有潛在的調(diào)控作用,還參與決定蜂群中雌性蜂的級(jí)型分化等,其中MRJP1在蜂王漿主蛋白中含量最豐富,所占比例達(dá)31%,占蜂王漿中水溶性蛋白的48%[40]。本研究發(fā)現(xiàn),TCONS_00020918可調(diào)控西方蜜蜂MRJP1編碼基因,表明該lncRNA可能在意蜂的營(yíng)養(yǎng)吸收、級(jí)型分化中發(fā)揮重要的調(diào)控功能。此外,多個(gè)DElncRNA可能調(diào)控未知蛋白編碼基因,具體調(diào)控關(guān)系的明確有賴于Nr數(shù)據(jù)庫(kù)蛋白功能注釋信息的完善及進(jìn)一步的實(shí)驗(yàn)驗(yàn)證。

    蜜蜂腸道是物質(zhì)和能量合成與代謝的主要場(chǎng)所。本研究中,DElncRNA上下游基因的KEGG富集分析結(jié)果顯示,共有221個(gè)上下游基因富集在多達(dá)65個(gè)物質(zhì)代謝通路上,包括甘油酯代謝(5 genes)和脂肪酸的生物合成(2 genes)等10條脂代謝相關(guān)通路,檸檬酸鹽循環(huán)(6 genes)和丙酮酸鹽代謝(2 genes)等13條碳水化合物代謝相關(guān)通路,賴氨酸降解(3 genes)和精氨酸生物合成(2 genes)等13條氨基酸代謝相關(guān)通路,以及嘌呤代謝(10 genes)和嘧啶代謝(1 gene),表明DElncRNA廣泛參與意蜂工蜂中腸發(fā)育過(guò)程中的物質(zhì)代謝調(diào)控。物質(zhì)合成與代謝往往伴隨著旺盛的能量代謝,本研究發(fā)現(xiàn)共有15個(gè)上下游基因富集在甲烷代謝(5 genes)、氧化磷酸化(5 genes)和硫代謝(2 genes)等能量代謝通路,表明相應(yīng)的DElncRNA在意蜂工蜂中腸的能量代謝方面具有重要的調(diào)控功能。蜜蜂腸道是物質(zhì)消化吸收的主要位置。本研究發(fā)現(xiàn)28個(gè)上下游基因富集在7個(gè)消化系統(tǒng)通路上,包括胰腺分泌(8 genes)、蛋白消化和吸收(2 genes)、胃酸分泌(7 genes)、脂肪體消化和吸收(4 genes)、唾液分泌(4 genes)、維生素消化和吸收(1 genes)和膽汁分泌(2 genes),表明DElncRNA參與調(diào)控意蜂工蜂中腸對(duì)物質(zhì)的消化和吸收。Hippo信號(hào)通路可以調(diào)節(jié)器官的大小[41-42],也能與其他信號(hào)通路相互作用共同調(diào)節(jié)中腸組織的穩(wěn)態(tài)[43]。Notch信號(hào)影響細(xì)胞正常形態(tài)發(fā)生的多個(gè)過(guò)程,包括多能細(xì)胞分化、細(xì)胞凋亡、細(xì)胞增殖及細(xì)胞邊界的形成[44]。本研究發(fā)現(xiàn),分別有36、19和8個(gè)上下游基因分別富集在Hippo、Wnt和Notch信號(hào)通路。此外,還發(fā)現(xiàn)7、3、3和7個(gè)上下游基因富集在TGF-beta、mTOR、Hedgehog和胰島素信號(hào)通路。結(jié)果表明相應(yīng)的DElncRNA可通過(guò)參與調(diào)控上述信號(hào)通路對(duì)意蜂工蜂中腸發(fā)育進(jìn)行調(diào)控。

    角質(zhì)層和圍食膜是昆蟲(chóng)免疫防御的第一道防線[45]。如果病原微生物突破第一道防線,包括蜜蜂在內(nèi)的昆蟲(chóng)將會(huì)啟動(dòng)細(xì)胞和體液免疫抵抗病原的入侵。蜜蜂的細(xì)胞免疫包括內(nèi)吞作用、黑化作用、吞噬作用、蛋白的酶促水解等,體液免疫主要為抗菌肽的合成與釋放[46]。本研究中,DElncRNA上下游基因的KEGG pathway富集分析結(jié)果顯示分別有11、9、7、5和4個(gè)上下游基因富集在溶酶體、內(nèi)吞作用、泛素蛋白水解、吞噬體和黑化作用等細(xì)胞免疫通路,此外,分別有13、3、1和1個(gè)上下游基因富集在MAPK、Jak-STAT、NF-Kappa B和Toll-like受體等體液免疫通路。上述結(jié)果表明相應(yīng)的DElncRNA同時(shí)參與了宿主的細(xì)胞和體液免疫的調(diào)控過(guò)程,推測(cè)其在宿主的免疫防御過(guò)程發(fā)揮關(guān)鍵作用。

    Salmena等提出內(nèi)源性競(jìng)爭(zhēng)RNA(ceRNA)假說(shuō)[38],該假說(shuō)認(rèn)為lncRNA和mRNA等RNA可以通過(guò)miRNA應(yīng)答原件(MRE)與miRNA競(jìng)爭(zhēng)性結(jié)合,從而對(duì)基因表達(dá)進(jìn)行調(diào)控。此后,ceRNA假說(shuō)已被越來(lái)越多的研究[47]所證明。為進(jìn)一步揭示DElncRNA的作用,本研究通過(guò)靶向關(guān)系預(yù)測(cè)DElncRNA結(jié)合的miRNA以及miRNA結(jié)合的mRNA,并構(gòu)建三者間的調(diào)控網(wǎng)絡(luò),發(fā)現(xiàn)部分上調(diào)與下調(diào)lncRNA位于調(diào)控網(wǎng)絡(luò)中心位置且結(jié)合較多的miRNA。其中,部分lncRNA可結(jié)合多個(gè)miRNA,如TCONS_00004891、XR_001704571.1和TCONS_00024939分別結(jié)合36、28和27個(gè)miRNA,表明這些lncRNA可能在意蜂工蜂中腸發(fā)育過(guò)程中發(fā)揮重要的調(diào)控功能。此外,也有部分lncRNA共同靶向結(jié)合一個(gè)miRNA,如有多達(dá)50個(gè)lncRNA能夠靶向結(jié)合ame-miR-6001-3p。COLLINS等在研究真社會(huì)昆蟲(chóng)等級(jí)分化過(guò)程中發(fā)現(xiàn),bte-miR-6001-3p在歐洲熊蜂()等級(jí)分化的幼蟲(chóng)中差異表達(dá)顯著,ame-miR-6001-3p在西方蜜蜂等級(jí)分化的幼蟲(chóng)中差異表達(dá)較弱,說(shuō)明miR-6001-3p影響蜜蜂卵巢發(fā)育[48]。上述結(jié)果表明DElncRNA可作為ceRNA吸附結(jié)合miRNA,從而對(duì)意蜂工蜂中腸發(fā)育過(guò)程進(jìn)行調(diào)控。目前,lncRNA的功能研究多集中在人類、哺乳動(dòng)物和重要疾病領(lǐng)域[17],在果蠅中也有少量研究[49],但包括蜜蜂[18-21]在內(nèi)的絕大多數(shù)昆蟲(chóng)的lncRNA研究仍處于初級(jí)階段,相關(guān)的lncRNA信息十分有限。本研究預(yù)測(cè)出的多數(shù)lncRNA,其具體的調(diào)控關(guān)系尚不明確,需要進(jìn)一步克隆出lncRNA全長(zhǎng)序列并深入探究其功能,這是未來(lái)的研究方向。

    本研究?jī)H對(duì)意蜂7和10日齡工蜂中腸進(jìn)行了測(cè)序和lncRNA相關(guān)分析,若要全面解析意蜂工蜂中腸的發(fā)育機(jī)理及調(diào)控機(jī)制,需要對(duì)更多日齡的工蜂中腸進(jìn)行測(cè)序,進(jìn)而在全局水平進(jìn)行更加深入的分析,例如對(duì)所有DElncRNA進(jìn)行基因權(quán)重共表達(dá)分析(WGCNA)。

    4 結(jié)論

    通過(guò)對(duì)意蜂工蜂中腸發(fā)育過(guò)程中的lncRNA及其功能進(jìn)行深入分析,提供了意蜂工蜂中腸發(fā)育過(guò)程中的DElncRNA信息,揭示了DElncRNA廣泛參與意蜂工蜂中腸的新陳代謝、細(xì)胞活動(dòng)和免疫調(diào)控并作為競(jìng)爭(zhēng)性內(nèi)源RNA(ceRNA)發(fā)揮作用,為關(guān)鍵lncRNA的篩選和功能研究提供了必要的數(shù)據(jù)支持。

    [1] PARK D, JUNG J W, CHIO B S, JAYAKODI M, LEE J, LIM J, YU Y, CHOI Y S, LEE M L, PARK Y, CHOI I Y, YANG T Y, EDWARDS O R, NAH G, KWON H W. Uncovering the novel characteristics of Asian honey bee,by whole genome sequencing., 2015, 16(1): 1.

    [2] National Research Council, Division on Earth and Life Studies, Board on Agriculture and Natural Resources, Board on Life Sciences. Committee on the Status of Pollinators in North America.. Washington, D.C: National Academies Press, 2007.

    [3] 周冰峰. 蜜蜂飼養(yǎng)管理學(xué). 廈門: 廈門大學(xué)出版社, 2002.

    ZHOU B F.. Xiamen: Xiamen University Publishing Company, 2002. (in Chinese)

    [4] DJEBALI S, DAVIS C A, MERKEL A, DOBIN A, LASSMANN T, MORTAZAVI A, TANZER A, LAGARDE J, LIN W, SCHLESINGER F, XUE C, MARINOV G K, KHATUN J, WILLIAMS B A, ZALESKI C, ROZOWSKY J, RODER M, KOKOCINSKI F, ABDELHAMID R F, ALIOTO T, ANTOSHECHKIN I, BAER M T, BAR N S, BATUT P, BELL K, BELL I, CHAKRABORTTY S, CHEN X, CHRAST J, CURADO J, DERRIEN T, DRENKOW J, DUMAIS E, DUMAIS J, DUTTAGUPTA R, FALCONNET E, FASTUCAM, FEJES-TOTH K, FERREIRA P, FOISSAC S, FULLWOOD M J, GAO H, GONZALEZ D, GORDON A, GUNAWARDENA H, HOWALD C, JHA S, JOHNSON R, KAPRANOV P, KING B, KINGSWOOD C, LUO O J, PARK E, PERSAUD K, PREALL J B, RIBECA P, RISK B, ROBYR D, SAMMETH M, SCHAFFER L, SEE L H, SHAHAB A, SKANCKE J, SUZUKI A M, TAKAHASHI H, TILGNER H, TROUT D, WALTERS N, WANG H, WROBEL J, YU Y, RUAN X, HAYASHIZAKI Y, HARROW J, GERSTEIN M, HUBBARD T, REYMOND A, ANTONARAKIS S E, HANNON G, GIDDINGS M C, RUAN Y, WOLD B, CARNINCI P, GUIGO R, GINGERAS T R.Landscape of transcription in human cells., 2012, 489(7414): 101-108.

    [5] GORODKIN J, HOFACKER I L. From structure prediction to genomic screens for novel non-coding RNAs., 2011, 7(8): e1002100.

    [6] KANDURI C. Kcnq1ot1: a chromatin regulatory RNA., 2011, 22(4): 343-350.

    [7] TRIPATHI V, SHEN Z, CHAKRABORTY A, GIRI S, FREIER S M, WU X, ZHANG Y, GOROSPE M, PRASANTH S G, LAL A, PRASANTH K V. Long noncoding RNA MALAT1 controls cell cycle progression by regulating the expression of oncogenic transcription factor B-MYB., 2013, 9(3): e1003368.

    [8] LI M, SUN X, CAI H, SUN Y, PLATH M, LI C, LAN X, LEI C, LIN F, BAI Y, CHEN H. Long non-coding RNA ADNCR suppresses adipogenic differentiation by targeting miR-204., 2016, 1859(7): 871-882.

    [9] ULITSKY I, BARTEL D P. LincRNAs: genomics, evolution, and mechanisms., 2013, 154(1): 26-46.

    [10] HUNG T, WANG Y, LIN M F, KOEGEL A K, KOTAKE Y, GRANT G D, HORLINGS H M, SHAH N, UMBRICHT C, WANG P, WANG Y, KONG B, LANEROD A, BORRESEN-DALE A L, KIM S K, VAN D V M, SUKUMAR S, WHITFIELD M L, KELLIS M, XIONG Y, WONG D J, CHANG H Y. Extensive and coordinated transcription of noncoding RNAs within cell cycle promoters., 2011, 43(7): 621-629.

    [11] YOON J H, ABDELMOHSEN K, SRIKANTAN S, YANG X, MARTINDALE J L, DE S, HUARTE M, BECKER K G, GOROSPE M. LincRNA-p21 suppresses target mRNA translation., 2012, 47(4): 648-655.

    [12] CHOONIEDASS-KOTHARI S, EMBERLEY E, HAMEDANI M K, TROUP S, WANG X, CZOSNEK A, HUBE F, MUTAWE M, WATSON P H, LEYGUE E. The steroid receptor RNA activator is the first functional RNA encoding a protein., 2004, 566(1/3): 43-47.

    [13] OGAWA Y, SUN B K, LEE J T. Intersection of the RNA interference and X-inactivation pathways., 2008, 320(5881): 1336-1341.

    [14] KENIRY A, OXLEY D, MONNIER P, KYBA M, DANDOLO L, SMITS G, REIK W. The H19 lincRNA is a developmental reservoir of miR-675 which suppresses growth and Igf1r., 2012, 14(7): 659-665.

    [15] MUDGE J M, HARROW J. Creating reference gene annotation for the mouse C57BL6/J genome assembly., 2015, 26(9/10): 366-378.

    [16] HON C C, RAMILOWSKI J A, HARSHBARGER J, BERTIN N, GOUGH J, DENISENKO E, SCHMEIER S, POULSEN T M, SEVERIN J, LIZIO M, KAWAJI H, KASUKAWA T, ITOH M, BURROUGHS A M, NOMA S, DJEBALI S, ALAM T, MEDVEDEVA Y A, TESTA A C, LIPOVICH L, YIP C W, ABUGESSAISA I, MENDEZ M, HASEGAWA A, TANG D, LASSMANN T, HEUTINK P, BABINA M, WELLS C A, KOJIMA S, NAKAMURA Y, SUZUKI H, DAUB C O, DE-HOON M J, ARNER E, HAYASHIZAKI Y, CARNINCI P, FORREST A R. An atlas of human long non-coding RNAs with accurate 5′ ends., 2017, 543(7644): 199-204.

    [17] 朱斌, 梁沛, 高希武. 長(zhǎng)鏈非編碼RNA (lncRNA)及其在昆蟲(chóng)學(xué)研究中的進(jìn)展. 昆蟲(chóng)學(xué)報(bào), 2016, 59(11): 1272-1281.

    ZHU B, LIANG P, GAO X W. Long noncoding RNAs (lncRNAs) and their research advances in entomology., 2016, 59(11): 1272-1281. (in Chinese)

    [18] 郭昱, 蘇松坤, 陳盛祿, 張少吾, 陳潤(rùn)生. LncRNA在蜜蜂級(jí)型分化中的功能研究. 生物化學(xué)與生物物理進(jìn)展, 2015, 42(8): 750-757.

    GUO Y, SU S K, CHEN C L, ZHANG S W, CHEN R S. The function of lncRNAs in the caste determination of the honeybee., 2015, 42(8): 750-757. (in Chinese)

    [19] HUMANN F C, TIBERIO G J, HARTFELDER K. Sequence and expression characteristics of long noncoding RNAs in honey bee caste development–potential novel regulators for transgressive ovary size., 2013, 8(10): e78915.

    [20] CHEN X, MA C, CHEN C, LU Q, SHI W, LIU Z G, WANG H H, GUO H K. Integration of lncRNA-miRNA-mRNA reveals novel insights into oviposition regulation in honey bees., 2017, 5: e3881.

    [21] JAYAKODI M, JUNG J W, PARK D, AHN Y J, LEE S C, SHIN S Y, CHIN C,YANG T J, KWON H W. Genome-wide characterization of long intergenic non-coding RNAs (lincRNAs) provides new insight into viral diseases in honey beesand., 2015, 16(1): 680.

    [22] BABENDREIER D, JOLLER D, ROMEIS J, BIGLER F, WIDMER F. Bacterial community structures in honeybee intestines and their response to two insecticidal proteins., 2007, 59(3): 600-610.

    [23] KOCH H, ABROL D P, LI J, SCHMID-HEMPEL P. Diversity and evolutionary patterns of bacterial gut associates of corbiculate bees., 2013, 22(7): 2028-2044.

    [24] ELLEGAARD K M, TAMARIT D, JAVELIND E, OLOFSSON T C, ANDERSSON S G, VASQUEZ A. Extensive intra-phylotype diversity in lactobacilli and bifidobacteria from the honeybee gut., 2015, 16(1): 284.

    [25] MOHR K I, TEBBE C C. Diversity and phylotype consistency of bacteria in the guts of three bee species (Apoidea) at an oilseed rape field., 2006, 8(2): 258-272.

    [26] ENGEL P, MARTINSON V G, MORAN N A. Functional diversity within the simple gut microbiota of the honey bee., 2012, 109(27): 11002-11007.

    [27] GILLIAM M, VALENTINE D K. Enterobacteriaceae isolated from foraging worker honey bees,., 1974, 23(1): 38-41.

    [28] LANGMEND B, TRAPNELL C, POP M, SALZBERG S L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome., 2009, 10(3): R25.

    [29] KIM D, PERTEA G, TRAPNELL C, PIMENTEL H, KELLEY R, SALZBERG S L. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions., 2013, 14(4): R36.

    [30] Honeybee Genome Sequencing Consortium.Insights into social insects from the genome of the honeybee., 2006, 443(7114): 931-949.

    [31] KONG L, ZHANG Y, YE Z Q, LIU X Q, ZHAO S Q, WEI L, GAO G. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine., 2007, 35(Web Server issue): 345-349.

    [32] SUN L, LUO H, BU D, ZHAO G, YU K, ZHANG C, LIU Y, CHEN R, ZHAO Y.Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts., 2013, 41(17): e166.

    [33] ROBINSON M D, MCCARTHY D J, SMYTH G K. EdgeR: a bioconductor package for differential expression analysis of digital gene expression data., 2010, 26(1): 139-140.

    [34] REHMSMEIER M, STEFFEN P, HOCHSMANN M, GIEGERICH R. Fast and effective prediction of microRNA/target duplexes., 2004, 10(10): 1507-1517.

    [35] BETEL D, WILSON M, GABOW A, MARKS D S, SANDER C. The microRNA.org resource: targets and expression., 2008, 36(Database issue): 149-153.

    [36] 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.

    [37] 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.

    [38] SALMENA L, POLISENO L, TAY Y KATS L, PANDOLFI P P. Ahypothesis: The Rosetta stone of a hidden RNA language?, 2011, 146(3): 353-358.

    [39] FLIRIAN K, JOSHUA T M. Functional classi?cation and experimental dissection of long noncoding RNAs., 2018, 172(3): 393-407.

    [40] 趙亞周, 田文禮, 胡熠凡, 彭文君. 蜜蜂蜂王漿主蛋白 (MRJPs)的研究進(jìn)展. 應(yīng)用昆蟲(chóng)學(xué)報(bào), 2012, 49(5): 1345-1353.

    ZHAO Y Z, TIAN W L, HU Y F, PENG W J. Research advances in major royal jelly proteins of honeybee., 2012, 49(5): 1345-1353. (in Chinese)

    [41] HALDER G, JOHNSON R L. Hippo signaling: growth control and beyond., 2011, 138(1): 9-22.

    [42] PAN D. The hippo signaling pathway in development and cancer., 2010, 19(4): 491-505.

    [43] CAMARGO F D, GOKHALE S, JOHNNIDIS J B, FU D, BELL G W, JAENISCH R, BRUMMELKAMP T R. YAP1 increases organ size and expands undifferentiated progenitor cells., 2007, 17(23): 2054-2060.

    [44] BITEAU B, HOCHMUTH C E, JASPER H.JNK activity in somatic stem cells causes loss of tissue homeostasis in the aginggut., 2008, 3(4): 442-455.

    [45] ORIHEL T C. The peritrophic membrane: its role as a barrier to infection of the arthropod host//. Academic Press, 1975: 65-73.

    [46] ARONSTEIN K A, MURRAY K D. Chalkbrood disease in honey bees., 2010, 103(Suppl. 1): 20-29.

    [47] 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.

    [48] COLLINS D H, MOHORIANU I, BECKERS M, MOULTON V, DALMAY T, BOURKE A F. MicroRNAs associated with caste determination and differentiation in a primitively eusocial insect., 2017, 7: 45674.

    [49] Li E H, Zhao X, Zhang C, LIU W. Fragile X mental retardation protein participates in non-coding RNA pathways., 2018, 40(2): 87-94.

    (責(zé)任編輯 岳梅)

    Differential expression analysis of long non-coding RNAs during the developmental process ofworker’s Midgut

    GUO Rui1, GENG Sihai1, XIONG Cuiling1, ZHENG Yanzhen1, FU Zhongmin1, WANG Haipeng1, DU Yu1, TONG Xinyu1, ZHAO Hongxia2, CHEN Dafu1

    (1College of Bee Science, Fujian Agriculture and Forestry University, Fuzhou 350002;2Guangdong Institute of Applied Biological Resources, Guangzhou 510260)

    【Objective】Long non-coding RNA (lncRNA) plays an important role in regulation of gene expression, epigenetics and cell cycle in eukaryotes. The objective of this study is to investigate the expression profile and role of lncRNAs in the developmental process ofworker’s midgut. 【Method】In this study, 7- and 10-day-old worker’s midguts of(Am7, Am10) were sequenced using RNA-seq technology and strand-specific library construction method. using Perl script, raw reads were filtered to obtain clean reads with high-quality. Bowtie tool was used to compare clean reads to the ribosome database, and TopHat2 software was employed to compare unmapped clean reads to the reference genome. CPC and CNCI softwares were utilized to predict coding capacity of the transcripts. RT-PCR was performed to identify partial lncRNAs. Investigation of differentially expressed lncRNAs (DElncRNAs) was carried out with edgeR, followed by prediction of upstream and downstream genes, for which GO and KEGG pathway enrichment analyses were performed. RNAhybrid, Miranda and TargetScan softwares were utilized together to predict target miRNAs of DElncRNAs and target genes of miRNAs, and DElncRNAs-miRNAs-mRNAs regulation networks were visualized via Cytoscape. Finally, RT-qPCR was conducted to verify reliability of the sequencing data.【Result】134 802 058 and 147 051 470 raw reads were gained from deep sequencing of Am7 and Am10, respectively, and after stringent filtration, 134 166 157 and 146 293 288 were obtained. In total, 6 353 lncRNAs were predicted, and 3 890 DElncRNAs were obtained based on expression calculation, including 2 005 up-regulated lncRNAs and 1 885 down-regulated lncRNAs. The result of RT-PCR suggested the expected signal bands could be amplified from 8 lncRNAs, implying their true existence. There were 1 793 upstream and downstream genes of DElncRNAs, which were involved in 42 GO terms, including metabolic processes, developmental processes, cellular processes, stress responses, immune system processes and so forth. They were also associated with 251 KEGG pathways, including material metabolism pathways such as carbon metabolism, purine metabolism and fatty acid biosynthesis; energy metabolism pathways such as sulfur metabolism, methane metabolism and oxidative phosphorylation; signaling pathways such as Hippo, Wnt and Notch signaling pathways; cellular immune pathways such as lysosome, endocytosis and ubiquitin mediated proteolysis; humoral immune pathways such as MAPK, Jak-STAT and NF-kappa B pathways, these results demonstrated the DElncRNAs were involved in the material and energy metabolism, cell life activity and immunity regulation in the developmental process ofworker’s midgut. Further analysis showed TCONS_00020918 might play a regulatory part in the nutrient absorption and caste differentiation in the worker’s midgut.Analysis of regulation networks demonstrated that complex networks existed between DElncRNAs and target miRNAs and mRNAs, partial DElncRNAs lie in the central of the networks and link many miRNAs, and partial miRNAs could be bound by many DElncRNAs, which indicated that these DElncRNAs might play an important role during the developmental process of the worker’s midgut. Finally, 5 DElncRNAs were randomly selected for RT-qPCR assay, and the result proved the reliability of sequencing data in this study.【Conclusion】DElncRNA is widely involved in the metabolism, cellular activity and immune regulation ofworker’s midgut, and plays a role as a competitive endogenous RNA (ceRNA). The results provide the necessary data support for the screening and functional study of key lncRNA.

    ; midgut; development; long non-coding RNA; upstream and downstream genes

    10.3864/j.issn.0578-1752.2018.18.016

    2018-03-17;

    2018-05-08

    國(guó)家自然科學(xué)基金(31702190)、國(guó)家現(xiàn)代農(nóng)業(yè)產(chǎn)業(yè)技術(shù)體系建設(shè)專項(xiàng)資金(CARS-44-KXJ7)、福建省教育廳中青年教師教育科研項(xiàng)目(JAT170158)、福建農(nóng)林大學(xué)科技創(chuàng)新專項(xiàng)基金(CXZX2017343)、福建省大學(xué)生創(chuàng)新創(chuàng)業(yè)訓(xùn)練計(jì)劃(201610389053)

    郭睿,E-mail:ruiguo@fafu.edu.cn。耿四海,E-mail:15737313592@163.com。郭睿和耿四海為同等貢獻(xiàn)作者。通信作者陳大福,E-mail:dfchen826@fafu.edu.cn

    猜你喜歡
    意蜂中腸工蜂
    工蜂甲(上)
    工蜂甲(下)
    小保姆成長(zhǎng)記
    如何處理意蜂盜取中蜂群
    蜜蜂雜志(2021年6期)2021-12-05 09:57:44
    勤勞的工蜂
    詳解意蜂盜劫中蜂之過(guò)程
    蜜蜂雜志(2020年6期)2020-12-02 08:07:09
    意蜂蜂蜜和中蜂蜂蜜的區(qū)別
    蜜蜂雜志(2019年3期)2019-12-30 10:25:52
    黃星天牛中腸中內(nèi)切葡聚糖酶的鑒定與酶活性測(cè)定
    杠柳新苷P和E對(duì)粘蟲(chóng)和小地老虎中腸3種解毒酶的影響
    大黑鰓金龜消化與解毒相關(guān)基因的組織表達(dá)研究
    国产男女超爽视频在线观看| 国产激情久久老熟女| 国产日韩一区二区三区精品不卡| 丁香六月欧美| 91精品国产国语对白视频| 黄色片一级片一级黄色片| 久久中文看片网| 丝瓜视频免费看黄片| 亚洲成人免费电影在线观看| 亚洲专区国产一区二区| 美女国产高潮福利片在线看| www.自偷自拍.com| 精品国产一区二区三区四区第35| 亚洲九九香蕉| 亚洲aⅴ乱码一区二区在线播放 | 国产精品成人在线| 精品国产乱码久久久久久男人| 亚洲精品国产精品久久久不卡| 少妇裸体淫交视频免费看高清 | 18禁黄网站禁片午夜丰满| 日本黄色视频三级网站网址 | 日韩视频一区二区在线观看| 午夜福利欧美成人| 免费在线观看视频国产中文字幕亚洲| 精品第一国产精品| 久久精品国产99精品国产亚洲性色 | 一级,二级,三级黄色视频| 亚洲av日韩精品久久久久久密| 久久久久久久国产电影| 国产成人精品久久二区二区免费| 成人18禁高潮啪啪吃奶动态图| 男女高潮啪啪啪动态图| 国产精品一区二区精品视频观看| 一进一出好大好爽视频| 久久九九热精品免费| 国产精品免费大片| 黄色怎么调成土黄色| 91在线观看av| 国产精品国产高清国产av | 亚洲国产欧美日韩在线播放| 在线观看www视频免费| 好男人电影高清在线观看| 亚洲国产欧美日韩在线播放| 村上凉子中文字幕在线| 国产成人精品无人区| 可以免费在线观看a视频的电影网站| 亚洲av片天天在线观看| 亚洲三区欧美一区| 国产在线一区二区三区精| 亚洲专区字幕在线| 国产欧美日韩一区二区精品| 国产高清激情床上av| 人妻 亚洲 视频| 国产精品香港三级国产av潘金莲| 久久中文字幕人妻熟女| bbb黄色大片| 亚洲国产中文字幕在线视频| 人人妻人人添人人爽欧美一区卜| www.精华液| 久久性视频一级片| 男人的好看免费观看在线视频 | 老汉色av国产亚洲站长工具| 人人妻人人澡人人爽人人夜夜| 欧美日韩乱码在线| 久久精品国产清高在天天线| 99re在线观看精品视频| av在线播放免费不卡| 午夜福利视频在线观看免费| 精品人妻在线不人妻| 亚洲av美国av| 十八禁网站免费在线| 91成年电影在线观看| 男女免费视频国产| 一进一出抽搐动态| 欧美日韩亚洲国产一区二区在线观看 | 亚洲专区中文字幕在线| 欧美激情久久久久久爽电影 | 国内毛片毛片毛片毛片毛片| 中亚洲国语对白在线视频| 国产精华一区二区三区| 久久精品人人爽人人爽视色| 午夜免费观看网址| 午夜福利,免费看| 亚洲视频免费观看视频| 亚洲精品国产精品久久久不卡| 久久精品国产亚洲av高清一级| 男人操女人黄网站| 在线观看免费午夜福利视频| 免费人成视频x8x8入口观看| 免费观看人在逋| 亚洲精品中文字幕一二三四区| 亚洲精华国产精华精| 可以免费在线观看a视频的电影网站| 中文字幕高清在线视频| 久久国产精品影院| 男人舔女人的私密视频| 下体分泌物呈黄色| 成年版毛片免费区| 午夜日韩欧美国产| netflix在线观看网站| 亚洲专区字幕在线| 久久久久久久国产电影| 另类亚洲欧美激情| 亚洲av美国av| 啦啦啦免费观看视频1| 国产又色又爽无遮挡免费看| 欧美黑人精品巨大| 宅男免费午夜| 久久中文字幕一级| 一进一出抽搐动态| 精品第一国产精品| 女性被躁到高潮视频| 又黄又粗又硬又大视频| 久久人人爽av亚洲精品天堂| 国产一区二区三区综合在线观看| 久久精品人人爽人人爽视色| 亚洲五月天丁香| 国产欧美日韩一区二区精品| 一级片'在线观看视频| 韩国av一区二区三区四区| 久久久久久久国产电影| 亚洲,欧美精品.| 亚洲欧美一区二区三区黑人| 中亚洲国语对白在线视频| 久久99一区二区三区| 国产男女超爽视频在线观看| 中国美女看黄片| 精品国产一区二区久久| 免费久久久久久久精品成人欧美视频| 老司机亚洲免费影院| 人妻 亚洲 视频| 欧美乱妇无乱码| 啦啦啦视频在线资源免费观看| 国产99白浆流出| 久久久精品免费免费高清| 19禁男女啪啪无遮挡网站| 91麻豆精品激情在线观看国产 | 亚洲专区中文字幕在线| 老熟女久久久| 国产91精品成人一区二区三区| 色精品久久人妻99蜜桃| 色婷婷av一区二区三区视频| 首页视频小说图片口味搜索| 欧美日韩一级在线毛片| 男女之事视频高清在线观看| 久久久国产精品麻豆| 狠狠狠狠99中文字幕| 国产熟女午夜一区二区三区| 夫妻午夜视频| 亚洲一区二区三区不卡视频| 成人永久免费在线观看视频| 免费观看精品视频网站| av国产精品久久久久影院| 热re99久久国产66热| 老司机亚洲免费影院| 搡老岳熟女国产| 最近最新中文字幕大全免费视频| 亚洲成a人片在线一区二区| 99精品在免费线老司机午夜| 欧美国产精品va在线观看不卡| 精品国产乱子伦一区二区三区| 国精品久久久久久国模美| 欧美日韩av久久| 精品人妻1区二区| 十八禁网站免费在线| 久久精品91无色码中文字幕| 精品国产国语对白av| avwww免费| tocl精华| 国产精品乱码一区二三区的特点 | svipshipincom国产片| 国产成+人综合+亚洲专区| 午夜视频精品福利| 男女午夜视频在线观看| 欧美在线黄色| 免费看a级黄色片| 99re在线观看精品视频| 欧美激情极品国产一区二区三区| 亚洲av欧美aⅴ国产| 女人精品久久久久毛片| 成人特级黄色片久久久久久久| 免费在线观看完整版高清| cao死你这个sao货| 久久午夜亚洲精品久久| 日韩有码中文字幕| 人妻久久中文字幕网| 久久ye,这里只有精品| 亚洲中文日韩欧美视频| 在线观看www视频免费| 18禁国产床啪视频网站| 女人被狂操c到高潮| 国产亚洲av高清不卡| 欧美激情极品国产一区二区三区| 777米奇影视久久| 在线观看免费视频日本深夜| 免费观看人在逋| 亚洲国产中文字幕在线视频| 女性被躁到高潮视频| 精品人妻1区二区| 国产精品国产av在线观看| 午夜福利,免费看| 精品久久久久久久毛片微露脸| 在线看a的网站| 99久久国产精品久久久| 欧美久久黑人一区二区| 建设人人有责人人尽责人人享有的| 美女高潮到喷水免费观看| 成年女人毛片免费观看观看9 | 在线天堂中文资源库| 国产熟女午夜一区二区三区| 夜夜躁狠狠躁天天躁| 欧美人与性动交α欧美软件| 亚洲午夜精品一区,二区,三区| av线在线观看网站| 两个人免费观看高清视频| videosex国产| 日韩有码中文字幕| 看黄色毛片网站| 村上凉子中文字幕在线| 亚洲欧美一区二区三区黑人| 国产亚洲欧美98| 黄色怎么调成土黄色| 中文字幕高清在线视频| a级毛片黄视频| 日韩免费av在线播放| 深夜精品福利| 青草久久国产| 精品少妇久久久久久888优播| 国产在线观看jvid| 最近最新中文字幕大全免费视频| 男人舔女人的私密视频| 久久久久久久午夜电影 | 亚洲第一青青草原| 国产欧美亚洲国产| 十八禁人妻一区二区| 女人被狂操c到高潮| 丝袜美足系列| 少妇粗大呻吟视频| 久久久国产精品麻豆| 日韩欧美在线二视频 | 亚洲成国产人片在线观看| 女性被躁到高潮视频| 黄色视频不卡| 91老司机精品| 两个人免费观看高清视频| 18禁黄网站禁片午夜丰满| 少妇 在线观看| 大陆偷拍与自拍| 自拍欧美九色日韩亚洲蝌蚪91| 99国产精品一区二区三区| 国产真人三级小视频在线观看| 国产亚洲精品久久久久5区| 国产精华一区二区三区| 亚洲综合色网址| 一级毛片高清免费大全| 欧美日韩亚洲高清精品| 国产精品一区二区在线观看99| 欧美乱妇无乱码| 精品国内亚洲2022精品成人 | 精品人妻在线不人妻| 视频在线观看一区二区三区| 两性午夜刺激爽爽歪歪视频在线观看 | 久久精品人人爽人人爽视色| 黄色怎么调成土黄色| 日本五十路高清| 久久精品国产亚洲av高清一级| 丰满饥渴人妻一区二区三| 久久国产精品影院| 亚洲欧美色中文字幕在线| 俄罗斯特黄特色一大片| 久久久精品国产亚洲av高清涩受| 中文亚洲av片在线观看爽 | 国产精品99久久99久久久不卡| 国产精品自产拍在线观看55亚洲 | 在线观看一区二区三区激情| 99热只有精品国产| 国内毛片毛片毛片毛片毛片| 日本wwww免费看| 亚洲精品av麻豆狂野| 亚洲av成人不卡在线观看播放网| 叶爱在线成人免费视频播放| 真人做人爱边吃奶动态| 日本黄色日本黄色录像| 两个人看的免费小视频| 国产又爽黄色视频| 久久天躁狠狠躁夜夜2o2o| 精品国产乱子伦一区二区三区| 成年版毛片免费区| 99re在线观看精品视频| 日韩视频一区二区在线观看| 高清毛片免费观看视频网站 | 亚洲精品国产精品久久久不卡| 国产精华一区二区三区| 亚洲精品国产色婷婷电影| 精品一区二区三区视频在线观看免费 | 国产一区在线观看成人免费| 免费人成视频x8x8入口观看| 欧美中文综合在线视频| 欧美 亚洲 国产 日韩一| 日韩制服丝袜自拍偷拍| 人妻 亚洲 视频| 国产精品久久久久久人妻精品电影| 亚洲第一欧美日韩一区二区三区| 欧美性长视频在线观看| 欧美在线黄色| 这个男人来自地球电影免费观看| 欧美 亚洲 国产 日韩一| 90打野战视频偷拍视频| 窝窝影院91人妻| 欧美老熟妇乱子伦牲交| 亚洲欧洲精品一区二区精品久久久| 国产成人系列免费观看| e午夜精品久久久久久久| 18禁观看日本| 久久久国产成人免费| 欧美成人午夜精品| 国产一区二区三区综合在线观看| aaaaa片日本免费| 91九色精品人成在线观看| 在线观看www视频免费| 亚洲 国产 在线| 久久天堂一区二区三区四区| 18在线观看网站| 两个人看的免费小视频| 老司机靠b影院| 俄罗斯特黄特色一大片| 亚洲精品久久午夜乱码| 老汉色av国产亚洲站长工具| 亚洲 欧美一区二区三区| 男女床上黄色一级片免费看| 自拍欧美九色日韩亚洲蝌蚪91| 久久香蕉精品热| 91精品三级在线观看| 18禁黄网站禁片午夜丰满| 精品国内亚洲2022精品成人 | 一进一出抽搐gif免费好疼 | a级片在线免费高清观看视频| 成人手机av| 精品卡一卡二卡四卡免费| 黄色视频不卡| 成人精品一区二区免费| 国产成人精品在线电影| 色播在线永久视频| 国产精品自产拍在线观看55亚洲 | 欧美在线一区亚洲| 老司机在亚洲福利影院| 最近最新免费中文字幕在线| 一边摸一边做爽爽视频免费| 久久性视频一级片| 曰老女人黄片| 黄色视频不卡| 99久久精品国产亚洲精品| 中文字幕制服av| 黄色视频,在线免费观看| 久久ye,这里只有精品| 999久久久国产精品视频| 久久久国产精品麻豆| 日本vs欧美在线观看视频| 欧美精品啪啪一区二区三区| 久久久精品免费免费高清| 中文欧美无线码| 91在线观看av| 亚洲第一青青草原| 大片电影免费在线观看免费| av欧美777| √禁漫天堂资源中文www| 国产高清激情床上av| 法律面前人人平等表现在哪些方面| 波多野结衣av一区二区av| 女人精品久久久久毛片| 日本欧美视频一区| 国产精品 欧美亚洲| av欧美777| 人妻久久中文字幕网| 伊人久久大香线蕉亚洲五| 丝袜人妻中文字幕| 在线国产一区二区在线| 国产不卡一卡二| 母亲3免费完整高清在线观看| 亚洲综合色网址| 国产精品1区2区在线观看. | 99国产精品一区二区三区| 色老头精品视频在线观看| 深夜精品福利| 中文字幕人妻丝袜制服| svipshipincom国产片| 不卡一级毛片| 欧美老熟妇乱子伦牲交| 热99re8久久精品国产| 日韩免费av在线播放| 国产麻豆69| 日本黄色视频三级网站网址 | 国产国语露脸激情在线看| 91九色精品人成在线观看| 精品一区二区三卡| 精品乱码久久久久久99久播| 满18在线观看网站| 1024视频免费在线观看| 色尼玛亚洲综合影院| 日韩中文字幕欧美一区二区| 亚洲视频免费观看视频| 日韩大码丰满熟妇| 亚洲午夜精品一区,二区,三区| 久久久久久久久久久久大奶| 曰老女人黄片| 叶爱在线成人免费视频播放| 成年版毛片免费区| 成在线人永久免费视频| 91精品国产国语对白视频| 亚洲第一av免费看| 久久精品91无色码中文字幕| 男女高潮啪啪啪动态图| 看片在线看免费视频| 亚洲精品在线美女| 狠狠婷婷综合久久久久久88av| 久久人妻av系列| 国产欧美日韩一区二区三区在线| 国产男女超爽视频在线观看| 精品欧美一区二区三区在线| 国产精品乱码一区二三区的特点 | 婷婷精品国产亚洲av在线 | 欧美日韩视频精品一区| 黄色a级毛片大全视频| 久久久久精品人妻al黑| 中文字幕色久视频| 国产亚洲精品第一综合不卡| 亚洲五月天丁香| 在线十欧美十亚洲十日本专区| 久久精品熟女亚洲av麻豆精品| 久久精品国产99精品国产亚洲性色 | www.999成人在线观看| 国产99久久九九免费精品| 免费在线观看完整版高清| svipshipincom国产片| 人人妻人人爽人人添夜夜欢视频| 国产aⅴ精品一区二区三区波| 欧美黑人精品巨大| 免费在线观看黄色视频的| 丝袜美足系列| 色在线成人网| 少妇裸体淫交视频免费看高清 | 极品教师在线免费播放| 黄色片一级片一级黄色片| 欧美人与性动交α欧美软件| 国产精品av久久久久免费| 久久国产精品大桥未久av| 中文字幕精品免费在线观看视频| 国产亚洲精品久久久久5区| 看黄色毛片网站| 亚洲熟妇熟女久久| 国产精品电影一区二区三区 | 18禁裸乳无遮挡动漫免费视频| 亚洲中文字幕日韩| 在线看a的网站| 亚洲情色 制服丝袜| 亚洲伊人色综图| 国产精品自产拍在线观看55亚洲 | 夜夜躁狠狠躁天天躁| 丝袜在线中文字幕| 欧美午夜高清在线| 成人黄色视频免费在线看| 久久香蕉精品热| 亚洲av第一区精品v没综合| 看免费av毛片| 国产91精品成人一区二区三区| 夜夜爽天天搞| 成人黄色视频免费在线看| 午夜日韩欧美国产| 国产xxxxx性猛交| 99国产精品一区二区蜜桃av | 999久久久精品免费观看国产| 亚洲av欧美aⅴ国产| 美女福利国产在线| 久久精品aⅴ一区二区三区四区| 久久性视频一级片| 成人三级做爰电影| 大型黄色视频在线免费观看| 国产一区二区激情短视频| 黄色片一级片一级黄色片| 搡老熟女国产l中国老女人| 精品人妻1区二区| 亚洲国产欧美一区二区综合| 怎么达到女性高潮| 大香蕉久久网| 老熟妇乱子伦视频在线观看| 亚洲精华国产精华精| 国产精品1区2区在线观看. | 80岁老熟妇乱子伦牲交| 狠狠狠狠99中文字幕| av天堂久久9| 大型黄色视频在线免费观看| 成人亚洲精品一区在线观看| 50天的宝宝边吃奶边哭怎么回事| 村上凉子中文字幕在线| 一级作爱视频免费观看| 手机成人av网站| 中文字幕人妻丝袜制服| 午夜福利影视在线免费观看| 国产一区二区激情短视频| 国产人伦9x9x在线观看| 12—13女人毛片做爰片一| 亚洲av成人不卡在线观看播放网| 久久久久久久精品吃奶| 免费av中文字幕在线| www日本在线高清视频| 精品乱码久久久久久99久播| 大片电影免费在线观看免费| 91字幕亚洲| 丝袜人妻中文字幕| 人妻丰满熟妇av一区二区三区 | 在线观看免费高清a一片| 国产又色又爽无遮挡免费看| 男人操女人黄网站| av不卡在线播放| 午夜福利,免费看| 亚洲成国产人片在线观看| 搡老岳熟女国产| 人人妻人人澡人人看| 国产亚洲一区二区精品| 欧美激情高清一区二区三区| 国产精品免费大片| 好男人电影高清在线观看| 18禁美女被吸乳视频| 免费在线观看日本一区| 午夜福利乱码中文字幕| 国产国语露脸激情在线看| 色老头精品视频在线观看| 99精国产麻豆久久婷婷| 搡老熟女国产l中国老女人| 狠狠狠狠99中文字幕| 美女视频免费永久观看网站| 欧美日韩一级在线毛片| 丝袜在线中文字幕| 18在线观看网站| 国产视频一区二区在线看| 精品国产乱码久久久久久男人| 亚洲情色 制服丝袜| 久久精品国产亚洲av香蕉五月 | 91精品国产国语对白视频| 新久久久久国产一级毛片| 国产成人av教育| 老司机在亚洲福利影院| 91成年电影在线观看| 精品国内亚洲2022精品成人 | 亚洲性夜色夜夜综合| 下体分泌物呈黄色| 人成视频在线观看免费观看| 老熟妇乱子伦视频在线观看| 成人手机av| 亚洲全国av大片| 亚洲人成电影免费在线| 免费黄频网站在线观看国产| 高清在线国产一区| 一二三四社区在线视频社区8| 国产精品电影一区二区三区 | 精品午夜福利视频在线观看一区| 黑人巨大精品欧美一区二区蜜桃| 久久这里只有精品19| 亚洲中文字幕日韩| 欧美国产精品一级二级三级| 女同久久另类99精品国产91| 大陆偷拍与自拍| 9色porny在线观看| 深夜精品福利| 免费观看精品视频网站| 99在线人妻在线中文字幕 | 国产高清激情床上av| 性少妇av在线| 老熟女久久久| 久久性视频一级片| 成年人午夜在线观看视频| 一个人免费在线观看的高清视频| 色尼玛亚洲综合影院| 村上凉子中文字幕在线| 成人手机av| 成年版毛片免费区| 最近最新中文字幕大全免费视频| www日本在线高清视频| 久久久国产精品麻豆| 脱女人内裤的视频| 亚洲熟妇中文字幕五十中出 | 久久久久久久精品吃奶| 国产成人精品久久二区二区免费| 丝袜人妻中文字幕| 欧美日韩亚洲高清精品| a在线观看视频网站| 国产高清视频在线播放一区| bbb黄色大片| 男女之事视频高清在线观看| 多毛熟女@视频| 国产欧美日韩一区二区三| 午夜老司机福利片| 中文亚洲av片在线观看爽 | 色综合欧美亚洲国产小说| 精品电影一区二区在线| 精品一区二区三区四区五区乱码| 国产一区二区三区视频了| 丝袜人妻中文字幕| 12—13女人毛片做爰片一| 色尼玛亚洲综合影院| 欧美日本中文国产一区发布| 一级黄色大片毛片| 久99久视频精品免费| 久久天堂一区二区三区四区| 丁香六月欧美| 大香蕉久久成人网| 成人国产一区最新在线观看| 9色porny在线观看| 18禁黄网站禁片午夜丰满| 国产精品香港三级国产av潘金莲| 亚洲欧美一区二区三区久久| 一区在线观看完整版| 精品久久久久久,| 91老司机精品| 亚洲久久久国产精品|