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

基于QTL和轉錄組測序鑒定甘藍型油菜耐旱候選基因

2024-03-28 05:46:18李陽陽許軍紅陳倬永徐昕媛徐金盼唐鐘林張婭茹嚴卓立周清元李加納劉列釗唐章林
作物學報 2024年4期

李陽陽 吳 丹 許軍紅 陳倬永 徐昕媛 徐金盼 唐鐘林 張婭茹 朱 麗 嚴卓立 周清元 李加納 劉列釗 唐章林,*

基于QTL和轉錄組測序鑒定甘藍型油菜耐旱候選基因

李陽陽1,2,3吳 丹2,3許軍紅2,3陳倬永1,2,3徐昕媛1,2,3徐金盼1,2,3唐鐘林1,2,3張婭茹1,2,3朱 麗1,2,3嚴卓立1,2,3周清元1,2,3李加納1,2,3劉列釗1,2,3唐章林1,2,3,*

1西部(重慶)科學城種質創制大科學中心, 重慶 401329;2西南大學農學與生物科技學院, 重慶 400715;3西南大學農業科學研究院, 重慶 400715

干旱脅迫嚴重限制了甘藍型油菜種植面積的擴大和產量的提升。耐旱性是由多基因控制的復雜數量性狀, 將QTL定位與轉錄組測序相結合, 是鑒定甘藍型油菜耐旱候選基因的有效手段。本研究對甘藍型油菜干旱敏感品系三六矮和耐旱品系科里納-2構建的F2:6和F2:8重組自交系群體幼苗進行正常灌溉和干旱脅迫處理, 測定地上部鮮重、地上部干重、葉片相對含水量、丙二醛和可溶性糖含量, 利用SSR和SNP多態性分子標記構建遺傳連鎖圖譜, 鑒定耐旱相關QTL和候選區間, 結合耐旱材料No11和干旱敏感材料No28的轉錄組測序, 篩選耐旱相關候選基因。研究結果表明: 干旱脅迫使甘藍型油菜幼苗地上部鮮重、地上部干重和葉片相對含水量下降, 使葉片丙二醛和可溶性糖含量上升; 耐旱相關QTL和候選區間分布于A01、A02、A06、A08、A09、A10、C02、C03、C04、C06和C09染色體; 對耐旱材料和干旱敏感材料正常灌溉、干旱24 h、36 h和48 h進行轉錄組分析, 主要差異表達基因顯著富集到光合作用、脂肪酸代謝、氨基酸代謝、植物激素信號轉導、核糖體、晝夜節律及角質、木栓素和蠟質的生物合成等相關途徑; 將QTL與轉錄組測序相結合, 鑒定到28個耐旱相關候選基因, 主要編碼FLC、bHLH105、TGA4、TEM1、ERF003、ACO3、CHLI1、LHCB6和PORC等, 具有轉錄因子活性、乙烯產生和信號傳導、葉綠素生物合成與結合、葉綠素氧化還原酶以及編碼核糖體相關蛋白等功能。這些結果可為揭示甘藍型油菜耐旱機理及分子標記輔助選育耐旱新品種奠定基礎。

甘藍型油菜; 耐旱; QTL; 轉錄組; 候選基因

油菜是世界主要的油料作物之一, 也是我國食用植物油的主要來源[1-2]。長江流域是我國油菜主產區[3], 但年內降雨分布不均, 秋、冬季降雨較少, 季節性干旱嚴重[4], 致使油菜苗期常遭受干旱脅迫。干旱脅迫導致氣孔關閉, 蒸騰速率降低, 光合速率下降, 嚴重影響有機物的積累, 導致苗小葉少, 難以形成壯苗越冬, 影響后期生長發育導致產量下降[5]。甘藍型油菜(L.)因產量和含油量高且穩定被廣泛種植, 而干旱脅迫的頻發嚴重限制了其產量的提高和種植面積的擴大。

甘藍型油菜耐旱性是非常復雜的數量性狀, 由多基因控制[6]。前期, 較多報道用AFLP或SSR標記構建遺傳連鎖圖譜, 定位到多個甘藍型油菜耐旱相關的QTL。例如, 李真[7]基于雙單倍體(doubled haploid, DH)群體, 利用183個SSR標記和157個AFLP標記構建遺傳連鎖圖譜, 在正常灌溉和干旱脅迫兩種水分條件下定位到50個QTL, 利用耐旱系數定位到19個QTL。王丹丹[8]利用重組自交系(recombinant inbred line, RIL)群體, 構建了一張包含120個SSR標記的遺傳連鎖圖譜, 定位到8個耐旱相關QTL。但是, 大部分研究未鑒定到耐旱相關候選基因。近年來, 隨著甘藍型油菜基因組測序的完成, SNP、Indel等標記相繼被開發, 通過QTL定位、全基因組關聯分析(genome wide association study, GWAS)方法, 鑒定到一些耐旱候選基因。Gad等[9]基于DH群體, 利用SNP和Indel標記鑒定到39個甘藍型油菜種子萌發期耐旱相關QTL和256個候選基因。Shahzad等[10]通過甘藍型油菜GWAS, 鑒定到139個失水速率相關SNP位點和4個候選基因。Khanzada等[11]對228份甘藍型油菜種質資源幼苗進行PEG模擬干旱脅迫處理, 通過GWAS鑒定到314個顯著關聯SNP標記和85個候選基因。Tan等[12]對甘藍型油菜鹽脅迫和PEG-6000模擬干旱脅迫下種子發芽率和發芽指數進行GWAS, 鑒定到43個顯著關聯SNP標記和37個候選基因。通過轉錄組測序, 可獲得物種在某種狀態下的所有轉錄本信息[13], 已被廣泛應用于甘藍型油菜耐旱相關研究。Liu等[14]對正常灌溉和干旱脅迫下甘藍型油菜根和葉進行轉錄組測序, 在根和葉中分別發現6018個和5377個差異表達基因。Wang等[15]通過甘藍型油菜耐旱和干旱敏感材料轉錄組測序, 篩選到169個差異表達基因。

已有QTL或GWAS研究所用群體及種植環境或脅迫處理方式不同, 檢測到的QTL和候選基因不盡相同; 轉錄組測序結果差異表達基因數量龐大, 難以篩選出進一步開展功能驗證的目標基因。因此, 將QTL或GWAS與轉錄組測序相結合, 可提高特定性狀候選基因的鑒定效率, 已成為預測候選基因的有效手段。Zhou等[16]將GWAS和轉錄組分析相結合, 鑒定出甘藍型油菜耐鋁性候選基因。Jian等[17]結合QTL與轉錄組分析, 鑒定了甘藍型油菜花期相關候選基因。Guo等[18]通過GWAS和轉錄組揭示了干旱脅迫下與玉米幼苗根長相關的候選基因。Sevanthi等[19]通過干旱和長期缺氮雙脅迫下轉錄組測序和主效QTL整合挖掘到水稻關鍵脅迫響應基因。但截至目前, 尚未見有關于將QTL定位與轉錄組測序相結合用于鑒定甘藍型油菜耐旱候選基因的報道。

