999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

谷子萌發期響應干旱脅迫的基因表達譜分析

2018-05-14 09:26:46許冰霞尹美強溫銀元裴帥帥柯貞進張彬原向陽
中國農業科學 2018年8期
關鍵詞:差異分析

許冰霞,尹美強,溫銀元,裴帥帥,柯貞進,張彬,原向陽

?

谷子萌發期響應干旱脅迫的基因表達譜分析

許冰霞,尹美強,溫銀元,裴帥帥,柯貞進,張彬,原向陽

(山西農業大學農學院,山西太谷 030801)

【目的】谷子是一種耐旱作物,通過二代測序技術獲得大量的谷子萌發期響應干旱脅迫的差異基因,進而挖掘谷子萌發期抵御干旱的關鍵基因及其相關的分子機制?!痉椒ā恳詴x谷45為材料,谷子萌發期分別用18%PEG-6000干旱脅迫(處理組)和蒸餾水(對照組)處理種子并測定1、10和18 h種子的SOD、POD和CAT活性。SOD活性用氮藍四唑(NBT)法測定,POD活性用愈創木酚法測定,CAT活性用比色法測定;對萌發10和18 h種子的對照組和處理組構建cDNA文庫并進行差異基因表達譜分析;利用Bowtie將reads比對到參考基因組,采用RSEM對bowtie的比對結果進行表達量估計;使用DESeq進行差異表達基因分析;利用NR、Swiss-Prot、KEGG、COG和GO在線數據庫對差異基因進行功能注釋,挖掘調控谷子萌發的關鍵基因;利用qRT-PCR驗證測序結果的可靠性?!窘Y果】處理組的SOD活性整體比對照組高,而POD活性和CAT活性與之相反;隨著萌發時間的變化,SOD活性在不斷地增加,但CAT和POD活性逐漸減小。基因表達譜序列與所選參考基因組序列高度一致,基因表達呈現出高度不均一性。通過高通量測序最后獲得35 470個基因,以RPKM≥0.01為篩選標準,對照樣本中分別篩選出24 030和24 486個表達基因,PEG干旱脅迫處理10和18 h的樣本分別篩選出24 019和23 877個表達基因;差異表達基因分析表明,谷子萌發10和18 h分別篩選出456和545個差異基因,其中87和267個上調表達基因,369和278個下調表達基因;GO功能顯著性富集分析表明,差異基因主要涉及代謝過程,細胞進程和響應刺激;KEGG富集分析表明,差異基因參與到苯丙烷代謝和植物激素信號轉導過程;通過qRT-PCR對5個差異基因在干旱脅迫下種子萌發時的表達分析表明,其表達趨勢與表達譜分析結果基本一致。【結論】差異表達基因廣泛涉及到糖、蛋白質、核酸等生物大分子代謝、次生代謝和能量代謝等過程;和可能在干旱脅迫下調節種子的萌發。

谷子;干旱脅迫;種子萌發;基因表達譜

0 引言

【研究意義】種子萌發率在作物生產中直接影響作物的出苗率以及作物生長和最終的產量。種子自身的休眠與萌發機制調節種子的萌發率,而外界的環境條件也會影響種子的萌發率。其中,水分是種子萌發的首要條件,種子萌發起始需要吸收大量水分啟動一系列的生理生化活動為種子發芽做準備。提高干旱環境下作物種子的萌發率,可以解決干旱環境下出苗率低、出苗不齊等問題,從而達到保障產量穩定的目的?!厩叭搜芯窟M展】谷子是一種古老的禾本科耐旱穩產作物,與象草、柳枝稷等一些生物燃料作物有很近的親緣關系,同時也是世界上許多干旱和半干旱地區的重要糧食和飼料作物[1],因此,已有許多關于谷子抗旱性的鑒定,以及形態指標和生理指標與抗旱性關系的研究報道[2-3],如裴帥帥等[4]和代小冬等[5]采用PEG滲透劑模擬干旱脅迫對不同品種谷子種子的抗旱性進行研究,并通過測定可溶性蛋白、脯氨酸、SOD和POD活性、MDA含量等生理指標以及相對發芽勢、相對發芽率、發芽指數和抗旱指數等指標來鑒定不同品種谷子的抗旱性。同時谷子也是進化研究、生理研究和基因組研究的模型[6],隨著栽培谷子和野生谷子基因組測序的完成及公布[7-8],使得基因組數據和基因資源的積累越來越多;因此有關谷子的抗旱性研究也在逐漸增加,許多研究是對谷子幼苗進行干旱脅迫,例如:谷子的鑒定[9]及脫水脅迫下差異表達基因的轉錄組分析[10]。盡管在谷子中已克隆出DnaJ蛋白(SiDnaJ)[11]、DREB轉錄因子[12]、C2H2型鋅指蛋白基因()[13]等抗旱基因,但更多的研究是將抗旱基因通過基因工程的手段轉入其他模式植物鑒定其抗逆性,如竇祎凝等[14]通過谷子NAC類轉錄因子基因在擬南芥中過表達分析,發現可能通過ABA信號途徑、氧化脅迫調控等途徑正向調控植物在干旱條件下的萌發過程;王一帆等[15]發現在NaCl、PEG、ABA、GA、MeJA等不同處理下都會響應,尤其在鹽脅迫條件下更為明顯,并將其轉化水稻驗證,進而推測該基因會對谷子抗鹽和抗逆起到一定作用?!颈狙芯壳腥朦c】目前,關于谷子的抗旱性研究較多集中于抗旱性鑒定和幼苗期,但是在分子水平對谷子種子萌發期的抗旱機制研究鮮見報道?!緮M解決的關鍵問題】本研究利用數字表達譜來探究萌發期抗旱性強的晉谷45在干旱脅迫10和18 h時基因表達種類及豐度,快速準確地挖掘谷子萌發期的特有抗旱基因,在分子水平深入研究谷子萌發期對干旱脅迫響應的分子機制,并為轉基因耐旱谷子提供參考和提高谷子萌發率提供理論依據。

1 材料與方法

1.1 試驗材料

以晉谷45(山西農業大學植物生理實驗室抗旱萌發鑒定中抗性最強的品種)為材料,挑選成熟飽滿大小均一且谷殼完好無病蟲害的種子,用0.1% HgCl2浸泡消毒10 min,蒸餾水沖洗5—6次,晾干,然后均勻擺放在培養皿中,每皿100粒。試驗組用18%(m/v)PEG-6000溶液進行模擬干旱脅迫處理,對照組使用蒸餾水處理,兩組均進行3次重復,將處理好的樣品置于25℃恒溫培養箱中暗培養1、10和18 h,用于SOD、POD和CAT活性的測定,暗培養10和18 h,液氮速凍-80℃保存備用。

1.2 SOD、POD和CAT活性的測定

SOD活性用氮藍四唑(NBT)法測定,POD活性用愈創木酚法測定,CAT活性用比色法測定。

1.3 總RNA的提取、cDNA文庫的構建和測序

總RNA的提取、檢測及cDNA文庫的構建由中國北京百邁克公司完成。

1.3.1 總RNA的檢測 分別采用Nanodrop、Qubit 2.0、Agilent 2100方法檢測RNA樣品的純度、濃度和完整性等,合格的RNA進行基因表達譜測序。

1.3.2 cDNA文庫構建 主要流程如下:

(1)用帶有Oligo(dT)的磁珠富集真核生物mRNA;

(2)加入fragmentation buffer將mRNA進行隨機打斷;

(3)以mRNA為模板,用六堿基隨機引物(random hexamers)合成第一條cDNA鏈,然后加入緩沖液、dNTPs、RNase H和DNA polymeraseⅠ合成第二條cDNA鏈,利用AMPure XP beads純化cDNA;

(4)純化的雙鏈cDNA再進行末端修復、加A尾并連接測序接頭,然后用AMPure XP beads進行片段大小選擇;

(5)最后通過PCR富集得到cDNA文庫。

共4個樣本進行高通量測序建庫,谷子萌發10 h對照組命名為10hCK,谷子萌發10hPEG脅迫處理組命名為10hPEG,谷子萌發18 h對照組命名為18hCK,谷子萌發18hPEG脅迫處理組命名為18hPEG。測序數據為3次重復。

1.3.3 測序 用HiSeq2500進行高通量測序,測序讀長為SE50。

1.4 基因表達分析及功能注釋

1.4.1 基因表達分析 采用Bowtie[16]將樣品的reads與參考基因組豫谷1號的基因組進行序列比對,使用RSEM[17]對bowtie的比對結果進行表達量估計,最后利用RPKM[16]值來反應基因的表達豐度。

1.4.2 差異基因表達分析 使用DESeq[18]進行差異表達分析,篩選出差異表達基因;最后采用Benjamini- Hochberg方法對基因的表達值進行統計檢驗,也就是用FDR(false discovery rate)作為差異表達基因篩選的指標。在篩選差異基因的過程中,將FDR<0.01且FC(fold change)≥2作為篩選標準。

1.4.3 功能注釋 使用BLAST[18]軟件進行搜索比對,找到每一條基因的描述。在NR[19]、Swiss-Prot[20]、KEGG[21]、COG[22]和GO[23]在線數據庫中進行比對,獲得差異表達基因的功能注釋。利用topGO[24]軟件對差異基因進行富集分析。

1.5 實時熒光定量PCR分析

隨機挑選5個響應干旱的差異表達基因,用實時熒光定量技術驗證表達譜測序數據的可靠性。通過Trizol法提取樣品總RNA,用反轉錄試劑盒PrimeScript?RT Reagent Kit With gDNA Eraser(TaKaRa)合成cDNA。熒光定量使用SYBR?Premix Ex Taq Ⅱ(Tli RNaseH Plus)試劑盒(TaKaRa)。使用Beacon Designer7.5設計引物(表1),長度為18—20 bp,擴增長度為80—200 bp,以谷子組成型表達基因作為內參基因,引物由北京六合華大基因科技股份有限公司合成。擴增程序為95℃30 s;95℃5 s,60℃34 s;40個循環。采用2-△△CT計算差異表達基因的相對表達量,每個樣品3個技術重復。

2 結果

2.1 谷子萌發過程中SOD、CAT和POD活性的差異分析

整體來看,種子萌發過程中CAT活性無顯著性差異,SOD活性有極顯著差異(<0.01),POD活性有顯著差異(<0.05)。處理組的SOD活性整體比對照組高,而POD活性和CAT活性與之相反;處理組比對照組的SOD活性分別高出20.47%、26.18%和21.06%,POD活性分別低3.02%、24.86%和12.81%,而CAT活性分別降低20.00%、12.90%和2.82%,隨著萌發時間的變化SOD活性在不斷地增加,但CAT和POD活性逐漸減?。▓D1)。

2.2 測序質量的分析

對Raw Data中的測序引物、接頭等人工序列進行截除,并將其中的低質量數據過濾,獲得高質量數據即Clean data(表2)。谷子萌發10 h的對照樣本和PEG干旱脅迫處理樣本的Total reads分別為其12 080 922和10 669 558,谷子萌發18 h的對照樣本和PEG干旱脅迫處理樣本的Total reads分別為12 248 194和10 405 536,reads長度為50—100 bp,表明數據深度足夠進行基因表達定量分析。GC含量在4組樣本中均大于50%,分別為56.38%、56.59%、56.72%和56.59%,同時堿基質量值(quality score或Q-score)Q30都在85%以上,表明測序質量非??煽浚梢赃M一步進行分析。

表1 引物序列信息

不同字母表示在對照組和處理組之間的差異顯著

表2 樣品序列與參考基因組比對分析

對照樣本和PEG干旱脅迫處理樣本的reads被注釋到參考基因組的情況相似,都大于86%,表明所選用的谷子材料與參考基因組的生物學分類關系較近,遺傳基礎差異較小,參考基因組滿足后續分析的要求。另一方面,比對位置多于1個的reads數量較少(對照樣本分別為14.67%和15.66%,PEG干旱脅迫處理樣本分別為15.57%和16.75%),大多數reads比對到1個位置(對照樣本分別為85.33%和84.34%,PEG干旱脅迫處理樣本分別為84.43%和83.25%),表明reads質量較好,相對比對率高。片段隨機性分析表明mRNA片段隨機性較高,樣本沒有嚴重的降解現象;測序數據飽和度分析(隨著數據量的增加,新增加的基因逐漸減少)表明樣本有效數據充足;這兩者從不同的角度說明了表達譜測序文庫質量較高可以用于后續的信息分析。

2.3 谷子萌發時期響應干旱脅迫的基因表達分析