本研究基于甘藍型油菜干旱敏感品系三六矮與耐旱品系科里納-2構建的F2:6和F2:8重組自交系群體, 利用多態性SSR和SNP標記分別構建遺傳連鎖圖譜, 對耐旱相關性狀進行QTL定位, 獲得候選區間。對實驗室篩選的甘藍型油菜耐旱材料No11和干旱敏感材料No28進行轉錄組測序和差異表達基因分析, 結合QTL定位候選區間篩選出耐旱相關候選基因, 以期為進行功能驗證揭示甘藍型油菜耐旱機理和分子標記輔助選育耐旱新品種奠定基礎。

1 材料與方法

1.1 試驗材料與處理

本研究以盆缽(直徑25 cm, 高30 cm)方式將甘藍型油菜干旱敏感品系三六矮和耐旱品系科里納-2構建的105個F2:6重組自交系家系種植于西南大學7號溫室, 126個F2:8重組自交系家系種植于重慶市油菜工程技術研究中心歇馬油菜基地旱棚。以均勻一致的大田表層土作為種植土, 定期間苗, 保持幼苗長勢一致。當幼苗長至四葉一心時, 進行干旱脅迫(drought stress, DS)處理, 以正常灌溉為對照(control check, CK)。干旱脅迫組土壤含水量保持在9%~12%, 正常灌溉組土壤含水量保持在22%~27%, 利用浙江托普儀器有限公司的多參數土壤水分記錄儀, 每隔3~5 d測定土壤含水量, 儀器探針入土深度為7 cm, 每盆測量3個不同位置取平均值, 根據裝盆土壤干重計算應澆水量, 用燒杯均勻澆于盆缽表面。

甘藍型油菜耐旱材料No11和干旱敏感材料No28以盆栽(長5 cm, 寬5 cm, 高5 cm)方式種植于RDN-300C-4 (寧波東南)人工氣候箱, 腐植土和蛭石按照3∶1比例混合均勻后作為種植土, 每盆種植3株幼苗。正常灌溉期間, 光合有效輻射為270mmol m–2s–1, 光照和黑暗時長分別為16 h和8 h, 溫度分別為24℃和22℃, 濕度為70%~80%。待幼苗長至三葉一心時, 選取長勢一致的盆缽進行正常灌溉和干旱脅迫處理。干旱脅迫處理期間, 濕度參數為30%, 其他培養條件與正常灌溉相同。

1.2 取材與性狀測定

干旱處理28 d后, 調查正常灌溉和干旱脅迫下甘藍型油菜重組自交系各家系地上部鮮重(shoot fresh weight, SFW)、地上部干重(shoot dry weight, SDW)、葉片相對含水量(relative water content, RWC)、可溶性糖(soluble sugar, SUG)含量和丙二醛 (malondialdehyde, MDA)含量。F2:6重組自交系調查2個生物學重復, F2:8重組自交系調查3個生物學重復。對耐旱材料No11和干旱敏感材料No28正常灌溉及干旱后24 h、36 h和48 h進行RWC測定, 并取材進行SUG和MDA含量測定以及轉錄組測序。

1.2.1 地上部鮮重和干重測定 將植株從子葉節處剪開, 稱其地上部鮮重, 然后于110℃殺青30 min后75℃烘至恒重, 稱其地上部干重。

1.2.3 可溶性糖和丙二醛含量測定 選取新鮮倒數第2片葉, 采用荊蓉蓉[20]所用的蒽酮-硫酸法和硫代巴比妥酸法分別測定可溶性糖和丙二醛含量, 略有改動。

1.3 表型數據處理