通過高通量測序最后獲得35 470個基因,這些基因的表達豐度用RPKM(Reads Per Kilobase per Million mapped reads)來反映。以RPKM≥0.01為篩選標準,對照樣本分別共篩選出24 030和24 486個表達基因,PEG干旱脅迫處理樣本分別篩選出24 019和23 877個表達基因,4組樣本中落在RPKM的不同區間的基因數目及比例大致相同,但是同一樣本中落在RPKM的不同區間的的基因數目及比例存在顯著差異,RPKM在0—5的基因達到35%以上,接近50%的基因的RPKM值在6—80(表3)。

表3 不同表達水平區間的基因數目

2.4 谷子萌發時期響應干旱脅迫的基因差異表達分析

以FDR<0.01且差異倍數FC(Fold Change)≥2篩選標準,對對照樣本和PEG處理樣本進行差異表達基因分析,發現在谷子萌發10 h和18 h分別獲得了456和545個差異表達基因,分別有87和267個基因上調表達,369和278個基因下調表達。其中,表達上(下)調2—4倍的基因分別占差異表達上(下)調基因總數的15.35%和38.16%(66.89%和39.63%);表達上(下)調4—8倍的基因分別占2.85%和7.7%(10.75%和8.62%);表達上(下)調8—16倍的基因分別占0.44%和2.20%(2.41%和2.01%);表達上(下)調16倍以上的基因分別占0.44%和0.73%(0.88%和0.92%)。從火山圖中可以直觀地看出對照樣本與PEG處理樣本之間的表達水平具有顯著的差異,同時也可以看出谷子萌發10 h和18 h時差異表達基因數目整體有增加,說明隨著脅迫時間的增加,干旱對谷子萌發的影響也在不斷增加,從而阻礙谷子萌發的進程(圖2)。由圖3和圖4可知,種子萌發10 h和18 h對照組有20個差異基因同樣在對照組和處理組之間也是差異表達,這20個基因在谷子萌發階段響應干旱脅迫,同一基因在處理組與對照組之間表達量有很大變化,在對照組之間隨著種子的萌發表達量也有很大變化(電子附表1)。

2.5 差異表達基因功能的注釋和富集分析

使用BLAST軟件分別將谷子萌發10和18 h 的456(電子附表2)和545(電子附表3)個差異基因與nr(non-redundant,NCBI)、Swiss-Prot、COG(Clusters of Orthologous Groups of proteins,NCBI)、GO(the Gene Ontology)和KEGG(Kyoto Encyclopedia of Genes and Genomes)5個蛋白質數據庫進行序列比對,獲得差異表達基因的注釋信息。

將比對到COG、GO和KEGG的差異表達基因再進行分類分析。這些差異表達基因注釋到GO數據庫中進行功能分類。結果顯示,谷子萌發10 h和18 h注釋到細胞組分(cellular component)、分子功能(molecular function)和生物過程(biological process)的差異基因分別有326和367、308和325及304和348個。以-value≤0.01為篩選標準,發現生物過程(圖5)中有128和140個基因顯著富集到20和15功能類別,其中氧化應激反應(response to oxidative stress,GO:0006979)、氧化還原過程(oxidation- reduction process,GO:0055114)、纖維素微纖絲組織(cellulose microfibril organization,GO:0010215)和有性生殖(sexual reproduction,GO:0019953)是種子萌發10 h和18 h共同富集的通路。分子功能(圖6)包含17個功能類別,分別有91和149個基因富集到這些功能類別中,其中鐵離子結合(iron ion binding,GO:0005506)、酪氨酸氨解酶活性(tyrosine ammonia- lyase activity,GO:0052883)、木葡聚糖轉移酶活性(xyloglucan:xyloglucosyl transferase activity,GO:0016762)和纖維二糖葡萄糖苷酶活性(cellobiose glucosidase activity,GO:0080079)是種子萌發10 h和18 h共同富集的通路。在這些共同富集的生物功能中氧化應激反應,氧化還原過程和鐵離子結合是有顯著差異的,隨著萌發時間的增加富集到這些功能類別的差異基因也在增加,再加上一個基因會有一個或者多個功能注釋,因此,可以推斷谷子萌發期響應干旱脅迫時,其代謝是復雜的。

A:10hCKvs10hPEG;B:18hCKvs18hPEG

圖3 差異基因維恩圖

COG分析發現在種子萌發10和18 h時差異基因涉及21個功能類別(圖7),在這些基因中基因數目最多的功能類別是一般功能(general function prediction only),基因數目達到52個;在種子萌發10 h時基因數目大于10的功能類別有轉錄(transcription)有18個基因,復制、重組和修復(replication, recombination and repair)有16個基因,信號轉導機制(signal transduction mechanisms)有26個基因,能量的產生和轉化(energy production and conversion)有19個基因,碳水化合物轉運和代謝(carbohydrate transport and metabolism)有27個基因以及氨基酸轉運與代謝(amino acid transport and metabolism)有24個基因,其他的功能類別涉及到的基因均等于或低于10;而在種子萌發18 h時基因數目大于10的功能類別分別為翻譯、核糖體結構與生物合成(translation, ribosomal structure and biogenesis)有15個基因,轉錄(transcription)有21個基因,復制、重組和修復(replication, recombination and repair)有16個基因,信號轉導機制(signal transduction mechanisms)有22個基因,胞壁/膜生物發生(cell wall/membrane/envelope biogenesis)有14個基因,蛋白質翻譯后修飾與轉運(posttranslational modification, protein turnover, chaperones)有15個基因,碳水化合物轉運和代謝(carbohydrate transport and metabolism)有26個基因,能量的產生和轉化(energy production and conversion)有16個基因,氨基酸轉運與代謝(amino acid transport and metabolism)為19個基因,次生代謝產物的生物合成、運輸和分解代謝(secondary metabolites biosynthesis, transport and catabolism)有15個基因。盡管其他的功能類別的差異基因數目并不客觀,但是囊括種類較多說明種子萌發階段受到PEG脅迫后,物質代謝及能量轉化受到較大的影響并且一直伴隨著種子萌發整個過程,而種子萌發到18 h時開始進入發芽,在發芽早期種子會進行DNA修復以及一些復合體如核糖核蛋白體的水解,因此,PEG脅迫以后對蛋白質翻譯后修飾與轉運和翻譯,核糖體與生物合成功能的影響顯著增加;同時隨著脅迫時間的增加,次生代謝物質也在種子萌發過程中不斷增加。

圖4 差異基因表達譜

圖5 谷子萌發10 h和18 h差異基因的生物過程分類

差異基因被注釋到98條KEGG代謝通路中,廣泛涉及到物質代謝、能量代謝、信號轉導及次生代謝(圖4)。以-value≤0.05為篩選標準,發現這些基因在干旱處理10 h顯著的富集到16個代謝通路中,處理18 h顯著富集到12個代謝通路中(表4)。

干旱脅迫10 h,富集到植物激素信號轉導、氧化磷酸化、糖酵解、苯丙烷類代謝及植物病源菌的相互作用的差異表達基因達到了51%,其中植物激素信號轉導的差異表達的基因數目最多為10個。干旱脅迫18h,富集到苯丙烷類代謝、苯丙氨酸代謝、植物激素信號轉導及氧化磷酸化的差異基因達到51%,其中富集到苯丙烷類代謝通路的差異表達基因數目最多為10個。

植物內源激素廣泛參與到種子萌發過程,對種子的萌發進程起到一定調節作用。在干旱脅迫下,通過分析KEGG代謝通路發現,谷子萌發過程中關鍵的2個時間點10 h及18 h分別有10個和8個基因發生了差異表達,涉及生長素、乙烯、脫落酸、茉莉酸、水楊酸和油菜素內酯6種激素(圖8)。其中上調表達的基因包括絲氨酸/蘇氨酸蛋白激酶SnRK2、ABA反應元件結合因子ABF、乙烯受體ETR、轉錄因子TGA、生長素反應蛋白SAUR家族的X10A蛋白;下調表達的基因包括生長素載體AUX1 LAX家族、生長素響應的GH3家族、生長素反應蛋白SAUR家族的X10A和ARG7蛋白、ABA受體PYR/PYL家族、茉莉酸ZIM結構域蛋白質、木葡聚糖轉移酶TCH4。其中在萌發10和18 h都上調表達(圖5),其通常作為一個正調控因子參與到ABA信號轉導途徑,且ABA的含量在種子萌發過程中決定著種子是否萌發,大量的研究表明其抑制種子的萌發[34-35]。而其他的激素對于種子萌發的影響都是通過調節ABA并且與ABA形成互作網絡來實現的[32],因此推斷在干旱脅迫下調節種子的萌發。

通過分析KEGG代謝通路,還發現有較多數量的差異表達基因富集到苯丙烷類代謝和苯丙氨酸代謝。其中苯丙烷類代謝在干旱脅迫10和18 h的谷子萌發進程中最顯著,而苯丙氨酸代謝在干旱脅迫10 h時顯著性較差,在干旱脅迫18 h時顯著性僅次于苯丙烷類代謝。苯丙烷類代謝途徑是植物體內非常重要的次生代謝途徑,分析苯丙烷類代謝途徑(ko00940)發現苯丙氨酸解氨酶(PAL)、肉桂酸-4-羥化酶(CA4H)基因(表5)在谷子萌發10和18 h均發生了下調表達。其中PAL是苯丙烷類次級代謝的第一個關鍵酶[25-26],不僅受到內部發育的調節(活性隨著植物的的生長發育不斷發生變化)而且也受到外部條件的的誘導(植物在受到寒冷、傷害和病原菌感染都會激活苯丙烷代謝)。通過對擬南芥4個的研究中發現在沒有給植株澆水2周后,pla1pla2雙重突變體植物仍然是亮綠色,而野生型的植物出現萎蔫現象,說明和的缺失會提高植物的抗旱性[27]。通過分析發現在干旱脅迫下種子萌發過程中表達水平下降,從而使PAL酶活性降低,直接影響肉桂酸的含量而使其他次生代謝產物合成減少,進而緩解PEG脅迫對種子影響,因此,可以推測在種子萌發過程中起著關鍵作用。

圖6 谷子萌發10 h和18 h差異基因的分子功能分類

表4 差異表達基因富集的KEGG通路

2.6 熒光定量PCR分析

經過對顯著富集的GO功能分類條目和KEGG代謝通路的分析,篩選出5個差異表達基因,其中Si035342m.g、Si035908m.g、Si001109m.g和Si016504m.g下調表達;Si022448m.g上調表達。采用qRT-PCR分析其表達譜,發現5個差異基因的表達模式與其在表達譜中的表達趨勢一致,但是表達量的變化上有一定的不同(圖9)。

表5 苯丙烷代謝相關的差異表達基因

1:翻譯、核糖體結構與生物合成Translation, ribosomal structure and biogenesis;2:RNA加工與修飾RNA processing and modification;3:轉錄Transcription;4:復制、重組與修復Replication, recombination and repair;5:細胞周期調控與分裂、染色體重排Cell cycle control, cell division, chromosome partitioning;6:防御機制Defense mechanisms;7:信號轉導機制Signal transduction mechanisms;8:胞壁/膜生物發生Cell wall/membrane/envelope biogenesis;9:細胞骨架Cytoskeleton;10:胞內分泌與膜泡運輸Intracellular trafficking, secretion, and vesicular transport;11:蛋白質翻譯后修飾與轉運、分子伴侶Posttranslational modification, protein turnover, chaperones;12:能量產生和轉化Energy production and conversion;13:碳水化合物運輸與代謝Carbohydrate transport and metabolism;14:氨基酸運輸與代謝Amino acid transport and metabolism;15:核酸運輸與代謝Nucleotide transport and metabolism;16:輔酶運輸與代謝Coenzyme transport and metabolism;17:脂類運輸與代謝Lipid transport and metabolism;18:無機離子運輸與代謝Inorganic ion transport and metabolism;19:次生產物合成、運輸及代謝Secondary metabolites biosynthesis, transport and catabolism;20:一般功能基因General function prediction only;21:功能未知Function unknown

3 討論

植物受到干旱脅迫時,為了保護自己適應環境脅迫,不僅在生理生化代謝方面作出響應,同時在基因表達水平作出響應防御干旱脅迫,因此,從分子水平上來研究干旱脅迫響應的分子基礎,尋找響應干旱的關鍵基因,為以后利用植物基因工程手段改良作物的抗旱性奠定基礎。

種子在萌發過程尤其是在胚根突破種皮時會產生大量的活性氧,因此,種子會發生氧化應激反應,種子通過產生抗氧化物質和酶來調節這種反應[28]。QIE等[29]和TANG等[30]通過QTL定位發現種子萌發與干旱相關,并將相關信息和RNA-seq數據的聯合分析發現4個基因與種子萌發相關,其中包括了氧化應激反應和氧化還原反應。本研究通過GO分析發現氧化應激反應和氧化還原反應功能類別中富集了大量的差異基因,并且種子萌發到18 h富集到這兩個通路中的基因明顯增加。在前期的生理試驗中發現SOD的活性在持續增加,而CAT活性和POD活性呈現相反的趨勢,說明在晉谷45萌發受到到干旱脅迫時,調控SOD的基因對干旱的應答不夠靈敏高效,而調控CAT和POD的基因可以作出迅速的應答,從而導致SOD的活性增加,CAT和POD活性降低;玉米種子受到低溫脅迫后SOD、POD和CAT的活性整體是先增加后降低,而基因的表達量則是先減小后增加[31]。