利用R-4.3.1 (https://cloud.r-project.org/)進行方差分析, Microsoft Excel 2016進行性狀描述統計分析, 用下式計算各耐旱相關性狀的耐旱系數(drought resistance index, DRI):

1.4 SSR標記帶型及SNP標記基因型分析

本研究所采用的SSR引物均由西南大學重慶市油菜工程技術研究中心實驗室保存提供, 先對親本進行聚丙烯酰胺凝膠電泳檢測, 將帶型清晰且有多態性的引物對F2:6重組自交系群體進行電泳檢測, 將SSR引物于BnTIR (https://yanglab.hzau.edu.cn/ BnTIR)網站以Darmor. v4.1為參考基因組進行序列比對, 獲取SSR標記物理位置; 以Darmor. v4.1為參考基因組序列, 采用中玉金標記生物有限公司(北京) 6K SNP芯片對親本及F2:8重組自交系群體進行檢測和帶型分析。用“a”代表與母本三六矮相同的帶型(基因型), “b”代表與父本科里納-2相同的帶型(基因型), “h”代表雜合帶型(基因型), “-“代表缺失[8]。

1.5 遺傳連鎖圖譜構建及QTL定位

采用JoinMap 4[21]分別進行SSR標記和SNP標記遺傳連鎖圖譜構建, 采用MapQTL 6[22]對耐旱相關性狀在正常灌溉和干旱脅迫下的表型值及耐旱系數進行QTL定位(取LOD值≥2.5), 利用MapChart 2.0[23]軟件繪制QTL在各連鎖群(染色體)的分布圖。QTL命名以斜體加性狀的英文縮寫加處理組或耐旱系數英文縮寫再加染色體命名, 例如“”[24]。

根據李陽陽等[25]連鎖不平衡分析的甘藍型油菜各染色體的衰減距離, 將QTL包含標記物理位置減去所在染色體衰減距離作為候選區間上游位點, 加上衰減距離作為候選區間下游位點, 以此獲得候選區間。

1.6 轉錄組分析

由廣州基迪奧生物科技有限公司完成葉片RNA提取和質量檢測、文庫制備、文庫質檢、上機測序、序列比對分析及基因表達量分析, 參考基因組為法國油菜(__v5)。

對基因表達量進行兩兩差異表達基因分析, 保留值≤0.01和| log2(FC) |≥1的基因。使用OmicShare云平臺韋恩圖在線軟件(https://www. omicshare.com/tools/Home/Soft/venn)繪制韋恩圖, 采用百邁客云平臺(https://www.biocloud.net/) GO、KEGG分類富集軟件進行KEGG富集分析, 使用OmicShare云平臺熱圖在線軟件(https://www. omicshare.com/tools/Home/Soft/heatmap)繪制熱圖。采用百邁客云平臺(https://www.biocloud.net/) WGCNA在線分析軟件進行WGCNA分析。

2 結果與分析

2.1 表型數據分析

對甘藍型油菜F2:6和F2:8重組自交系群體的5個耐旱相關性狀進行方差分析和描述統計分析發現,各性狀在水分間、家系間以及水分×家系互作間均達到極顯著差異(表1); 干旱處理后, 葉片相對含水量、地上部鮮重和地上部干重較對照極顯著降低, 可溶性糖含量和丙二醛含量極顯著升高(表2)。頻數分布直方圖顯示各性狀均呈連續近似正態分布(圖1)。

表1 甘藍型油菜F2:6和F2:8重組自交系群體耐旱相關性狀的方差分析

RWC: 葉片相對含水量; SFW: 地上部鮮重; SDW: 地上部干重; SUG: 可溶性糖含量; MDA: 丙二醛含量; DF: 自由度;-value:值;**:< 0.01。

RWC: the relative water content; SFW: shoot fresh weight; SDW: shoot dry weight; SUG: soluble sugar content; MDA: malondialdehyde content; DF: degree of freedom;value:-value of-test;**:< 0.01.

表2 甘藍型油菜親本及F2:6和F2:8重組自交系群體耐旱相關性狀描述統計分析

(續表2)

CK: 正常灌溉; DS: 干旱脅迫。其他縮寫同表1。

CK: well water; DS: drought stress. The other abbreviations are the same as those given in Table 1.

圖1 甘藍型油菜F2:6和F2:8重組自交系群體耐旱相關性狀的頻數分布圖

縮寫同表1和表2。Abbreviations are the same as those given in Tables 1 and 2.

測定耐旱和干旱敏感材料(No11和No28)葉片相對含水量、丙二醛和可溶性糖含量發現, 干旱24 h、36 h和48 h, 干旱敏感材料葉片相對含水量極顯著低于耐旱材料; 干旱36 h和48 h, 干旱敏感材料丙二醛含量極顯著高于耐旱材料, 而干旱24 h兩材料間差異不顯著; 干旱24 h和48 h, 干旱敏感材料可溶性糖含量極顯著高于耐旱材料, 而干旱36 h差異未達到顯著水平(圖2)。說明干旱敏感材料較耐旱材料有更嚴重的失水和膜脂過氧化, 也積累了更多的滲透調節物質。

2.2 遺傳連鎖圖譜構建及QTL定位

本研究對親本及F2:6重組自交系群體進行SSR多態性標記檢測, 共獲得301個多態性SSR標記, 對其進行遺傳連鎖圖譜構建, 獲得了包含14個連鎖群和181個多態性SSR標記的甘藍型油菜遺傳連鎖圖譜, 該圖譜總長717.88 cM, 平均密度為5.69 cM (圖3)。利用6K SNP芯片對親本及F2:8重組自交系群體進行檢測, 得到1074個多態性SNP位點, 構建了總長為2577.03 cM的甘藍型油菜遺傳連鎖圖譜, 該圖譜覆蓋19條染色體, 包含1041個多態性SNP位點, 平均密度為2.48 cM (圖4)。

利用SSR標記遺傳連鎖圖譜對F2:6重組自交系群體的耐旱相關性狀進行QTL檢測, 共獲得6個QTL。其中正常灌溉下地上部干重檢測到1個QTL, 包含2個A02染色體的SSR標記; 正常灌溉下丙二醛含量和干旱處理下葉片相對含水量檢測到置信區間重疊的2個QTL, 包含1個A09染色體的SSR標記; 干旱處理下葉片相對含水量、葉片相對含水量耐旱系數和干旱處理下丙二醛含量檢測到置信區間重疊的3個QTL, 包含4個SSR標記, 其中3個位于C02染色體, 另1個位于C02或A02染色體(圖3和表3)。

圖2 甘藍型油菜耐旱和干旱敏感材料耐旱相關性狀

**:< 0.01。縮寫同表1。**:< 0.01. Abbreviations are the same as those given in Table 1.

圖3 SSR標記遺傳連鎖圖譜及QTL分布

DRI: 耐旱系數。其他縮寫同表1和表2。

DRI: drought resistance index. The other abbreviations are the same as those given in Tables 1 and 2.

表3 甘藍型油菜耐旱相關性狀QTL及候選區間

(續表3)

利用SNP標記遺傳連鎖圖譜對F2:8和F2:6重組自交系群體進行QTL檢測, 共檢測到18個QTL, 分布于A01、A02、A06、A08、A09、A10、C02、C03、C04、C06和C09染色體, 包含44個SNP位點。其中, 基于F2:6重組自交系群體正常灌溉下地上部干重和地上部鮮重在A09染色體檢測到相同的QTL位點, 基于F2:8重組自交系群體葉片相對含水量耐旱系數和干旱脅迫下葉片相對含水量在A02染色體檢測到的QTL位點置信區間相重合(圖4和表3)。

根據各染色體衰減距離, 共獲得52個候選區間,其中QTL位點、、、、所鑒定到的多個候選區間存在重疊(表3)。基于SSR標記遺傳連鎖圖譜定位到的QTL位點、和與基于SNP標記遺傳連鎖圖譜定位到的QTL位點存在候選區間重疊。

2.3 轉錄組分析及候選基因鑒定

對耐旱和干旱敏感材料正常灌溉及干旱24 h、36 h、48 h的幼苗進行轉錄組測序發現, 2個材料的干旱處理與對照間均存在至少7000個差異表達基因(圖5-A); 其中干旱48 h時, 兩材料間差異表達基因最多, 耐旱材料較干旱敏感材料上調和下調表達的基因分別有4631個和5825個(圖5-B)。在所有干旱處理條件下, 兩材料間均存在的表達差異基因有1723個(圖5-C), 2個材料的3個干旱時期均與對照存在表達差異的基因有1459個(圖5-D), 對此3182個核心差異表達基因進行KEGG富集分析, 顯著富集在賴氨酸降解及角質、木栓素和蠟質的生物合成途徑(圖5-E)。在核心差異表達基因的轉錄因子編碼基因和顯著富集的KEGG通路基因中, 有8個位于QTL定位鑒定的候選區間, 可認定為耐旱相關候選基因。其中, 編碼ERF003的和在干旱處理下較正常灌溉下調表達; 干旱脅迫下, 耐旱材料的()、()和()基因表達水平低于干旱敏感材料,()和()基因表達水平高于干旱敏感材料;()在干旱處理24 h和36 h, 耐旱材料表達水平低于干旱敏感材料, 在干旱處理48 h高于干旱敏感材料(圖5-F)。

基于耐旱和干旱敏感材料的表型和差異基因表達水平進行WGCNA分析(圖6)發現, 各差異表達基因被分為14個模塊, 與丙二醛含量最為正相關的模塊為MEfirebrick2, 最為負相關的模塊為MEdeeppink1; 與可溶性糖含量最為正相關的模塊為MEplum, 最為負相關的模塊為MEdeeppink1; 與葉片相對含水量最為正相關的模塊為MEindianred4, 最為負相關的模塊為MEfirebrick2。

圖5 正常灌溉與干旱脅迫處理下甘藍型油菜耐旱和干旱敏感材料差異表達基因分析

A和D: 兩材料各個干旱處理時期與對照間的差異基因數目和韋恩圖; B和C: 各處理時期2個材料間的差異基因數目和韋恩圖; E: 核心差異表達基因KEGG富集分析; F: 耐旱相關候選基因表達水平熱圖; No11: 耐旱材料; No28: 干旱敏感材料; CK: 正常灌溉; DS24h: 干旱處理24 h; DS36h: 干旱處理36 h; DS48h: 干旱處理48 h。

A and D: the number of differentially expressed genes and Venn diagram between each time under drought stress and well water in two materials; B and C: the number of differentially expressed genes and Venn diagrams between two materials; E: KEGG enrichment pathway of core differentially expressed genes; F: heatmap of candidate genes related to drought tolerance; No11: drought resistance material; No28: drought sensitivity material; CK: well water; DS24h: 24 h under drought stress; DS36h: 36 h under drought stress; DS48h: 48 h under drought stress.

圖6 WGCNA模塊聚類樹狀圖和各模塊與性狀的相關性

縮寫同表1。Abbreviations are the same as those given in Table 1.

對MEdeeppink1、MEfirebrick2、MEindianred4和MEplum模塊差異表達基因進行KEGG富集分析, MEdeeppink1模塊顯著富集到戊糖和葡萄糖醛酸相互轉化、光合作用-天線蛋白、檸檬烯和松萜降解、氨基酸生物合成、淀粉和蔗糖代謝以及角質、木栓素和蠟質的生物合成途徑(圖7-A); MEfirebrick2模塊顯著富集到自噬調節、脂肪酸降解、脂肪酸代謝及纈氨酸、亮氨酸和異亮氨酸降解途徑(圖7-B); MEindianred4模塊顯著富集到光合作用、光合作用-天線蛋白、卟啉與葉綠素代謝、植物激素信號轉導、核糖體和光合生物固碳作用途徑(圖7-C); MEplum模塊顯著富集到晝夜節律和萜類骨架生物合成途徑(圖7-D)。

4個模塊所包含的轉錄因子編碼基因和顯著富集的KEGG通路基因中, 有23個基因位于QTL定位鑒定的候選區間, 可認定為耐旱相關候選基因(圖7-E)。其中,()、()、()、()和()在干旱脅迫下較正常灌溉下調表達, 且在干旱敏感材料中表達水平更高, 其余基因隨干旱脅迫下調表達。()、()、()、()、()和()表達水平在正常灌溉下, 耐旱材料顯著高于干旱敏感材料;()、()、()、()表達水平在正常灌溉下, 耐旱材料低于干旱敏感材料;()、()、()、()、()、()、()和()表達水平在正常灌溉和干旱脅迫處理下, 耐旱材料低于干旱敏感材料, 且正常灌溉和干旱脅迫48 h最為明顯(圖7-E)。

綜上所述, 將3182個核心差異表達基因以及MEdeeppink1、MEfirebrick2、MEindianred4和MEplum模塊差異表達基因視為耐旱材料No11與干旱敏感材料No28響應干旱脅迫的主要差異表達基因, 顯著富集到光合作用、脂肪酸代謝、氨基酸代謝、植物激素信號轉導、核糖體、晝夜節律及角質、木栓素和蠟質的生物合成等相關途徑。基于3182個核心差異表達基因和WGCNA分析, 分別鑒定到8個和23個候選基因, 其中位于和重疊QTL區間的()、位于和重疊QTL區間的()和位于區間的()被重復檢測到(表4)。

圖7 各性狀顯著相關模塊的KEGG富集分析和耐旱相關候選基因表達量熱圖

A: MEdeeppink1模塊; B: MEfirebrick2模塊; C: MEindianred4模塊; D: MEplum模塊; E: 熱圖。縮寫同圖5。

A: MEdeeppink1; B: MEfirebrick2; C: MEindianred4; D: MEplum; E: heatmap. Abbreviations are the same as those given in Fig. 5.

表4 甘藍型油菜耐旱相關候選基因

(續表4)

3 討論

王瑞霞等[26]研究發現干旱脅迫使芥菜型春油菜地上部干重顯著下降, 丙二醛含量顯著上升。王丹丹[8]和石鵬程[27]對甘藍型油菜進行干旱脅迫處理, 發現地上部鮮重和干重及葉片相對含水量顯著下降,丙二醛和可溶性糖含量顯著上升。謝小玉等[28]研究發現, 隨著干旱脅迫時間延長和脅迫程度增加, 丙二醛和可溶性糖含量上升, 葉片相對含水量下降。本研究中, 干旱脅迫下地上部鮮重、地上部干重和葉片相對含水量顯著低于正常灌溉, 丙二醛和可溶性糖含量顯著高于正常灌溉, 與前人研究結果一致。

本研究通過SSR和SNP構建的甘藍型油菜遺傳連鎖圖譜, 定位到的耐旱相關QTL和候選區間分布于A01、A02、A06、A08、A09、A10、C02、C03、C04、C06和C09染色體。蒙姜宇等[6]和Gad等[9]定位到的甘藍型油菜干旱脅迫相關QTL, 主要分布在A01、A02、A03、A05、A06、A07、A09、A10、C01、C02、C03和C09染色體; 洪雙[29]、吳金鋒[30]、石鵬程[27]和Khanzada等[11]通過GWAS, 在A01、A02、A03、A04、A05、A07、A08、A09、A10、C01、C02、C03、C04、C05、C06、C07和C09染色體檢測到與甘藍型油菜耐旱性狀相關的SNP位點。由此可見, 甘藍型油菜耐旱性是由多基因控制的復雜數量性狀, 候選基因分布于A、C亞組的多條染色體。

Zhu等[31]對兩種鉀水平下耐旱材料和干旱敏感材料進行轉錄組比較分析, 差異基因主要富集到光合作用、谷胱甘肽生物合成、生長素信號轉導、氨基酸生物合成、半胱氨酸和蛋氨酸代謝以及硫代葡萄糖苷生物合成途徑; Fang等[32]對嚴重干旱和正常培養條件下甘藍型油菜差異表達基因進行KEGG富集分析, 植物激素信號轉導、脂肪酸降解、次生代謝產物生物合成和代謝、光合作用、卟啉和葉綠素代謝、碳代謝、氮代謝以及角質、木栓和蠟的生物合成等途徑被顯著富集到; Xiong等[33]研究發現, 干旱脅迫下秦優8號較對照下調的基因富集到核糖體途徑, 上調的基因富集到次生代謝產物生物合成和植物激素信號轉導途徑; 本研究與上述研究結果基本一致, 對耐旱材料No11和干旱敏感材料No28進行轉錄組測序和分析, 主要差異表達基因顯著富集到光合作用、脂肪酸代謝、氨基酸代謝、植物激素信號轉導、核糖體、晝夜節律及角質、木栓素和蠟質的生物合成途徑。

本研究基于QTL定位和轉錄組數據分析, 鑒定到的甘藍型油菜耐旱候選基因涉及多種類型, 包括轉錄因子FLC、bHLH105、TGA4和TEM1編碼基因, 乙烯產生和信號轉導相關的ACO3和ERF003編碼基因, 葉綠素相關蛋白CHLI1、LHCB6和PORC編碼基因, CYP家族蛋白CYP71B22編碼基因, E2亞基E2-OGDH編碼基因以及核糖體相關蛋白編碼基因等。Chen等[34]的研究結果揭示了擬南芥FRIGIDA在干旱脅迫下通過P5CS1途徑增強耐旱性,且依賴于靶標基因FLC; Sun等[35]認為玉米bHLH105可能通過調節抗氧化機制介導的活性氧清除和Mn/Fe相關轉運蛋白的表達提高錳脅迫耐受性, 可用于開發耐旱玉米新品種。Zhong等[36]發現過表達bZIP轉錄因子編碼基因可提高擬南芥耐旱性。擬南芥基因是RAV轉錄因子亞家族成員, 過表達轉基因擬南芥表現出相對電導率和丙二醛含量增加和脯氨酸含量降低, 而突變體植株耐旱性增強[37]。乙烯在非生物脅迫耐受機制中起著關鍵的信號傳導作用, 編輯和的矮牽牛突變體乙烯的產生降低, 對鹽和干旱脅迫更敏感[38]。乙烯響應因子包含單一AP2/ERF結構域, 在植物發育和生物與非生物脅迫響應中發揮著重要作用, 在鹽和干旱脅迫下, 16個剛毛怪柳乙烯響應因子編碼基因均在葉、莖或根中存在差異表達[39]。在小麥中, TaERF3通過激活脅迫相關基因積極調控對鹽和干旱脅迫的響應[40]。光捕獲葉綠素結合蛋白(LHCB)在保護細胞對ABA的信號傳導中發揮積極作用,的過表達可增強氣孔對ABA的敏感性[41]。當暴露于強光下時,過表達植物中單線態氧和丙二醛的產生均最少[42]。突變體表現出部分ABA不敏感[43]。Condori-Apfata等[44]研究表明E2亞基亞型2-OGDH功能缺失使植物生長加快, 暗呼吸降低, 碳水化合物代謝改變。Niu等[45]發現蘋果MdSUVH1是耐旱性正調控因子, Seraj等[46]發現在乳薊干旱脅迫下上調表達, Pandian等[47]認為CYP在植物非生物和生物脅迫響應中發揮著重要作用, Pal等[48]認為核糖體蛋白編碼基因在生物和非生物脅迫下發揮著重要作用。由此可見, 本研究鑒定的甘藍型油菜耐旱相關候選基因具有可靠性, 但還需進一步開展功能驗證。

4 結論

干旱脅迫下, 甘藍型油菜地上部鮮重、地上部干重和葉片相對含水量顯著低于正常灌溉, 丙二醛和可溶性糖含量顯著高于正常灌溉。通過SSR和SNP構建的甘藍型油菜遺傳連鎖圖譜, 在A01、A02、A06、A08、A09、A10、C02、C03、C04、C06和C09等11條染色體定位到耐旱相關QTL和候選區間。對耐旱材料No11和干旱敏感材料No28進行干旱脅迫處理和轉錄組測序分析, 主要差異表達基因顯著富集到光合作用、脂肪酸代謝、氨基酸代謝、植物激素信號轉導、核糖體、晝夜節律及角質、木栓素和蠟質的生物合成等相關途徑。基于QTL定位和轉錄組數據分析, 鑒定到轉錄因子、乙烯產生和信號轉導、葉綠素生物合成與結合、葉綠素氧化還原酶以及核糖體相關蛋白編碼基因等28個耐旱候選基因, 編碼FLC、bHLH105、TGA4、TEM1、ERF003、ACO3、CHLI1、LHCB6和PORC等。這些候選基因可為進一步進行基因功能驗證、揭示甘藍型油菜耐旱機制奠定基礎, 將有助于耐旱新品種的培育。

[1] Batool M, El-Badri A M, Wang Z K, Mohamed I A A, Yang H Y, Ai X Y, Salah A, Hassan M U, Sami R, Kuai J, Wang B, Zhou G S. Rapeseed morpho-physio-biochemical responses to drought stress induced by PEG-6000., 2022, 12: 579.

[2] 周廣生, 王晶, 蒯婕, 汪波. 專輯導讀: 加強大田經濟作物栽培措施與環境/資源配置的互作研究、推動產業高效優質發展. 作物學報, 2021, 47: 1633–1638. Zhou G S, Wang J, Kuai J, Wang B. Editorial: strengthening the research on the interaction between cultivated measures and environment/resource allocation of field economic crops to promote the development of industry with high efficiency and high quality., 2021, 47: 1633–1638 (in Chinese with English abstract).

[3] 寧寧, 莫嬌, 胡冰, 李大雙, 婁洪祥, 王春云, 白晨陽, 蒯婕, 汪波, 王晶, 徐正華, 李曉華, 賈才華, 周廣生. 長江流域不同生態區油菜籽關鍵品質比較研究. 作物學報, 2023, 49: 3315–3327. Ning N, Mo J, Hu B, Li D S, Lou H X, Wang C Y, Xiang C Y, Kuai J, Wang B, Wang J, Xu Z H, Li X H, Jia C H, Zhou G S. Comparative study on the processing quality of winter rape in different ecological zones of the Yangtze River valley., 2023, 49: 3315–3327 (in Chinese with English abstract).

[4] 張佳運, 馬淑梅, 余常兵, 王淑彬, 魏亞鳳, 楊文鈺, 王小春. 長江流域旱地多熟模式水分供需平衡特征與水分生產效益. 作物學報, 2022, 48: 2891–2907. Zhang J Y, Ma S M, Yu C B, Wang S B, Wei Y F, Yang W Y, Wang X C. Characteristics of water supply-demand equilibrium and water production benefits of the dryland multiple cropping patterns in the Yangtze Rive basin., 2022, 48: 2891–2907 (in Chinese with English abstract).

[5] 萬林, 李張開, 李素, 劉麗欣, 馬霓, 張春雷. 外源獨腳金內酯對油菜苗期干旱脅迫的緩解效應. 中國油料作物學報, 2020, 42: 461–471. Wan L, Li Z K, Li S, Liu L X, Ma N, Zhang C L. Alleviation effects of exogenous strigolactone on growth ofL. seedling under drought stress., 2020, 42: 461–471 (in Chinese with English abstract).

[6] 蒙姜宇, 梁光偉, 賀亞軍, 錢偉. 甘藍型油菜耐鹽和耐旱相關性狀的QTL分析. 作物學報, 2021, 47: 462–471. Meng J Y, Liang G W, He Y J, Qian W. QTL mapping of salt and drought tolerance related traits inL., 2021, 47: 462–471 (in Chinese with English abstract).

[7] 李真. 甘藍型油菜苗期耐濕性和抗旱性相關QTL分析. 華中農業大學碩士學位論文, 湖北武漢, 2008. pp 46–48. Li Z. Study on QTL Associated with Waterlogging Tolerance and Drought Resistance during Seedling Stage inL. MS Thesis of Huazhong Agricultural University, Wuhan, Hubei, China, 2008. pp 46–48 (in Chinese with English abstract).

[8] 王丹丹. 甘藍型油菜遺傳圖譜構建及苗期耐旱相關性狀的QTL定位. 西南大學碩士學位論文, 重慶, 2013. pp 32–33. Wang D D. Mapping and QTL Analysis of Genes to Drought Tolerance inL. MS Thesis of Southwest University, Chongqing, China, 2013. pp 32–33 (in Chinese with English abstract).

[9] Mahmoud G, Chao H B, Li H X, Zhao W G, Lu G Y, Li M T. QTL mapping for seed germination response to drought stress in., 2021, 11: 629970.

[10] Shahzad A, Qian M C, Sun B Y, Mahmood U, Li S T, Fan Y H, Chang W, Dai L S, Zhu H, Li J N, Qu C M, Lu K. Genome-wide association study identifies novel loci and candidate genes for drought stress tolerance in rapeseed., 2021, 6: 12–22.

[11] Khanzada H, Wassan G M, He H H, Mason A S, Keerio A A, Khanzada S, Faheem M, Solangi A M, Zhou Q H, Fu D H, Huang Y J, Rasheed A. Differentially evolved drought stress indices determine the genetic variation ofat seedling traits by genome-wide association mapping., 2020, 24: 447–461.

[12] Tan M, Liao F, Hou L T, Wang J, Wei L J, Jian H J, Xu X F, Li J N, Liu L Z. Genome-wide association analysis of seed germination percentage and germination index inL. under salt and drought stresses., 2017, 213: 40.

[13] 王浩. 基于GWAS和轉錄組測序鑒定谷子硒響應相關候選基因. 山西農業大學碩士學位論文, 山西太原, 2022. p 9. Wang H. Identification of Selenium-Responsive Candidate Genes in Foxtail Millet Based on GWAS and Transcriptome Sequencing. MS Thesis of Shanxi Agricultural University, Taiyuan, Shanxi, China, 2022. p 9 (in Chinese with English abstract).

[14] Liu C Q, Zhang X K, Zhang K, An H, Hu K N, Wen J, Shen J X, Ma C Z, Yi B, Tu J X, Fu T D. Comparative analysis of theroot and leaf transcript profiling in response to drought stress., 2015, 16: 18752–18777.

[15] Wang P, Yang C L, Chen H, Song C P, Zhang X, Wang D J. Transcriptomic basis for drought-resistance inL., 2017, 7: 40532.

[16] Zhou H W, Xiao X J, Asjad A, Han D P, Zheng W, Xiao G B, Huang Y J, Zhou Q H. Integration of GWAS and transcriptome analyses to identify SNPs and candidate genes for aluminum tolerance in rapeseed (L.)., 2022, 22: 130.

[17] Jian H J, Zhang A X, Ma J Q, Wang T Y, Yang B, Shuang L S, Liu M, Li J N, Xu X F, Paterson A H, Liu L Z. Joint QTL mapping and transcriptome sequencing analysis reveal candidate flowering time genes inL., 2019, 20: 21.

[18] Guo J, Li C H, Zhang X Q, Li Y X, Zhang D F, Shi Y S, Song Y C, Li Y, Yang D G, Wang T Y. Transcriptome and GWAS analyses reveal candidate gene for seminal root length of maize seedlings under drought stress., 2020, 292: 110380.

[19] Sevanthi A M, Sinha S K, Sureshkumar V, Rani M, Saini M R, Kumari S, Kaushik M, Prakash C, Venkatesh K, Singh G P, Mohapatra T, Mandal P K. Integration of dual stress transcriptomes and major QTL from a pair of genotypes contrasting for drought and chronic nitrogen starvation identifies key stress responsive genes in rice., 2021, 14: 49.

[20] 荊蓉蓉. 甘藍型油菜苗期耐濕相關性狀的全基因組關聯分析. 西南大學碩士學位論文, 重慶, 2017. p 16. Jing R R. Genome-wide Assocaiton Mapping of Water Logging Traits inL. MS Thesis of Southwest University, Chongqing, China, 2017. p 16 (in Chinese with English abstract).

[21] Van Ooijen J W. JoinMap 4: Software for the Calculation of Genetic Linkage Maps in Experimental Populations, 2006. Wageningen: Kyazma B V. pp 1–63.

[22] Van Ooijen J W. MapQTL 6: Software for the Mapping of Quantitative Trait Loci in Experimental Populations of Diploid Species, 2009. Wageningen: Kyazma B V. pp 1–64.

[23] Voorrips R E, MapChart: software for the graphical presentation of linkage maps and QTL., 2002, 93: 77–78.

[24] 李星, 楊會, 駱璐, 李華東, 張昆, 張秀榮, 李玉穎, 于海洋, 王天宇, 劉佳琪, 王瑤, 劉風珍, 萬勇善. 栽培種花生單仁重QTL定位分析. 作物學報, 2023, 49: 2160–2170. Li X, Yang H, Luo L, Li H D, Zhang K, Zhang X R, Li Y Y, Yu H Y, Wang T Y, Liu J Q, Wang Y, Liu F Z, Wan Y S. QTL mapping for single-seed weight of cultivated peanut., 2023, 49: 2160–2170 (in Chinese with English abstract).

[25] 李陽陽, 荊蓉蓉, 呂蓉蓉, 石鵬程, 李欣, 王芹, 吳丹, 周清元,李加納, 唐章林. 甘藍型油菜濕害脅迫響應性狀的全基因組關聯分析及候選基因預測. 作物學報, 2019, 45: 1806–1821. Li Y Y, Jing R R, Lyu R R, Shi P C, Li X, Wang Q, Wu D, Zhou Q Y, Li J N, Tang Z L. Genome-wide association analysis and candidate genes prediction of waterlogging-responding traits in., 2019, 45: 1806–1821 (in Chinese with English abstract).

[26] 王瑞霞, 杜海平, 田宏先. 干旱脅迫對不同生育期芥菜型春油菜幼苗生長的影響. 南方農業, 2019, 13(26): 140–143. Wang R X, Du H P, Tian H X. Effects of drought stress on the growth of seedlings at different growth stages in spring., 2019, 13(26): 140–143 (in Chinese).

[27] 石鵬程. 甘藍型油菜苗期干旱及旱后復水相關性狀的全基因組關聯分析. 西南大學碩士學位論文, 重慶, 2019. p 15. Shi P C. Genome-wide Association Mapping for Drought and Rewatering Related Traits at Seedling Stage inL. MS Thesis of Southwest University, Chongqing, China, 2019. p 15 (in Chinese with English abstract).

[28] 謝小玉, 張霞, 張兵. 油菜苗期抗旱性評價及抗旱相關指標變化分析. 中國農業科學, 2013, 46: 476–485. Xie X Y, Zhang X, Zhang B. Evalution of drought resistance and analysis of variation of relevant parameters at seedling stage of rapeseed (L.)., 2013, 46: 476–485 (in Chinese with English abstract).

[29] 洪雙. 全基因組關聯分析挖掘甘藍型油菜耐旱候選基因. 中國農業科學院碩士學位論文, 北京, 2018. pp 31–32. Hong S. Genome-wide Association Study Identifies Candidate Genes for Drought Tolerance in. MS Thesis of Chinese Academy of Agricultural Sciences, Beijing, China, 2018. pp 31–32 (in Chinese with English abstract).

[30] 吳金鋒. 甘藍型油菜SNP與SSR分析及耐旱性狀的全基因組關聯分析. 中國農業科學院博士學位論文, 北京, 2014. pp 27–28. Wu J F. SNP and SSR Analysis and Genome-wide Association Mapping of Drought Tolerance Trait in. PhD Dissertation of Chinese Academy of Agricultural Sciences, Beijing, China, 2014. pp 27–28 (in Chinese with English abstract).

[31] Zhu B, Xu H X, Guo X, Lu J X, Liu X Y, Zhang T. Comparative analysis of drought responsive transcriptome ingenotypes with contrasting drought tolerance under different potassium levels., 2023, 219: 25.

[32] Fang S, Zhao P M, Tan Z D, Peng Y, Xu L T, Jin Y T, Wei F, Guo L, Yao X. Combining physio-biochemical characterization and transcriptome analysis reveal the responses to varying degrees of drought stress inL., 2022, 23: 8555.

[33] Xiong J L, Dai L L, Ma N, Zhang C L. Transcriptome and physiological analyses reveal that AM1 as an ABA-mimicking ligand improves drought resistance in., 2018, 85: 73–90.

[34] Chen Q, Zheng Y, Luo L D, Yang Y P, Hu X Y, Kong X X. Functional FRIGIDA allele enhances drought tolerance by regulating the P5CS1 pathway in., 2018, 495: 1102–1107.

[35] Sun K L, Wang H Y, Xia Z L. The maize bHLH transcription factor bHLH105 confers manganese tolerance in transgenic tobacco., 2019, 280: 97–109.

[36] Zhong L, Chen D D, Min D H, Li W W, Xu Z S, Zhou Y B, Li L C, Chen M, Ma Y Z. AtTGA4, a bZIP transcription factor, confers drought resistance by enhancing nitrate transport and assimilation in., 2015, 457: 433–439.

[37] Luo G Y, Liu A L, Zhou X Y, Zhang X W, Peng Y, Chen X B.TEMPRANILLO1 transcription factor AtTEM1 negatively regulates drought tolerance., 2017, 83: 119–127.

[38] Naing A H, Campol J R, Kang H, Xu J P, Chung M Y, Kim C K. Role of ethylene biosynthesis genes in the regulation of salt stress and drought stress tolerance in., 2022, 13: 844449.

[39] Liu W J, Wang Y C, Gao C Q. The ethylene response factor (ERF) genes fromrespond to salt, drought and ABA treatment., 2013, 28: 317–327.

[40] Rong W, Qi L, Wang A Y, Ye X G, Du L P, Liang H X, Xin Z Y, Zhang Z Y. The ERF transcription factor TaERF3 promotes tolerance to salt and drought stresses in wheat., 2014, 12: 468–479.

[41] Xu Y H, Liu R, Yan L, Liu Z Q, Jiang S C, Shen Y Y, Wang X F, Zhang D P. Light-harvesting chlorophyll/-binding proteins are required for stomatal response to abscisic acid in., 2012, 63: 1095–1106.

[42] Pattanayak G K, Tripathy B C. Overexpression of protochlorophyllide oxidoreductase C regulates oxidative stress in., 2011, 6: e26532.

[43] Tomiyama M, Inoue S I, Tsuzuki T, Soda M, Morimoto S, Okigaki Y, Ohishi T, Mochizuki N, Takahashi K, Kinoshita T. Mg-chelatase I subunit 1 and Mg-protoporphyrin IX methyltransferase affect the stomatal aperture in., 2014, 127: 553–563.

[44] Condori-Apfata J A, Batista-Silva W, Medeiros D B, Vargas J R, Valente L M L, Pérez-Díaz J L, Fernie A R, Araújo W L, Nunes- Nesi A. Downregulation of the E2 subunit of 2-Oxoglutarate dehydrogenase modulates plant growth by impacting carbon- nitrogen metabolism in., 2021, 62: 798–814.

[45] Niu C D, Jiang L J, Cao F G, Liu C, Guo J X, Zhang Z T, Yue Q Y, Hou N, Liu Z Y, Li X W, Tahir M M, He J Q, Li Z X, Li C, Ma F W, Guan Q M. Methylation of a MITE insertion in the MdRFNR1-1 promoter is positively associated with its allelic expression in apple in response to drought stress., 2022, 34: 3983–4006.

[46] Seraj R G M, Tohidfar M, Irani M A, Esmaeilzadeh-Salestani K, Moradian T, Ahmadikhah A, Behnamian M. Metabolomics analysis of milk thistle lipids to identify drought-tolerant genes., 2022, 12: 12827.

[47] Pandian B A, Sathishraj R, Djanaguiraman M, Prasad P V V, Jugulam M. Role of cytochrome P450 enzymes in plant stress response., 2020, 9: 454.

[48] Pal G, Bakade R, Deshpande S, Sureshkumar V, Patil S S, Dawane A, Agarwal S, Niranjan V, PrasannaKumar M K, Vemanna R S. Transcriptomic responses under combined bacterial blight and drought stress in rice reveal potential genes to improve multi-stress tolerance., 2022, 22: 349.

Identification of candidate genes associated with drought tolerance based on QTL and transcriptome sequencing inL.

LI Yang-Yang1,2,3, WU Dan2,3, XU Jun-Hong2,3, CHEN Zhuo-Yong1,2,3, XU Xin-Yuan1,2,3, XU Jin-Pan1,2,3, TANG Zhong-Lin1,2,3, ZHANG Ya-Ru1,2,3, ZHU Li1,2,3, YAN Zhuo-Li1,2,3, ZHOU Qing-Yuan1,2,3, LI Jia-Na1,2,3, LIU Lie-Zhao1,2,3, and TANG Zhang-Lin1,2,3,*

1Integrative Science Center of Germplasm Creation in Western China (Chongqing) Science City, Chongqing 401329, China;2College of Agronomy and Biotechnology, Southwest University, Chongqing 400715, China;3Academy of Agricultural Sciences, Southwest University, Chongqing 400715, China

Drought stress severely limits planting promotion and yield increase inL. Drought tolerance is a complex quantitative trait controlled by multiple genes. Combining QTL and transcriptome is an effective method for identifying candidate genes associated with drought tolerance in. In this study, the seedlings of F2:6and F2:8recombinant inbred lines, constructed by Sanliu’ai (drought sensitivity line) and Kelina-2 (drought tolerance line), were treated with drought stress and well watering at seedling stage. Shoot fresh weight, shoot dry weight, leaf relative water content, malondialdehyde content, and soluble sugar content were measured. The QTL and candidate intervals were identified based on genetic linkage maps, which were constructed using SSR and SNP markers with polymorphism. Subsequently, candidate genes associated with drought tolerance were screened by combining transcriptome sequencing of No11 (drought tolerance material) and No28 (drought sensitivity material). Drought stress decreased shoot fresh weight, shoot dry weight, and leaf relative water content, and increased the contents of malondialdehyde and soluble sugar. QTL and candidate intervals related to drought tolerance were distributed on chromosome A01, A02, A06, A08, A09, A10, C02, C03, C04, C06, and C09. By transcriptome analysis of drought tolerance and sensitivity materials under well water, drought stress for 24, 36, and 48 h, the major different expression genes were enriched in the pathways associated with photosynthesis, fatty acid metabolism, amino acid metabolism, plant hormone signal transduction, ribosome, circadian rhythm and biosynthesis of keratin, cork and wax. A total of 28 candidate genes related to drought tolerance were identified by combining QTL and transcriptome. They coded FLC, bHLH105, TGA4, TEM1, ERF003, ACO3, CHLI1, LHCB6, PORC, etc., which had transcription factor activity, ethylene production and signal transduction, chlorophyll biosynthesis and binding, chlorophyll oxidoreductase and encoding ribosome proteins. These results could provide a basis for revealing drought tolerance mechanism and molecular breeding of drought tolerance variety in.

L.; drought tolerance; QTL; transcriptome; candidate genes

10.3724/SP.J.1006.2024.34144

本研究由國家重點研發計劃項目(2022YFD1201600)和重慶市技術創新與應用發展專項重點項目(cstc2021jscx-cylhX0003)資助。

This study was supported by the National Key Research and Development Program of China (2022YFD1201600) and the Chongqing Technology Innovation and Application Development Key Project (cstc2021jscx-cylhX0003).

唐章林, E-mail: tangzhlin@swu.edu.cn

E-mail: yyli1993swu@163.com

2023-08-24;

2023-10-23;

2023-11-13.

URL: https://link.cnki.net/urlid/11.1809.S.20231110.1316.008

This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

主站蜘蛛池模板: 高清欧美性猛交XXXX黑人猛交 | 日韩欧美国产区| 丁香婷婷在线视频| 91青草视频| 亚洲色欲色欲www在线观看| 日本亚洲成高清一区二区三区| 高清不卡毛片| 亚洲Av激情网五月天| 在线欧美国产| 久久国产精品国产自线拍| 在线播放精品一区二区啪视频| 又爽又大又黄a级毛片在线视频| 无码内射中文字幕岛国片| 欧美亚洲香蕉| 91无码人妻精品一区| 69国产精品视频免费| 久热中文字幕在线观看| 亚洲天天更新| 国产精品55夜色66夜色| 人妻少妇乱子伦精品无码专区毛片| 亚洲中文字幕久久精品无码一区 | 国产亚洲精品97AA片在线播放| 97se亚洲综合在线天天| 国产亚洲日韩av在线| 久久网欧美| 国产精品亚洲va在线观看| a毛片免费在线观看| 国产精品吹潮在线观看中文| 欧美一区二区福利视频| 亚洲国产成人在线| 亚洲无码高清免费视频亚洲| 国产剧情国内精品原创| 超碰免费91| 国产精品亚洲片在线va| 婷婷色狠狠干| 日本免费福利视频| 国产亚洲精品无码专| 91尤物国产尤物福利在线| 成人在线观看不卡| 亚洲成人网在线观看| 久久久久久久97| 日韩免费中文字幕| 亚洲无码视频一区二区三区| 国产欧美日韩精品综合在线| www.国产福利| 四虎影视国产精品| 国产尤物视频网址导航| 伊人精品视频免费在线| 午夜精品区| 欧美中文字幕一区二区三区| 亚洲人成网7777777国产| 亚洲一区网站| 国产经典免费播放视频| 无码日韩精品91超碰| 国产精品女主播| 久久亚洲国产视频| 国产一区二区三区免费观看| 欧洲精品视频在线观看| 成人在线不卡| 国产亚洲视频在线观看| 亚洲国产欧洲精品路线久久| 五月六月伊人狠狠丁香网| 亚洲国产黄色| 亚洲无限乱码一二三四区| 欧美啪啪视频免码| 久久亚洲日本不卡一区二区| 久久性视频| 片在线无码观看| 亚洲男女在线| 高清久久精品亚洲日韩Av| 99免费在线观看视频| 日韩av无码精品专区| 高清欧美性猛交XXXX黑人猛交| 久久成人18免费| 欧洲av毛片| 精品三级在线| 国产激爽大片高清在线观看| 欧洲av毛片| 亚洲区第一页| 国产男人的天堂| 99久久精品免费视频| 国产在线麻豆波多野结衣|