A:10hCKvs10PEG;B:18hCKvs18hPEG

Tryptophan metabolism:色氨酸代謝;Auxin:生長素;Ubiquitin mediated proteolysis:泛素介導的蛋白水解;Cell enlargement:細胞增大;Plant growth:植物生長;Zeatin biosynthesis:玉米素的生物合成;Cytokinine:細胞分裂素;Cell division:細胞分裂;Shoot initiation:根的發生;Diterpenoid biosynthesis:二萜類化合物的生物合成;Gibberellin:赤霉素;Stem growth:莖的生長;Induced germination:誘導萌發;Carotenoid biosynthesis:類胡蘿卜素的生物合成;Abscisic acid:脫落酸;Stomatal closure:氣孔;Seed dormancy:種子休眠;Cysteine and methionine metabolism:半胱氨酸和蛋氨酸代謝;Ethylene:乙烯;Endoplasmic reticulum(ER):內質網;Fruit ripening:果實成熟;Senescence:衰老;Brassinosteroid biosynthesis:油菜素內酯的生物合成;Brassinosteroid:油菜素甾醇;Proteasomal degradation:蛋白酶體降解;Cell elongation:細胞伸長;Alpha-Linolenic acid metabolism:α-亞麻酸代謝;Jasmonic acid:茉莉酸;Monoterpenoid biosynthesis:類單萜生物合成;Indole alkaloid biosynthesis:吲哚生物堿的生物合成;Stress response:應激反應;Phenylalanine metabolism:苯丙氨酸代謝;Salicylic acid:水楊酸;Disease resistance:抗病性

圖8 植物激素信號轉導

Fig. 8 Plant hormone signal transduction under seed germinate for 10h and18h between control sample and PEG-stress sample

圖9 谷子萌發10 h和18 h差異基因的qRT-PCR驗證

激素是調節植物生長和發育的小分子,在種子萌發的過程中脫落酸、赤霉素、乙烯、生長素、水楊酸、茉莉酸等激素都會影響種子的萌發[32-33]。其中脫落酸可以抑制種子的萌發,而赤霉素可以解除這種抑制促進種子的萌發,其他激素單獨的調節種子的萌發的機制還不清楚,但是它們可以與脫落酸、赤霉素形成工作網絡一起來控制種子的萌發,尤其是在響應環境脅迫的過程中[34-35]。ABA與GA的比值決定著種子是否萌發。蔗糖不發酵相關蛋白激酶2(sucrose non- fermenting-1-related protein kinase2 SnRK2)參與環境脅迫及ABA等的信號轉導。在對擬南芥中SnRK2s的二重突變體SnRK2.2/3、三重突變體SnRK2.2/3/6和十重突變體SnRK2.1/2/3/4/5/6/7/8/9/10 研究中表明擬南芥中SnRK2基因家族有利于種子的萌發及提高擬南芥在干旱脅迫的適應能力[36-38]。本研究通過KEGG分析發現植物激素信號轉導途徑中ABA信號轉導中的下游基因的表達上調,推斷這將會導致ABA對種子萌發的抑制作用增強。同時現在已經證明當植物處于冷、鹽和滲透脅迫時,降低GA的水平和信號會有利于打破逆境對植物生長的限制[39]。影響GA合成的關鍵基因是基因家族,因此,在環境脅迫下,基因家族的表達就會影響GA的合成。GA生物合成過程中的、和基因家族受到GA含量的調節,當合成基因下調表達時,上調表達[39]。在本文表達譜分析中發現上調表達,推斷在干旱脅迫下下調表達,導致GA的含量下降,對種子的萌發的促進作用減弱。通過表達譜的初步分析可以推斷出在干旱條件下,ABA作用增強,而GA的作用減弱,它們的拮抗作用導致種子萌發延遲。

以上這些代謝過程或者基因,基本上對植物響應環境脅迫起到一定的作用,但是有些還未證實能夠調節種子萌發,那么這些基因能否在干旱條件下去幫助種子萌發也是有待進一步的驗證。同時通過對表達譜分析發現不只有已驗證的5個差異基因,而這5個基因也不只參與到一條代謝通路或者是信號轉導通路中去,所以谷子萌發期對干旱的響應是多種代謝途徑和信號轉導通路共同作用的結果。對谷子萌發期對干旱響應的研究,不是通過研究單個基因的功能就可以理解干旱脅迫下各個代謝途徑之間的相互關系。表達譜的分析為我們研究谷子萌發期的抗旱機制,篩選出起關鍵作用的基因提供理論基礎,需要利用基因工程手段對這些基因的功能進行進一步的驗證,從而為改良谷子的抗旱性提供依據。

4 結論

谷子萌發10和18 h分別篩選出456和545個差異基因,87和267個在干旱脅迫下上調表達,369和278個在干旱脅迫下下調表達。差異基因的表達廣泛涉及到糖、蛋白質、核酸等生物大分子代謝、次生代謝和能量代謝等過程。

[1] LATA C, PRASAD M. Setaria genome sequencing: an overview., 2013, 22(3): 257-260.

[2] 張雁明, 劉曉東, 馬建萍, 溫琪汾, 韓淵懷. 谷子抗旱研究進展. 山西農業科學, 2013, 41(3): 282-285.

ZHANG Y M, LIU X D, MA J P, WEN Q F, HAN Y H. Research progress on drought resistance in foxtail millet(L.)., 2013, 41(3): 282-285. (in Chinese)

[3] 崔紀菡, 范佳興, 李順國, 趙宇, 劉猛, 宋世佳, 任曉利, 劉斐, 南春梅, 夏雪巖. 谷子抗旱性鑒定研究進展. 東北農業大學學報, 2017, 48(1): 89-96.

CUI J H, FAN J X, LI S G, ZHAO Y, LIU M, SONG S J, REN X L, LIU P, NAN C M, XIA X Y. Research advance on evaluation of drought resistance of., 2017, 48(1): 89-96. (in Chinese)

[4] 裴帥帥, 尹美強, 溫銀元, 黃明鏡, 張彬, 郭平毅, 王玉國, 原向陽. 不同品種谷子種子萌發期對干旱脅迫的生理響應及其抗旱性評價. 核農學報, 2014, 28(10): 1897-1904.

PEI S S, YIN M Q, WENG Y Y, HUANG M J, ZHANG B, GUO P Y, WANG Y G, YUAN X Y. Physiological response of foxtail millet to drought stress during seed germination and drought resistance evaluation., 2014, 28(10): 1897-1904. (in Chinese)

[5] 代小冬, 楊育峰, 朱燦燦, 魯曉民, 王春義, 楊曉平, 楊國紅, 李君霞. 谷子萌芽期對干旱脅迫的響應及抗旱性評價. 華北農學報, 2015, 30(4): 139-144.

DAI X D, YANG Y F, ZHU C C, LU X M, WANG C Y, YANG X P, YANG G H, LI J X. Seed germination response to drought stress and drought resistance evaluation of foxtail millet., 2015, 30(4): 139-144. (in Chinese)

[6] LATA C, GUPTA S, PRASAD M. Foxtail millet: a model crop for genetic and genomic studies in bioenergy grasses., 2012, 33(3): 328-343.

[7] BENNETZEN J L, SCHMUTZ J, WANG H, PERCIFIELD R, HAWKINS J, PONTAROLI A C, ESTEP M, FENG L, VAUGHN J N, GRIMWOOD J, JENLINS J, BARRY K, LINDQUIST E, HELLSTEN U, DESHPANDE S, WANG X, WU X, MITROS T, TRIPLETT J, YANG X, YE C Y, MAURO-HERRERA M, WANG L, LI P, SHARMA M, SHARMA R, RONALD P C, PANAUD O, KELLOGG E A, BRUTNELL T P, DOUST A N, TUSKAN G A, ROKHSAR D, DEVOS K M R. eference genome sequence of the model plant setaria., 2012, 30(6): 555-561.

[8] ZHANG G, LIU X, QUAN Z, CHENGG S, XU X, PAN S, XIE M, ZENG P, YUE Z, WANG W, TAO Y, BIAN C, HAN C, XIA Q, PENG X, CAO R, YANG X, ZHAN D, HU J, ZHANG Y, LI H, LI H, LI N, WANG J, WANG C, WANG R, GUO T, CAI Y, LIU C, XIANG H, SHI Q, HUANG P, CHEN Q, LI Y, WANG J, ZHAO Z, WANG J. Genome sequence of foxtail millet () provides insights into grass evolution and biofuel potential., 2012, 30: 549-554.

[9] 趙晉鋒, 余愛麗, 田崗, 杜艷偉, 郭二虎, 刁現民. 谷子CBL基因鑒定及其在干旱、高鹽脅迫下的表達分析. 作物學報, 2013, 39(2): 360-367.

ZHAO J F, YU A L, TIAN G, DU Y W, GUO E H, DIAO X M. Identification of CBL genes from foxtail millet ([L.] Beauv.) and its expression under drought and salt stresses., 2013, 39(2): 360-367. (in Chinese)

[10] LATA C, SAHU P P, PRASAD M. Comparative transcriptome analysis of differentially expressed genes in foxtail millet (L.) during dehydration stress., 2010, 393(4): 720-727.

[11] 崔潤麗, 智慧, 王永芳, 李偉, 李海權, 黃占景, 刁現民. 谷子蛋白基因的克隆. 華北農學報, 2007, 22(4): 9-13.

CUI R L, ZHI H, WANG Y F, LI W, LI H Q, HUANG Z J, DIAO X M. Cloning of-like protein gene from foxtail millet., 2007, 22(4): 9-13. (in Chinese)

[12] 楊希文, 胡銀崗. 谷子轉錄因子基因的克隆及其在干旱脅迫下的表達模式分析. 干旱地區農業研究, 2011, 29(5): 69-74.

YANG X W, HU Y G. Cloning of agene from foxtail millet (L.) and its expression during drought stress., 2011, 29(5): 69-74. (in Chinese)

[13] 張彬, 禾璐, 候蕊, 路陽, 馬芳芳, 王興春, 楊致榮, 韓淵懷, 李紅英. 谷子C2H2型鋅指蛋白基因的克隆及表達分析. 中國農業大學學報, 2015, 20(5): 9-15.

ZHANG B, HE L, HOU R, LU Y, MA F F, WANG X C, YANG Z R, HAN Y H, LI H Y. Isolation and functional analysis of a C2H2-type zinc finger protein gene, 2015, 20(5): 9-15. (in Chinese)

[14] 竇祎凝, 秦玉海, 閔東紅, 張小紅, 王二輝, 刁現民, 賈冠清, 徐兆師, 李連城, 馬有志, 陳明.谷子轉錄因子8通過ABA信號途徑正向調控干旱條件下的種子萌發. 中國農業科學, 2017, 50(16): 3071-3081.

DOU Y N, QIN Y H, MIN D H, ZHANG X H, WANG E H, DIAO X M, JIA G Q, XU Z S, LI L C, MA Y Z, CHEN M. Transcription factorpositively regulates seed germination under drought stress through ABA signaling pathway in Foxtail Millet (L.)., 2017, 50(16): 3071-3081. (in Chinese)

[15] 王一帆, 李臻, 潘教文, 李穎秀, 王慶國, 管延安, 劉煒. 谷子基因克隆及功能分析.遺傳, 2017, 39(5): 413-422.

WANG Y F, LI Z, PAN J W, LI Y X, WANG Q G, GUAN Y A, LIU W. Cloning and functional analysis of thegene inL.., 2017, 39(5): 413-422. (in Chinese)

[16] LANGMEAD B, TRAPNELL C. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome., 2009, 10(3): R25.

[17] MORTAZAVI A, WILLIAMS B A, MCCUE K, SCHAEFFER L, WOLD B. Mapping and quantifying mammalian transcriptomes by RNA-Seq., 2008, 5(7): 621-628.

[18] Anders S, Huber W. Differential expression analysis for sequence count data., 2010, 11(10): R106.

[19] 鄧泱泱, 荔建琦, 吳松鋒, 朱云平, 陳耀文, 賀福初. nr數據庫分析及其本地化. 計算機工程, 2006, 32(5): 71-74.

DENG Y Y, LI J Q, WU S F, ZHU Y P, CHEN Y W, HE F Z. Integrated nr database in protein annotation system and its localization., 2006, 32(5): 71-74. (in Chinese)

[20] APWEILER R, BAIROCH A, WU C H, BARKER W C, BOECKMANN B, FERRO S, GASTEIGER E, HUANG H, LOPEZ R, MAGRANE M, MARTIN M J, NATALE D A, O'DONOVAN C, REDASCHI N, YEH L S. UniProt: the universal protein knowledgebase., 2004, 32: D115-D119.

[21] ASHBURNER M, BALL C A, BLAKE J A, BOTSTEIN D, BUTLER H, CHERRY J M, DAVIS A P, DOLINSKI K, DWIGHT S S, EPPIG J T, HARRIS M A, HILL D P, ISSEL-TARVER L, KASARSKIS A, LEWIS S, MATESE J C, RICHARDSON J E, RINGWALD M, RUBIN G M, SHERLOCK G. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium.,2000,25(1): 25-29.

[22] TATUS R L, GALPERIN M Y, NATALE D A, KOONIN E V. The COG database: a tool for genome-scale analysis of protein functions and evolution.,2000,28(1): 33-36.

[23] KANEHISA M, GOTO S, KAWASHIMA S, OKUNO Y, HATTORI M. The KEGG resource for deciphering the genome., 2004(32): D277-D280.

[24] ALEXA A, RAHNENFUHRER J. topGO: enrichment analysis for gene ontology. R package version 2.8, 2010.

[25] 李莉, 趙越, 馬君蘭. 苯丙氨酸代謝途徑關鍵酶: PAL、C4H4、CL研究新進展. 生物信息學, 2007, 5(4): 187-189.

LI L, ZHAO Y, MA J L. Recent progress on key enzymes: PAL, C4H, 4CL of phenylalanine metabolism pathway., 2007, 5(4): 187-189. (in Chinese)

[26] 徐曉梅, 楊署光. 苯丙氨酸解氨酶研究進展. 安徽農業科學, 2009, 37(31): 15115-15119.

XU X M, YANG S G. Advances in the studies of phenylalanine ammonialyase., 2009, 37(31): 15115-15119. (in Chinese)

[27] HUANG J, GU M, LAI Z, FAN B, SHI K, ZHOU Y H, YU J Q, CHEN Z. Functional analysis of the Arabidopsis PAL gene family in plant growth,development, and response to environmental stress., 2010, 153(4): 1526.

[28] Bailly C. Active oxygen species and antioxidants in seed biology., 2004, 14: 93-107.

[29] QIE L F, JIA G Q, ZHANG W Y, JAMES S, SHANG Z L, LI W, LIU B H, LI M Z, CHAI Y, ZHI H, DIAO X M. Mapping of quantitative trait locus (QTLs) that contribute to germination and early seedling drought tolerance in the interspecific cross setaria italic×setaria viridis., 2014, 9(7): e101868 .

[30] TANG S, LI L, WANG Y Q, CHEN Q N, ZANG W Y, JIA G Q, ZHI H, ZHAO B H, DIAO X M. Genotype specific physiological and transcriptomic responses to drought stress in(an emerging model for Panicoideae grasses)., 2017, 7(1): 10009.

[31] 楊靜, 毛笈華, 于永濤, 李春艷, 王永飛, 胡建廣. 低溫對甜玉米種子氧化酶活性的影響及相關基因表達分析. 核農學報, 2016, 30(9): 1840-1847.

YANG J, MAO J H, YU Y T, LI C Y, WANG Y F, HU J G. Effects of chilling on antioxidant enzyme activity and related gene expression levels during seed germination., 2016,30(9): 1840-1847. (in Chinese)

[32] MIRANSARI M, SMITH D L. Plant hormones and seed germination., 2014, 99(3): 110-121.

[33] SHU K, LIU X D, XIE Q, HE Z H. Two faces of one seed: Hormonal regulation of dormancy and germination., 2016, 9(1): 34-45.

[34] LOIC R, MANUAL D, KARINE G, JULIE C, JULIA B, CLAUDETTE J, DOMINIQUE J.Seed germination and vigor., 2012, 63: 507-533.

[35] LEONIE B, MAARTEN K. Seed Dormancy and Germination., 2008, 12(30): e0119.

[36] FUJII H, ZHU J K, JAGENDORF A T.mutant deficient in 3 abscisic acid-activated protein kinases reveals critical roles in growth, reproduction, and stress., 2009, 106(20): 8380-8385.

[37] FUJII H, VERSLUES P E, ZHU J K. Identification of two protein kinases Required for abscisic acid regulation of seed germination, root growth and gene expression in., 2007, 19(2): 485-494.

[38] FUJII H, ZHU J K.decuple mutant reveals the importance of SnRK2 kinases in osmotic stress responses., 2011, 108(4): 1717-1722.

[39] COLEBROOK E H, THOMAS S G, PHILLIPS A L, HEDDEN P. The role of gibberellin signalling in plant responses to abiotic stress., 2014, 217(1): 67-75.

附表1 PEG脅迫處理下谷子種子萌發10 h和18 h時均與對照共表達的20個差異基因的表達情況

Supplemental table 1 Expression level of 20 common DEGs among 10hCKvs10PEG, 18hCKvs18hPEG, and 10hCKvs18CK

附表1 PEG脅迫處理下谷子種子萌發10 h和18 h時均與對照共表達的20個差異基因的表達情況

Supplemental table 1 Expression level of 20 common DEGs among 10hCKvs10PEG, 18hCKvs18hPEG, and 10hCKvs18CK

附表2 PEG脅迫處理下谷子種子萌發10 h時與對照差異基因的表達情況(10hCKvs10hPEG)

Supplemental table 2 DEGs under seed germinate for 10 h between control sample and PEG-stress sample

附表2 PEG脅迫處理下谷子種子萌發10 h時與對照差異基因的表達情況(10hCKvs10hPEG)

Supplemental table 2 DEGs under seed germinate for 10 h between control sample and PEG-stress sample

附表3 PEG脅迫處理下谷子種子萌發18 h時與對照差異基因的表達情況(18hCKvs18hPEG)

Supplemental table 3 DEGs under seed germinate for 18 h between control sample and PEG-stress sample

附表3 PEG脅迫處理下谷子種子萌發18 h時與對照差異基因的表達情況(18hCKvs18hPEG)

Supplemental table 3 DEGs under seed germinate for 18 h between control sample and PEG-stress sample

請根據該文DOI:10.3864/j.issn.0578-1752.2018.08.002 或掃描右側二維碼,從相應網頁“附件”處下載查閱以上信息(手機或不兼容)。

Gene Expression Profiling of Foxtail millet (L.) under Drought Stress during Germination

XU BingXia, YIN MeiQiang, WEN YinYuan, PEI ShuaiShuai, KE ZhenJin, ZHANG Bin, YUAN XiangYang

(College of Agriculture, Shanxi Agricultural University, Taigu 030801, Shanxi)

【Objective】Foxtail millet (L.) is a drought tolerant crop. The objective of this research is to get a lot of differently expressed genes during germination in response to drought stress by high-throughput sequencing, then to obtain the key gene and the related molecular mechanism at seed germination stage in foxtail millet under drought stress.【Method】Seed of JinGu45 was treated with 18%PEG-6000(PEG-stress)and distilled water germination(control sample) at 1 h, 10 h and 18 h as a test material, and the activities of SOD, POD and CAT were measured,respectively. SOD activity was assayed by nitro blue tetrazolium (NBT) method. POD activity was determined by guaiacol method, and the activity of CAT was measured by colorimetric method. Sample of control and PEG-stress that germinated for 10h and 18h were used to construct cDNA library by gene expression profiling technology. We compared reads to the reference genome by using Bowtie and analysed the result by using RSEM. Differential expression analysis used DESeq. The functional annotation of differently expresse genes were obtained by using NR, Swiss-Prot, KEGG, COG and GO online databases. The key genes that regulate germination in foxtail millet was obtained through analyse DEGs. The reliability of sequencing results was comfirmed by qRT-PCR. 【Result】The SOD activity of PEG-stress sample was higher than that in the control sample, but the activity of POD and CAT were lower than that in the control group. With the time of germination changes, the activity of SOD was increased, but the activity of CAT and POD was gradually decreased. The sequences of gene expression profile was highly consistent with the selected reference genome sequence, and the gene expression was highly heterogeneous. Expression analysis showed a total of 35 470 genes, and with the selection criteria of RPKM ≥0.01, there were 24 030 and 24 486 genes in the control samples and 24 019 and 23 877 genes in the samples under PEG drought stress, respectively. 456 and 545 DEGs were screened out during millet germination at 10h and 18h under drought stress, in which 87 and 267 DEGs were up-regulated and 369 and 278 DEGs were down-regulated. GO enrichment analysis showed that these DEGs were mainly relating to metabolism process, cell stimulation and response process. The KEGG enrichment analysis showed that these DEGs were associated with phenylpropanoid metabolism and plant hormone signal transduction. The results obtained from five genes tested by RT-PCR agreed with the trend of regulation identified by gene expression profile. 【Conclusion】DEGs were widely involved in the metabolism of biomacromolecule such as sugar, protein, nucleic acid, secondary metabolism and energy metabolism.andgenes may regulate seed germination in foxtail millet under drought stress.

foxtail millet (L.); drought stress; seed germination; gene expression profile

(責任編輯 李莉)

10.3864/j.issn.0578-1752.2018.08.002

2017-11-03;

2017-12-26

山西省重點研發計劃(201603D221003-2)、作物生態與旱作栽培生理山西省重點實驗室項目(201705D111007)、山西省科技重點研發項目(2015-TN-09)

許冰霞,Tel:18734450502;E-mail:1986984710@qq.com。

尹美強,Tel:15935411700;E-mail:yinmq999@163.com

猜你喜歡
差異分析
相似與差異
音樂探索(2022年2期)2022-05-30 21:01:37
隱蔽失效適航要求符合性驗證分析
找句子差異
DL/T 868—2014與NB/T 47014—2011主要差異比較與分析
電力系統不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
生物為什么會有差異?
電力系統及其自動化發展趨勢分析
M1型、M2型巨噬細胞及腫瘤相關巨噬細胞中miR-146a表達的差異
中西醫結合治療抑郁癥100例分析
收入性別歧視的職位差異
主站蜘蛛池模板: 亚洲色偷偷偷鲁综合| 亚洲日本中文字幕天堂网| 亚洲国产黄色| 成年免费在线观看| 美女扒开下面流白浆在线试听| 欧美日韩午夜| 亚洲av色吊丝无码| 日韩av无码DVD| 国产亚洲男人的天堂在线观看 | 久99久热只有精品国产15| 2022国产91精品久久久久久| 欧美亚洲国产日韩电影在线| 久久久国产精品无码专区| 国产激情无码一区二区免费| 麻豆国产精品视频| 久久五月天国产自| 无码精油按摩潮喷在线播放| www.亚洲一区| 高清大学生毛片一级| 亚洲aⅴ天堂| 国产午夜无码专区喷水| 国产手机在线观看| 91网在线| 欧美日韩成人在线观看| 亚洲天堂免费在线视频| 午夜精品区| 久久综合九色综合97网| 国产精选小视频在线观看| 精品国产免费观看一区| 毛片三级在线观看| 高清无码手机在线观看| 亚洲av无码成人专区| 制服丝袜一区| 日韩精品少妇无码受不了| 亚洲国产成人在线| 久久久久久久久久国产精品| 91精品伊人久久大香线蕉| 国产全黄a一级毛片| 国产乱子伦精品视频| 国内精品久久久久久久久久影视 | 国产人人射| 永久免费av网站可以直接看的| 亚洲精品无码久久毛片波多野吉| 免费一级毛片在线观看| 亚洲第一黄片大全| 日韩 欧美 小说 综合网 另类| www.狠狠| 国产精品免费露脸视频| 免费人成网站在线高清| 精品综合久久久久久97| 亚洲中文字幕97久久精品少妇| 国产午夜一级毛片| 国产永久免费视频m3u8| 91福利在线看| 一本久道热中字伊人| 国产精品自在在线午夜区app| 91欧洲国产日韩在线人成| 中文字幕啪啪| 黄色国产在线| 中文字幕波多野不卡一区| 91免费国产在线观看尤物| 国产亚洲精品无码专| 91小视频在线观看免费版高清| 亚洲欧美日本国产综合在线 | 高清欧美性猛交XXXX黑人猛交 | 国产久草视频| 国产自在线播放| 啊嗯不日本网站| 国产一区二区三区免费| 欧美精品影院| 91在线国内在线播放老师| 国产激情第一页| 国产白丝av| 欧美翘臀一区二区三区| 亚洲天堂精品在线观看| 欧美色视频在线| a毛片基地免费大全| 成人福利在线看| 亚洲一区二区精品无码久久久| 女人18毛片水真多国产| 日韩人妻无码制服丝袜视频| 日韩午夜片|