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

    大菱鲆體重和體尺性狀聯(lián)合GWAS分析*

    2021-03-06 03:03:48楊潤(rùn)清
    漁業(yè)科學(xué)進(jìn)展 2021年2期
    關(guān)鍵詞:大菱鲆體尺表型

    高 進(jìn) 楊潤(rùn)清

    大菱鲆體重和體尺性狀聯(lián)合GWAS分析*

    高 進(jìn)1,2楊潤(rùn)清1,2①

    (1. 南京農(nóng)業(yè)大學(xué)無錫漁業(yè)學(xué)院 無錫 214081;2. 中國(guó)水產(chǎn)科學(xué)研究院生物技術(shù)研究中心 北京 100141)

    為了揭示大菱鲆()體重和體尺性狀的分子遺傳機(jī)制,探尋用于改良目標(biāo)性狀的分子標(biāo)記及候選基因,本研究以大菱鲆育種群體為研究對(duì)象,分別測(cè)量其體重、體長(zhǎng)、體寬和尾柄寬性狀的表型值,利用簡(jiǎn)化基因組測(cè)序技術(shù)(2b-RAD)獲得相應(yīng)基因型數(shù)據(jù),進(jìn)行全基因組關(guān)聯(lián)研究(Genome-wide association study, GWAS),篩選與大菱鲆體重和體尺性狀顯著關(guān)聯(lián)的數(shù)量性狀核苷酸(Quantitative trait nucleotides, QTNs)遺傳位點(diǎn)。結(jié)果顯示,以多性狀線性混合模型(mvLMM)對(duì)體重–體長(zhǎng)和體長(zhǎng)–體寬–尾柄寬2個(gè)性狀組合進(jìn)行多性狀GWAS分析,分別檢測(cè)到9個(gè)和2個(gè)一因多效QTNs;以單一性狀線性混合模型(LMM)對(duì)各個(gè)性狀進(jìn)行GWAS分析,在體重性狀中檢測(cè)到4個(gè)與之顯著關(guān)聯(lián)的QTNs,在體長(zhǎng)和體寬性狀中各檢測(cè)到1個(gè)QTN,而在尾柄寬性狀中則沒有檢測(cè)到顯著的遺傳位點(diǎn)。比較2種模型的結(jié)果,發(fā)現(xiàn)mvLMM相較于LMM能夠檢測(cè)到更多QTNs,且檢測(cè)到的QTNs為更具生物學(xué)意義的一因多效QTNs。本研究首次利用mvLMM和LMM對(duì)大菱鲆體重和體尺性狀進(jìn)行聯(lián)合GWAS分析,共篩選到17個(gè)顯著的QTNs,其中,有4個(gè)QTNs被重復(fù)檢測(cè)到。以這些檢測(cè)到的QTNs為探針,在大菱鲆全基因組上找到了距離其最近的12個(gè)候選基因,它們可能是影響大菱鲆體重和體尺性狀的重要候選標(biāo)記和功能基因,本研究為大菱鲆體重和體尺性狀的分子標(biāo)記輔助選育提供了理論素材和參考。

    大菱鲆;體重和體尺性狀;全基因組關(guān)聯(lián)分析;多性狀線性混合模型

    大菱鲆(),俗稱“歐洲比目魚”或“多寶魚”,是中國(guó)北部沿海地區(qū)重要的海水養(yǎng)殖經(jīng)濟(jì)魚類。大菱鲆具有生長(zhǎng)速度快、肉質(zhì)鮮美以及抗逆性強(qiáng)等諸多優(yōu)點(diǎn),其養(yǎng)殖區(qū)域迅速擴(kuò)展、產(chǎn)量產(chǎn)值不斷攀升(雷霽霖等, 2002、2005)。然而,隨著工廠化養(yǎng)殖規(guī)模和集約化程度的擴(kuò)大,特別是種質(zhì)資源的退化,成活率低、病害嚴(yán)重和生產(chǎn)緩慢等問題逐漸顯現(xiàn)(李杰等, 2019; 王嵐等, 2017)。因此,開展大菱鲆相關(guān)經(jīng)濟(jì)性狀的選育研究,培育具有強(qiáng)抗逆性和高生長(zhǎng)性能的新品種對(duì)大菱鲆產(chǎn)業(yè)發(fā)展具有重要意義。傳統(tǒng)的育種方法是基于目標(biāo)性狀的表型值來開展選育工作(Hulata, 2001),對(duì)于高生長(zhǎng)性能品種選育,通常挑選符合標(biāo)準(zhǔn)的個(gè)體進(jìn)行群體混養(yǎng)或是建立多個(gè)家系并選擇優(yōu)良的家系作為親本用于進(jìn)一步選育。然而,分子生物學(xué)和基因組學(xué)的快速發(fā)展則促使魚類育種研究由傳統(tǒng)選育和雜交育種向基于基因組信息的分子育種轉(zhuǎn)變(Xu, 2015),如以基因定位和數(shù)量性狀位點(diǎn)(Quantitative trait locus, QTL)作圖為基礎(chǔ)的分子標(biāo)記輔助選擇(Marker-assisted selection, MAS)育種(Wang, 2015; 劉曉菲等, 2019)。利用目標(biāo)性狀與分子標(biāo)記間的顯著關(guān)聯(lián)性,依據(jù)標(biāo)記對(duì)后代個(gè)體進(jìn)行全基因組選擇或是綜合標(biāo)記與表型信息構(gòu)建選擇指數(shù)篩選個(gè)體,加速育種過程,提高選育效率。

    隨著測(cè)序成本的下降與高通量基因分型技術(shù)的迅猛發(fā)展,基于海量單核苷酸多態(tài)性(Single-nucleotide polymorphisms, SNP)標(biāo)記的全基因組關(guān)聯(lián)分析(Genome-wide association studies, GWAS)已取代利用稀疏遺傳標(biāo)記的QTL定位,成為鑒定與目標(biāo)性狀顯著關(guān)聯(lián)的遺傳變異的主要方法。早期GWAS研究主要針對(duì)單一性狀表型(Yu, 2006; Burton, 2007),而在進(jìn)行復(fù)雜性狀的遺傳解析且記錄有多個(gè)相關(guān)性狀表型值時(shí),考慮多個(gè)性狀聯(lián)合信息的多性狀GWAS策略通常要優(yōu)于逐個(gè)性狀與遺傳位點(diǎn)間的關(guān)聯(lián)分析(Kim, 2009; Bolormaa, 2010)。相較單一性狀GWAS分析,表型具有潛在相關(guān)性的多性狀GWAS分析更具優(yōu)勢(shì),主要是因?yàn)椋寒?dāng)不同性狀間存在遺傳相關(guān)時(shí),因考慮了單一性狀分析中忽略的性狀間協(xié)方差,多性狀分析提高了檢驗(yàn)功效(Stephens, 2013)和參數(shù)估計(jì)的精度(Zhu, 2009);在多性狀GWAS分析過程中,數(shù)量性狀核苷酸(Quantitative trait nucleotides, QTNs)的統(tǒng)計(jì)推斷只進(jìn)行一次統(tǒng)計(jì)檢驗(yàn),與逐個(gè)性狀單獨(dú)分析相比,降低了多重檢驗(yàn)造成的誤差(Klei, 2008);存在基因多效性時(shí),單個(gè)遺傳變異與多個(gè)性狀相關(guān)聯(lián),使得多性狀GWAS分析更具生物學(xué)意義(Chavali, 2010);Porter等(2017)綜合比較了當(dāng)下流行的多性狀GWAS分析方法,發(fā)現(xiàn)大多數(shù)現(xiàn)存的多性狀GWAS分析方法具有明顯相似的統(tǒng)計(jì)檢驗(yàn)功效,與單一性狀GWAS相比,能夠大幅增加顯著遺傳變異的檢測(cè)效率。

    本研究采用實(shí)驗(yàn)成本較低的簡(jiǎn)化基因組(2b- RAD)測(cè)序技術(shù)進(jìn)行SNP分型,基于多性狀線性混合模型對(duì)體重–體長(zhǎng)和體長(zhǎng)–體寬–尾柄寬共2種多性狀相關(guān)表型組合進(jìn)行全基因組關(guān)聯(lián)分析,關(guān)聯(lián)定位控制這些性狀的一因多效QTNs,為大菱鲆的體重和體尺性狀改良提供理論基礎(chǔ)。此外,對(duì)逐個(gè)性狀進(jìn)行了單一性狀線性混合模型GWAS分析,并將其結(jié)果與多性狀GWAS分析結(jié)果進(jìn)行比較。

    1 材料與方法

    1.1 實(shí)驗(yàn)群體及表型測(cè)量

    實(shí)驗(yàn)用魚共585條,取自29個(gè)全同胞家系組成的大菱鲆育種群體,4月齡時(shí),對(duì)群體中所有大菱鲆個(gè)體注射電子標(biāo)記,并提取鰭條組織DNA用于后續(xù)的2b-RAD測(cè)序。隨后,將標(biāo)記后的大菱鲆隨機(jī)混養(yǎng)在2個(gè)6 m×6 m×1.5 m具有循環(huán)水養(yǎng)殖系統(tǒng)的水泥池中,控制水溫在5℃~24℃,每天定時(shí)投喂商用餌料2次至飽食。實(shí)驗(yàn)用魚按孵化日期至表型測(cè)量時(shí)的生長(zhǎng)天數(shù)為275~1001 d不等。對(duì)實(shí)驗(yàn)用魚進(jìn)行11次不定期表型測(cè)量,每次測(cè)量前利用50 mg/ml MS-222魚用安定劑麻醉待測(cè)大菱鲆個(gè)體,避免應(yīng)激反應(yīng)造成的魚體損傷。利用電子秤稱量每個(gè)個(gè)體的體重(Body weight, g, BM)性狀,同時(shí)在統(tǒng)一的參考標(biāo)尺下,用數(shù)碼相機(jī)由固定高度向下垂直拍攝相應(yīng)個(gè)體的體尺性狀。根據(jù)拍攝圖形,利用ImageJ軟件標(biāo)定每條大菱鲆的體長(zhǎng)(Body length, cm, BL)、體寬(Body width, cm, BW)和尾柄寬(Caudal peduncle width, cm, CPW)共3個(gè)體尺性狀的表型值。采用簡(jiǎn)化基因組2b-RAD高通量標(biāo)記分型技術(shù)對(duì)具有表型測(cè)量記錄的個(gè)體進(jìn)行SNP分型,共獲得30049個(gè)多態(tài)SNP分子標(biāo)記,參考大菱鲆基因組(ASM318616v1)建立多態(tài)標(biāo)記物理圖譜。挑選生長(zhǎng)周期為473 d左右(前后皆不超過5 d)大菱鲆個(gè)體的表型觀測(cè)值作為多性狀GWAS分析的表型值,使用PLINK v1.9 (http://www.cog- genomics.org/plink2)對(duì)相應(yīng)樣本個(gè)體的基因型數(shù)據(jù)進(jìn)行質(zhì)量控制。剔除低于90%最小檢出率的個(gè)體以及最小檢出頻率低于95%、最小哈代溫伯格平衡為1.0×10–6、最小等位基因頻率小于3%和方差變異大于0.05的SNPs,最終得到441個(gè)樣本的23988個(gè)SNP標(biāo)記用于全基因組關(guān)聯(lián)分析。

    1.2 關(guān)聯(lián)統(tǒng)計(jì)分析

    本研究使用GEMMA軟件(Zhou, 2014)中的多性狀線性混合模型(Multivariate linear mixed model, mvLMM)進(jìn)行多個(gè)性狀表型值和SNPs標(biāo)記的全基因組關(guān)聯(lián)分析,所用模型為:

    =+++

    式中,為×維表型矩陣,是樣本個(gè)數(shù),是分析性狀個(gè)數(shù);為×維協(xié)變量(非遺傳固定效應(yīng),如性別、年齡等)矩陣,為維相應(yīng)系數(shù)行向量,是包含截距項(xiàng)在內(nèi)的協(xié)變量個(gè)數(shù);為維當(dāng)前檢驗(yàn)SNP的基因型指示變量列向量,為維當(dāng)前檢驗(yàn)SNP的加性遺傳效應(yīng)行向量;為×維不包括當(dāng)前檢驗(yàn)SNP的剩余多基因效應(yīng)矩陣,~(0, V,是由全基因組SNP標(biāo)記構(gòu)建的基因組親緣關(guān)系矩陣,V是×維剩余多基因方差–協(xié)方差矩陣;為×維誤差項(xiàng)矩陣,~(0,VI),V是×維誤差方差–協(xié)方差矩陣,是單位矩陣。

    最終,全基因組上每個(gè)SNP都能得到1個(gè)一因多效Wald統(tǒng)計(jì)量和相對(duì)應(yīng)的統(tǒng)計(jì)概率值。此外,還利用因子譜分解線性混合模型(Factored spectrally transformed linear mixed model, FaST-LMM) (Lippert, 2011)逐個(gè)性狀進(jìn)行單一性狀的GWAS分析,檢測(cè)與各個(gè)性狀相關(guān)聯(lián)的QTNs。使用R語(yǔ)言繪制曼哈頓和Quantile-Quantile (QQ)圖,同時(shí),統(tǒng)計(jì)用于判別群體分層影響大小的膨脹系數(shù)或稱基因組控制(Genomic control, GC)值,在實(shí)際研究中GC值被定義為所有SNP標(biāo)記的卡方統(tǒng)計(jì)量均值(Price, 2010)。

    1.3 遺傳參數(shù)估計(jì)和基因注釋

    將質(zhì)量控制后的基因型數(shù)據(jù)用PLINK軟件(Chang, 2015)處理為分析所需的格式,再利用單性狀約束最大法估計(jì)法(Restricted maximum likelihood, REML)逐個(gè)估計(jì)各性狀的遺傳參數(shù)。查找關(guān)聯(lián)分析得到的顯著SNP位點(diǎn)在全基因組上所處的物理位置,于大菱鲆全基因組(GCA_003186165.1)上選擇距離其最近的候選基因進(jìn)行注釋分析。

    2 結(jié)果

    2.1 表型值的描述性統(tǒng)計(jì)

    選取生長(zhǎng)天數(shù)為473 d左右的441條大菱鲆的體重、體長(zhǎng)、體寬和尾柄寬4個(gè)性狀進(jìn)行聯(lián)合GWAS分析。4個(gè)性狀的表型頻數(shù)分布見圖1。從圖1可以看出,各個(gè)性狀的表型值基本都符合正態(tài)分布,具有一定的可靠性,適合后續(xù)的GWAS分析。各性狀的原始表型的統(tǒng)計(jì)分析和遺傳力估計(jì)值見表1。此外,還對(duì)原始以及校正了性別、池子和測(cè)定日期等非遺傳固定效應(yīng)后的表型值進(jìn)行了性狀間的相關(guān)分析,發(fā)現(xiàn)各性狀表型值間具有較強(qiáng)的相關(guān)性(表2)。

    2.2 GWAS分析結(jié)果

    根據(jù)質(zhì)量控制結(jié)果,441條大菱鲆的23988個(gè)SNPs被用于體重和體尺性狀的聯(lián)合GWAS分析,分別進(jìn)行了BM-BL和BL-BW-CPW兩個(gè)性狀組合的多性狀線性混合模型以及逐個(gè)性狀的單一性狀線性混合模型GWAS分析。對(duì)BM-BL兩個(gè)性狀的多性狀GWAS分析后(圖2和表3)發(fā)現(xiàn),超過基因組顯著水平(5% Bonferroni校正閾值,即2.084×10–6)的QTN位點(diǎn)共9個(gè),其中,各有1個(gè)顯著QTN位于1、10和20號(hào)染色體上,即SNP_1_31495825、SNP_10_ 6893888和SNP_20_18773114,分別有2個(gè)顯著QTNs位于3 (SNP_3_21601589與SNP_3_24007357)、5(SNP_5_25635891與SNP_5_26888833)和22(SNP_22_ 6190492與SNP_22_6380029)號(hào)染色體上。而在相同顯著水準(zhǔn)下,BL-BW-CPW 3個(gè)體尺性狀的GWAS分析則在10號(hào)染色體上檢測(cè)到2個(gè)顯著的QTN位點(diǎn)SNP_10_2724296和SNP_10_6893888。2組多性狀GWAS分析的GC值分別為BM-BL(1.035)和BL- BW-CPW(1.044)。

    圖1 4個(gè)性狀的頻數(shù)分布

    表1 大菱鲆體重和體尺性狀的描述性統(tǒng)計(jì)

    Tab.1 Descriptive statistics of body weight and morphological traits in turbot

    表2 大菱鲆體重和體尺性狀相關(guān)程度

    Tab.2 Degree of correlation for body mass and morphological traits in turbot

    注:下三角中的值為性狀原始表型間相關(guān)程度,上三角中的值為校正了非遺傳固定效應(yīng)后表型間的相關(guān)程度

    Note: Values in lower-triangular matrix represent degrees of correlation among original phenotypes, and values in up-triangular matrix mean degrees of correlation among phenotypes which were adjusted by non-genetic fixed effects

    對(duì)體重和體尺的逐個(gè)性狀單一性狀GWAS分析結(jié)果(圖3和表3)顯示,在BM性狀中,于3和5號(hào)染色體上各檢測(cè)到1個(gè)顯著的QTN位點(diǎn),即SNP_3_21601589和SNP_5_26888833,10號(hào)染色體上檢測(cè)到2個(gè)顯著的QTN位點(diǎn)SNP_10_2724296和SNP_10_6893888。BL和BW性狀分別在62和8號(hào)染色體上檢測(cè)到1個(gè)顯著的QTN位點(diǎn)SNP_6_ 6047375和SNP_8_3870447。而在CPW性狀中則沒有檢測(cè)到顯著的遺傳變異位點(diǎn)。各性狀單一性狀GWAS分析的GC值分別為1.026 (BM)、1.021 (BL)、1.034 (BW)和1.031 (CPW)。

    圖2 大菱鲆體重與體尺多性狀GWAS曼哈頓和QQ圖

    圖3 大菱鲆體重與體尺單一性狀GWAS曼哈頓和QQ圖

    表3 大菱鲆體重和體尺性狀顯著關(guān)聯(lián)的QTNs信息

    Tab.3 Information of QTNs significantly associated with body mass and morphological traits in turbot

    注: *表示未在大菱鲆全基因組中進(jìn)行功能注釋的新基因

    Note: * represents novel gene which has not found any gene function annotation in genome of

    2.3 候選基因

    以大菱鲆多性狀與單一性狀GWAS獲得的顯著SNP位點(diǎn)為探針,根據(jù)其在大菱鲆全基因組的位置,向上下游尋找距離最近的候選基因。在本研究檢測(cè)到的17個(gè)QTNs附近共找到12個(gè)候選基因(表3),分別為GRIK2 (Glutamate receptor, ionotropic, kainate 2)、cyp21a2 (Cytochrome P450, family 21, subfamily A, polypeptide 2)、tshz3b (Teashirt zinc finger homeobox 3b)、adamts17 (ADAM metallopeptidase withthrombo- spondin type 1 motif, 17)、SYT8 (Synaptotagmin 8)、myt1la (Myelin transcription factor 1-like, a)、psmd4a (Proteasome 26S subunit, non-ATPase 4a)、grm4 (Glutamate receptor, metabotropic 4)、BFAR (Bifunctional apoptosis regulator)和3個(gè)新基因。對(duì)這些候選基因進(jìn)行GO注釋分析發(fā)現(xiàn),細(xì)胞組分主要涉及核膜、細(xì)胞質(zhì)膜及細(xì)胞核等;分子功能主要涉及離子通道活性、激素活性、氧化還原酶活性、肽酶活性、核苷酸結(jié)合、蛋白結(jié)合及G蛋白偶聯(lián)受體活性等;而生物學(xué)進(jìn)程則主要涉及離子轉(zhuǎn)運(yùn)、糖皮質(zhì)激素生物合成過程、氧化還原過程、細(xì)胞信號(hào)轉(zhuǎn)導(dǎo)、基因表達(dá)調(diào)控、蛋白酶解、胞外分泌、轉(zhuǎn)錄調(diào)控和G蛋白偶聯(lián)受體信號(hào)通路等。其中,3個(gè)新基因在大菱鲆基因組中未有相關(guān)注釋信息。

    3 討論

    GWAS作為一種定位影響重要經(jīng)濟(jì)性狀分子標(biāo)記的有效方法,是以全基因組SNPs直接關(guān)聯(lián)復(fù)雜性狀表型,廣泛應(yīng)用于動(dòng)植物分子標(biāo)記開發(fā)和輔助選擇育種等方面(Santana, 2015; He, 2017)。隨著諸多水產(chǎn)動(dòng)物基因組測(cè)序工作相繼完成以及測(cè)序成本的下降,水產(chǎn)動(dòng)物生長(zhǎng)和抗病等復(fù)雜性狀的GWAS分析多有研究報(bào)道。Gutierrez等(2015)對(duì)大西洋鮭()在性成熟時(shí)的生長(zhǎng)速率性狀進(jìn)行GWAS分析,檢測(cè)到多個(gè)在大西洋鮭代謝和生長(zhǎng)發(fā)育過程中發(fā)揮重要作用的候選基因。Zhou等(2019)利用GWAS揭示了半滑舌鰨()抗病的遺傳機(jī)制,共定位了與抗病性狀顯著關(guān)聯(lián)的33個(gè)SNP位點(diǎn),結(jié)合基因表達(dá)和甲基化分析,檢測(cè)到數(shù)個(gè)影響半滑舌鰨抗病力的候選基因。Jiang等(2019)則對(duì)羅非魚()的耐鹽性狀進(jìn)行GWAS和QTL定位研究,為羅非魚耐鹽性遺傳機(jī)理揭示及進(jìn)一步的功能分析奠定了基礎(chǔ)。目前,大菱鲆生長(zhǎng)和體尺相關(guān)性狀的GWAS分析鮮有應(yīng)用,但QTL定位研究已在大菱鲆生長(zhǎng)(Enrique, 2011)和氣單胞菌()耐藥性(Silvia, 2011)等相關(guān)性狀中開展。

    現(xiàn)有的絕大多數(shù)GWAS研究都是關(guān)聯(lián)全基因組上SNP位點(diǎn)與單個(gè)性狀表型值,即使當(dāng)多個(gè)相關(guān)表型性狀存在時(shí),往往也會(huì)選擇逐個(gè)性狀分別進(jìn)行單一性狀GWAS分析(Stephens, 2013)。在多個(gè)相關(guān)表型性狀條件下,單一性狀關(guān)聯(lián)定位方法忽略了QTNs對(duì)多個(gè)表型性狀的共同影響。相較于單一性狀關(guān)聯(lián)分析方法,多性狀GWAS方法具有更高的統(tǒng)計(jì)功效和參數(shù)估計(jì)及定位精確度。通過比較大菱鲆體重和體尺多性狀GWAS與各個(gè)性狀逐個(gè)進(jìn)行單一性狀GWAS的定位結(jié)果發(fā)現(xiàn),在相同的顯著水準(zhǔn)下,前者檢測(cè)到的QTN數(shù)多于后者(表3)。在BM-BL性狀組合的多性狀GWAS分析中,共檢測(cè)到9個(gè)QTNs,而在相對(duì)應(yīng)的BM和BL單一性狀GWAS分析中,BM中檢測(cè)到4個(gè)QTNs,BL中僅定位到1個(gè)QTN,且BM中檢測(cè)到的4個(gè)QTNs中有3個(gè)與多性狀GWAS分析得到的QTNs重復(fù)。同樣,BL-BW-CPW性狀組合多性狀GWAS分析檢測(cè)到的QTNs也多于3個(gè)性狀各自進(jìn)行單一性狀GWAS分析得到的QTNs。除了能夠檢測(cè)到更多的QTNs,由于多性狀GWAS分析檢測(cè)到的QTN是一因多效QTN,相較于單一性狀GWAS分析檢測(cè)到的QTN在生物學(xué)和實(shí)際應(yīng)用層面更有意義。

    在大菱鲆全基因組上尋找距離每個(gè)QTN最近的基因,共找到12個(gè)候選基因,其中,有9個(gè)已知的功能基因及3個(gè)新基因。這些影響大菱鲆體重和體尺性狀的候選基因在全基因組分布較為分散,沒有明顯集中區(qū)域。本研究為首次對(duì)大菱鲆體重和體尺性狀進(jìn)行全基因組關(guān)聯(lián)分析,檢測(cè)到在大菱鲆相關(guān)研究中尚未見報(bào)道的顯著的已知功能基因。相關(guān)研究表明,GRIK2為離子型受體,它偶聯(lián)離子通道并形成受體通道復(fù)合物,且與細(xì)胞缺血缺氧損傷有關(guān)(張冬梅等, 2010)。Eachus等(2017)和Weger等(2018)的研究均表明,cyp21a2在斑馬魚()幼魚的糖皮質(zhì)激素生物合成中發(fā)揮著重要作用。Erickson等(2011)以斑馬魚為實(shí)驗(yàn)動(dòng)物模型,發(fā)現(xiàn)tshz3b可能調(diào)控后腦中的Hox功能,而Hox功能則被認(rèn)為與形態(tài)發(fā)生和器官形成有關(guān)聯(lián)(Gair, 2003)。由此可見,該基因可能參與大菱鲆體形態(tài)發(fā)育過程。Myt1la可能在有絲分裂結(jié)束階段發(fā)揮重要作用(Nakajima, 2008),BFAR則能與p75NTR的蛋白相互作用并抑制p75NTR信號(hào)轉(zhuǎn)導(dǎo),而p75NTR可與高親和力受體TrkA協(xié)同作用促進(jìn)細(xì)胞增殖或與細(xì)胞內(nèi)配體結(jié)合誘導(dǎo)細(xì)胞凋亡 (李紅梅等, 2011)。本研究找到的候選基因可能都是大菱鲆體重和體尺性狀的重要候選功能基因,它們?cè)诖罅怫疑L(zhǎng)發(fā)育過程中的影響有待進(jìn)一步的功能實(shí)驗(yàn)驗(yàn)證。

    Allison DB, Thiel B, St Jean P,. Multiple phenotype modeling in gene-mapping studies of quantitative traits: Power advantages. American Journal of Human Genetics, 1998, 63(4): 1190–1201

    Burton PR, Clayton DG, Cardon LR,. Genome-wide association study of 14,000 cases of seven common diseases and 3,000 shared controls. Nature, 2007, 447(7145): 661–678

    Chang CC, Chow CC, Tellier LC,. Second-generation PLINK: Rising to the challenge of larger and richer datasets. GigaScience, 2015, 4(7): 1–16

    Chavali S, Barrenas F, Kanduri K,Network properties of human disease genes with pleiotropic effects. BMC Systems Biology, 2010, 4(78): 1–11

    Eachus H, Zaucker A, Oakes JA,. Genetic disruption of 21-hydroxylase in zebrafish causes interrenal hyperplasia. Endocrinology, 2017, 158(12): 4165–4173

    Enrique SM, Alex C, Miguel AT,. Detection of growth- related QTL in turbot (). BMC Genomics, 2011, 12(1): 473

    Erickson T, Pillay LM, Waskiewicz AJ. Zebrafish Tshz3b negatively regulates hox function in the developing hindbrain. Genesis, 2011, 49(9): 725–742

    Gair HJC, Lovegrove B. Beyond homeosis-HOX function in morphogenesis and organogenesis. Differentiation, 2003, 71(8): 461–476

    在施肥之前,應(yīng)該結(jié)合測(cè)土配方施肥結(jié)果以及不同整地模式,適當(dāng)施入農(nóng)家肥、有機(jī)肥和化肥。施肥要堅(jiān)持以有機(jī)肥為主、化肥為輔的施肥原則。此外,確保泥水能夠有效分層,沉淀時(shí)間長(zhǎng)短應(yīng)該根據(jù)當(dāng)?shù)厮咎飰K的實(shí)際情況綜合確定,一般沉淀時(shí)間維持在1-3天。沉淀好的秧田,應(yīng)該確保沉淀不板結(jié),水清澈不渾濁,然后才能栽插秧苗。機(jī)械化秧苗殘?jiān)饕捎帽∷疁\灌,大泥腳深度維持在30 cm以上,水層控制在1-3 cm。

    Gutierrez AP, José MY, Fukui S,. Genome-wide association study (GWAS) for growth rate and age at sexual maturation in Atlantic salmon (). PLoS One, 2015, 10(3): e0119730

    He JB, Meng S, Zhao TJ,. An innovative procedure of genome-wide association analysis fits studies on germplasm population and plant breeding. Theoretical and Applied Genetics, 2017, 130(11): 2327–2343

    Hulata G. Genetic manipulations in aquaculture: A review of stock improvement by classical and modern technologies. Genetica, 2001, 111(1–3): 155–173

    Jiang DL, Gu XH, Li BJ,. Identifying a long QTL cluster across chrLG18 associated with salt tolerance in Tilapia using GWAS and QTL-seq. Marine Biotechnology, 2019, 21(2): 250–261

    Kim S, Sohn KA, Xing EP. A multivariate regression approach to association analysis of a quantitative trait network. Bioinformatics, 2009, 25(12): 204–212

    Klei L, Luca D, Devlin B,. Pleiotropy and principal components of heritability combine to increase power for association analysis. Genetic Epidemiology, 2008, 32(1): 9–19

    Lei JL, Ma AJ, Chen C,. The present status and sustainable development of turbot (L.) culture in China. Engineering Sciences, 2005(5): 34–38 [雷霽霖, 馬愛軍, 陳超, 等. 大菱鲆(L.)養(yǎng)殖現(xiàn)狀與可持續(xù)發(fā)展. 中國(guó)工程科學(xué), 2005(5): 34–38]

    Lei JL, Men Q, Wang YG,. The factory farming mode of Dalingyu “greenhouse shed + deep well seawater”. Marine Fisheries Research, 2002, 23(4): 1–7 [雷霽霖, 門強(qiáng), 王印庚, 等. 大菱鲆“溫室大棚+深井海水”工廠化養(yǎng)殖模式. 海洋水產(chǎn)研究, 2002, 23(4): 1–7]

    Li HM, Shi HL, Huo KK. p75NTR signal transduction suppressed by BFAR and p75NTR interactions. Scientia Sinica Vitae, 2011, 41(12): 1140–1147 [李紅梅, 施慧莉, 霍克克. BFAR與p75NTR的蛋白相互作用抑制p75NTR信號(hào)轉(zhuǎn)導(dǎo). 中國(guó)科學(xué): 生命科學(xué), 2011, 41(12): 1140–1147]

    Li J, Liu YK, Bai L,. Isolation and identification ofassociated with splenic and renal granuloma disease of cultured turbot (). Progress in Fishery Sciences, 2019, 40(5): 195–199 [李杰, 劉耀寬, 白露, 等. 大菱鲆脾腎結(jié)節(jié)病病原菌的分離和鑒定. 漁業(yè)科學(xué)進(jìn)展, 2019, 40(5): 195–199]

    Lippert C, Listgarten J, Liu Y,. FaST linear mixed models for genome-wide association studies. Nature Methods, 2011, 8(10): 833–835

    Liu XF, Ma AJ, Huang ZH,. Expression characteristics analysis of major QTL candidate genes in response to high temperature stress in turbot (). Journal of Fisheries of China, 2019, 43(6): 1407–1409 [劉曉菲, 馬愛軍, 黃智慧, 等. 大菱鲆高溫脅迫應(yīng)答主效QTL候選基因的表達(dá)特性分析. 水產(chǎn)學(xué)報(bào), 2019, 43(6): 1407–1409]

    Nakajima H, Yonemura S, Murata M,. Myt1 protein kinase is essential for Golgi and ER assembly during mitotic exit. Journal of Cell Biology, 2008, 181(1): 89–103

    Porter HF, O'Reilly PF. Multivariate simulation framework reveals performance of multi-trait GWAS methods. Scientific Reports, 2017, 7: 38837

    Price AL, Zaitlen NA, Reich D,. New approaches to population stratification in genome-wide association studies. Nature Reviews Genetics, 2010, 11(7): 459–463

    Santana MHA, Ventura RV, Utsunomiya YT,. A genomewide association mapping study using ultrasound-scanned information identifies potential genomic regions and candidate genes affecting carcass traits in Nellore cattle. Journal of Animal Breeding and Genetics, 2015, 132(6): 420–427

    Silvia TRR, Toro MA, Bouza C,. QTL detection forresistance related traits in turbot (). BMC Genomics, 2011, 12(1): 541

    Stephens M. A unified framework for association analysis with multiple related phenotypes. PLoS One, 2013, 8(7): e65245.

    Wang L, Wang YG, Zhang Z,. Diversity and drug resistance of bacterial pathogens isolated from bacterial ascetic disease in cultured turbot. Progress in Fishery Sciences, 2017, 38(4): 17–24 [王嵐, 王印庚, 張正, 等. 養(yǎng)殖大菱鲆()腹水病的病原多樣性及其耐藥性分析. 漁業(yè)科學(xué)進(jìn)展, 2017, 38(4): 17–24]

    Wang WJ, Hu YL, Ma Y,. High-density genetic linkage mapping in turbot (L.) based on SNP markers and major sex- and growth-related regions detection. PLoS One, 2015, 10(3): e0120410

    Weger M, Diotel N, Weger BD,. Expression and activity profiling of the steroidogenic enzymes of glucocorticoid biosynthesis and the FDX1 co-factors in zebrafish. Journal of Neuroendocrinology, 2018, 30(4): e12586

    Xu K, Duan W, Xiao J,. Development and application of biological technologies in fish genetic breeding. Science China, Life Sciences, 2015, 58(2): 187–201

    Yu JM, Pressoir G, Briggs WH,. A unified mixed-model method for association mapping that accounts for multiple levels of relatedness. Nature Genetics, 2006, 38(2): 203–208

    Zhang DM, Wu AM, Lou LX,. Alteration of ion channel gene expression profile in rat model of post-myocardial infarction heart failure. Chinese Journal of Integrated Traditional and West Medicine, 2010, 30(1): 53–57 [張冬梅, 吳愛明, 婁利霞, 等. 心肌梗死后心力衰竭模型大鼠離子通道基因表達(dá)譜的變化. 中國(guó)中西醫(yī)結(jié)合雜志, 2010, 30(1): 53–57]

    Zhou Q, Su ZC, Li YZ,. Genome-wide association mapping and gene expression analyses reveal genetic mechanisms of disease resistance variations in. Frontiers in Genetics, 2019, 10: 1167

    Zhou X, Stephens M. Efficient multivariate linear mixed model algorithms for genome-wide association studies. Nature Methods, 2014, 11(4): 407–409

    Zhu W, Zhang H. Why do we test multiple traits in genetic association studies? Journal of Korean Statistics Society, 2009, 38(1): 25–27

    Joint Genome-Wide Association Study of Body Mass and Morphological Traits in Turbot ()

    GAO Jin1,2, YANG Runqing1,2①

    (1. Wuxi Fisheries College, Nanjing Agricultural University, Wuxi 214081; 2. Research Centre for Aquatic Biotechnology, Chinese Academy of Fishery Sciences, Beijing 100141)

    To reveal the molecular genetic mechanisms of body mass and morphological traits in turbot () and scan molecular markers and candidate genes, which can be used to improve the target traits, a genome-wide association study (GWAS) was carried out using specific-locus amplified fragment technology (restriction site-associated DNA, 2b-RAD). First, body mass (BM), body length (BL), body width (BW), and caudal peduncle width (CPW) of 441 individuals were measured at about 473 days of growth period in a turbot breeding population. Second, all individuals were genotyped using 2b-RAD, and 23,988 SNPs were obtained after strict quality control. Using a multivariate linear mixed model (mvLMM) for GWAS of traits of BM-BL and BL-BW-CPW, 9 and 2 pleiotropic QTNs were detected for each phenotypic combination, respectively. However, a single-trait linear mixed model (LMM) based on the FaST-LMM algorithm was used for the association analysis of each trait, and the results showed that 4 QTNs were detected in the BM trait, 1 QTN was associated with BL and BW traits, respectively, and no significant locus was found in the CPW trait. A comparison between results of mvLMM and LMM found that mvLMM could detect more QTNs than LMM in GWAS, and the pleiotropic QTNs detected by mvLMM were more biologically meaningful. This study applied mvLMM and LMM to the joint GWAS of body mass and morphological traits in turbot, 17 significant QTNs were detected both using mvLMM and LMM, and 4 of them were detected repeatedly. Furthermore, 12 candidate genes were found by searching the nearest gene of each detected QTN on the whole turbot genome. All of them might be important candidate markers and functional genes, which could influence turbot body mass and morphology. Our study also provided the theory and a reference for marker-assisted selection of body mass and morphological traits in turbot.

    Turbot; Body mass and morphological traits; Genome-wide Association Analysis; Multivariate linear mixed models

    YANG Runqing, E-mail: runqingyang@cafs.ac.cn

    S917.4

    A

    2095-9869(2021)02-0063-08

    10.19663/j.issn2095-9869.20200202001

    http://www.yykxjz.cn/

    高進(jìn), 楊潤(rùn)清. 大菱鲆體重和體尺性狀聯(lián)合GWAS分析. 漁業(yè)科學(xué)進(jìn)展, 2021, 42(2): 63–70

    Gao J, Yang RQ. Joint genome-wide association study of body mass and morphological traits in turbot (). Progress in Fishery Sciences, 2021, 42(2): 63–70

    *中國(guó)水產(chǎn)科學(xué)研究院基本科研業(yè)務(wù)費(fèi)項(xiàng)目“鲆鰈魚分子標(biāo)記輔助配套選育技術(shù)”(2017A001)資助 [This work was supported by the Special Scientific Research Funds for Central Non-Profit Institutes, Chinese Academy of Fishery Sciences (2017A001)].高 進(jìn),E-mail: gaojin427@126.com

    楊潤(rùn)清, 研究員, E-mail: runqingyang@cafs.ac.cn

    2020-02-02,

    2020-02-19

    (編輯 馮小花)

    猜你喜歡
    大菱鲆體尺表型
    家畜體尺自動(dòng)測(cè)量技術(shù)研究進(jìn)展
    基于Kinect相機(jī)的豬彎曲體尺測(cè)量算法研究
    黃海水產(chǎn)研究所“一種大菱鲆油乳化疫苗及其應(yīng)用”獲國(guó)家發(fā)明專利授權(quán)
    肉羊體尺測(cè)量 用上“智慧眼”
    遼寧大菱鲆養(yǎng)殖產(chǎn)業(yè)發(fā)展形勢(shì)分析
    建蘭、寒蘭花表型分析
    體外培養(yǎng)法探討不同蛋白源對(duì)大菱鲆腸道菌群的影響
    GABABR2基因遺傳變異與肥胖及代謝相關(guān)表型的關(guān)系
    慢性乙型肝炎患者HBV基因表型與血清學(xué)測(cè)定的臨床意義
    72例老年急性白血病免疫表型分析
    精品久久久久久久久亚洲 | 国产探花极品一区二区| 亚洲av成人不卡在线观看播放网| 国产三级黄色录像| 在线免费观看不下载黄p国产 | 麻豆久久精品国产亚洲av| 国产爱豆传媒在线观看| 亚洲,欧美精品.| 日韩av在线大香蕉| 校园春色视频在线观看| 特级一级黄色大片| 国产男靠女视频免费网站| 国产三级黄色录像| 99热6这里只有精品| 免费在线观看日本一区| 免费在线观看影片大全网站| 长腿黑丝高跟| 日韩亚洲欧美综合| 啦啦啦观看免费观看视频高清| 国产乱人伦免费视频| 美女免费视频网站| 极品教师在线免费播放| 国产单亲对白刺激| 亚洲午夜理论影院| 亚洲av成人不卡在线观看播放网| 舔av片在线| 美女 人体艺术 gogo| 夜夜躁狠狠躁天天躁| 久久精品91蜜桃| 宅男免费午夜| 国产精华一区二区三区| 久久人妻av系列| 久久香蕉精品热| 婷婷六月久久综合丁香| 欧美xxxx黑人xx丫x性爽| 欧洲精品卡2卡3卡4卡5卡区| 在线观看66精品国产| 搡女人真爽免费视频火全软件 | 亚洲美女视频黄频| 嫩草影院入口| 别揉我奶头 嗯啊视频| 国产白丝娇喘喷水9色精品| eeuss影院久久| 国产高清有码在线观看视频| 欧美高清性xxxxhd video| 在线十欧美十亚洲十日本专区| 亚洲成人精品中文字幕电影| 91久久精品国产一区二区成人| 亚洲欧美精品综合久久99| 久久精品国产亚洲av涩爱 | 久久国产精品人妻蜜桃| 日韩欧美 国产精品| 免费高清视频大片| 日韩有码中文字幕| 在线播放无遮挡| 精品人妻1区二区| 国产精品自产拍在线观看55亚洲| 亚洲一区高清亚洲精品| 白带黄色成豆腐渣| 黄片小视频在线播放| 桃红色精品国产亚洲av| 精品国内亚洲2022精品成人| 日韩精品青青久久久久久| 90打野战视频偷拍视频| 免费看日本二区| 天天一区二区日本电影三级| 久久婷婷人人爽人人干人人爱| 日本 欧美在线| 直男gayav资源| 在线a可以看的网站| 久久婷婷人人爽人人干人人爱| 色精品久久人妻99蜜桃| 国产在线精品亚洲第一网站| 在线播放无遮挡| 91狼人影院| 午夜免费激情av| 午夜精品在线福利| 好男人在线观看高清免费视频| 三级毛片av免费| 国产精品永久免费网站| 精品国内亚洲2022精品成人| 亚洲色图av天堂| 两个人的视频大全免费| 亚洲自拍偷在线| 男女床上黄色一级片免费看| 听说在线观看完整版免费高清| 国产精品久久视频播放| 国产中年淑女户外野战色| 亚洲国产精品久久男人天堂| 91在线观看av| 成年女人看的毛片在线观看| 天堂动漫精品| 无人区码免费观看不卡| 日韩精品中文字幕看吧| 国产精品免费一区二区三区在线| 一级作爱视频免费观看| 搞女人的毛片| 我的女老师完整版在线观看| 久久久国产成人免费| 国产激情偷乱视频一区二区| 少妇人妻一区二区三区视频| 中亚洲国语对白在线视频| 老女人水多毛片| 精品一区二区三区视频在线观看免费| 九九热线精品视视频播放| 97超级碰碰碰精品色视频在线观看| 欧美日韩国产亚洲二区| 亚洲综合色惰| 久久午夜亚洲精品久久| 亚洲无线观看免费| 九九久久精品国产亚洲av麻豆| 午夜福利视频1000在线观看| 中文资源天堂在线| 色噜噜av男人的天堂激情| 国模一区二区三区四区视频| 欧美最黄视频在线播放免费| 亚洲欧美日韩无卡精品| 成人无遮挡网站| 精品一区二区三区av网在线观看| 人人妻人人澡欧美一区二区| 欧美xxxx黑人xx丫x性爽| 成人美女网站在线观看视频| av天堂在线播放| 国产伦人伦偷精品视频| 精品久久久久久久久久久久久| 午夜福利在线观看免费完整高清在 | 欧美日韩国产亚洲二区| www日本黄色视频网| 国产精品久久视频播放| 赤兔流量卡办理| 久久国产乱子伦精品免费另类| 欧美zozozo另类| 日韩欧美国产在线观看| 国内少妇人妻偷人精品xxx网站| 色综合站精品国产| 日韩欧美在线二视频| 亚洲精品在线美女| 我要搜黄色片| 尤物成人国产欧美一区二区三区| 欧美+亚洲+日韩+国产| 国内毛片毛片毛片毛片毛片| 国产色爽女视频免费观看| 国产精品综合久久久久久久免费| 观看美女的网站| 久久伊人香网站| 国产一区二区三区视频了| 国产一区二区亚洲精品在线观看| 色av中文字幕| 成年女人看的毛片在线观看| 免费观看人在逋| 国产久久久一区二区三区| 国内精品美女久久久久久| 亚洲最大成人av| 人妻久久中文字幕网| 久久精品91蜜桃| 综合色av麻豆| 午夜久久久久精精品| ponron亚洲| 少妇裸体淫交视频免费看高清| 直男gayav资源| 99久久久亚洲精品蜜臀av| 在线观看美女被高潮喷水网站 | 午夜影院日韩av| 久久亚洲真实| 9191精品国产免费久久| 搡女人真爽免费视频火全软件 | a级毛片免费高清观看在线播放| 成年女人毛片免费观看观看9| 变态另类成人亚洲欧美熟女| av在线天堂中文字幕| 日本三级黄在线观看| 成人av一区二区三区在线看| 搡老妇女老女人老熟妇| 美女被艹到高潮喷水动态| 亚洲精品日韩av片在线观看| 在现免费观看毛片| 日本熟妇午夜| 亚洲精品粉嫩美女一区| 国产v大片淫在线免费观看| 中文亚洲av片在线观看爽| a级毛片免费高清观看在线播放| 99久国产av精品| 亚洲国产高清在线一区二区三| 欧美日韩中文字幕国产精品一区二区三区| av黄色大香蕉| 国产真实伦视频高清在线观看 | 亚洲在线自拍视频| 亚洲精品成人久久久久久| 久久午夜亚洲精品久久| 亚洲国产精品合色在线| 国产一区二区亚洲精品在线观看| 国产av一区在线观看免费| 成人国产一区最新在线观看| 嫁个100分男人电影在线观看| 欧美高清成人免费视频www| 国内揄拍国产精品人妻在线| 午夜精品久久久久久毛片777| 黄色配什么色好看| 久久久久久国产a免费观看| 亚洲avbb在线观看| 又粗又爽又猛毛片免费看| 毛片一级片免费看久久久久 | 国产精品乱码一区二三区的特点| 变态另类丝袜制服| 又爽又黄a免费视频| 大型黄色视频在线免费观看| 免费人成视频x8x8入口观看| 国产伦精品一区二区三区四那| 一a级毛片在线观看| 搡老妇女老女人老熟妇| 91麻豆av在线| 日韩国内少妇激情av| 亚洲色图av天堂| 国产乱人视频| 国产精品一区二区三区四区久久| 男女视频在线观看网站免费| 午夜免费男女啪啪视频观看 | 精品人妻视频免费看| 午夜免费激情av| 成人永久免费在线观看视频| 亚洲内射少妇av| 成熟少妇高潮喷水视频| 嫩草影视91久久| 97热精品久久久久久| 欧美高清成人免费视频www| 丰满的人妻完整版| 很黄的视频免费| 少妇高潮的动态图| 偷拍熟女少妇极品色| 国产精品一区二区免费欧美| 男女下面进入的视频免费午夜| 18禁黄网站禁片午夜丰满| 国语自产精品视频在线第100页| 久久精品综合一区二区三区| 国产成人a区在线观看| 搡老岳熟女国产| 亚洲人成电影免费在线| 亚洲成人免费电影在线观看| 久久久久久久精品吃奶| 亚洲无线观看免费| 色综合站精品国产| 免费在线观看日本一区| 精品乱码久久久久久99久播| 欧洲精品卡2卡3卡4卡5卡区| 一级a爱片免费观看的视频| 久9热在线精品视频| 白带黄色成豆腐渣| 日本a在线网址| 丁香六月欧美| 天堂动漫精品| 中文在线观看免费www的网站| 色在线成人网| 美女 人体艺术 gogo| 亚洲中文字幕日韩| 最近最新免费中文字幕在线| 亚洲最大成人中文| 欧美性猛交╳xxx乱大交人| 夜夜躁狠狠躁天天躁| 亚洲片人在线观看| 性欧美人与动物交配| 免费在线观看影片大全网站| 免费观看的影片在线观看| 日韩国内少妇激情av| 久久久久久久午夜电影| 老女人水多毛片| 3wmmmm亚洲av在线观看| 欧美色欧美亚洲另类二区| 午夜福利在线观看免费完整高清在 | 窝窝影院91人妻| 黄色女人牲交| 99热这里只有精品一区| 搡老岳熟女国产| 午夜激情福利司机影院| 国产爱豆传媒在线观看| 最近在线观看免费完整版| 1000部很黄的大片| 狠狠狠狠99中文字幕| 中文字幕高清在线视频| 变态另类丝袜制服| 亚洲三级黄色毛片| 美女高潮喷水抽搐中文字幕| 中文字幕人妻熟人妻熟丝袜美| 成人三级黄色视频| 成年免费大片在线观看| 国产乱人视频| 最近在线观看免费完整版| 亚洲最大成人中文| 搡老妇女老女人老熟妇| 欧美+日韩+精品| 欧美一区二区国产精品久久精品| 啦啦啦韩国在线观看视频| 免费av毛片视频| 国产精品久久久久久久久免 | 国产淫片久久久久久久久 | 欧美成人性av电影在线观看| 欧美日本视频| 国产精品爽爽va在线观看网站| 欧美一级a爱片免费观看看| 国产免费av片在线观看野外av| 久久精品国产99精品国产亚洲性色| 亚洲专区中文字幕在线| 国产成人aa在线观看| 亚洲电影在线观看av| 亚洲精品影视一区二区三区av| 国产免费男女视频| 一个人免费在线观看电影| 赤兔流量卡办理| 99在线人妻在线中文字幕| 国产 一区 欧美 日韩| 久久亚洲精品不卡| 精品久久久久久成人av| 一区福利在线观看| 精品99又大又爽又粗少妇毛片 | 最近中文字幕高清免费大全6 | 精品不卡国产一区二区三区| 久久午夜福利片| 成人高潮视频无遮挡免费网站| 俄罗斯特黄特色一大片| 国产伦一二天堂av在线观看| a级毛片a级免费在线| 欧美最黄视频在线播放免费| 99在线人妻在线中文字幕| 日本一本二区三区精品| 国产成年人精品一区二区| 精品人妻一区二区三区麻豆 | 九九在线视频观看精品| 国产成人影院久久av| 一个人看的www免费观看视频| 精品免费久久久久久久清纯| 少妇熟女aⅴ在线视频| av在线天堂中文字幕| 午夜免费男女啪啪视频观看 | 美女高潮喷水抽搐中文字幕| 淫妇啪啪啪对白视频| 黄片小视频在线播放| 国产av一区在线观看免费| 久久亚洲精品不卡| www.熟女人妻精品国产| 香蕉av资源在线| 亚洲美女搞黄在线观看 | 午夜日韩欧美国产| 欧美潮喷喷水| 日本黄色视频三级网站网址| 国产黄色小视频在线观看| 国内久久婷婷六月综合欲色啪| 久久久久亚洲av毛片大全| 国产精华一区二区三区| 欧美不卡视频在线免费观看| 国产精品日韩av在线免费观看| 亚洲精品乱码久久久v下载方式| 91久久精品电影网| 天天躁日日操中文字幕| 国产午夜精品论理片| 国产精品免费一区二区三区在线| 桃色一区二区三区在线观看| 身体一侧抽搐| 嫩草影院精品99| 国产久久久一区二区三区| 国产一级毛片七仙女欲春2| 国产一区二区亚洲精品在线观看| 日本三级黄在线观看| 中文字幕av成人在线电影| 天天躁日日操中文字幕| 一级黄色大片毛片| 久久久久九九精品影院| 日本黄大片高清| 十八禁人妻一区二区| 欧美日韩综合久久久久久 | 免费搜索国产男女视频| 在线a可以看的网站| 欧美成人一区二区免费高清观看| 中文资源天堂在线| 97超视频在线观看视频| 日本与韩国留学比较| 淫秽高清视频在线观看| 蜜桃久久精品国产亚洲av| 亚洲专区中文字幕在线| 午夜激情欧美在线| 亚洲自拍偷在线| 国产精品伦人一区二区| 久久精品久久久久久噜噜老黄 | 午夜免费男女啪啪视频观看 | 哪里可以看免费的av片| 怎么达到女性高潮| 午夜福利免费观看在线| 成人永久免费在线观看视频| 91久久精品国产一区二区成人| 老熟妇乱子伦视频在线观看| 18禁裸乳无遮挡免费网站照片| 在线观看免费视频日本深夜| 国产欧美日韩一区二区精品| 久久欧美精品欧美久久欧美| 久久草成人影院| 午夜福利欧美成人| 国产av一区在线观看免费| 国产私拍福利视频在线观看| 欧美最黄视频在线播放免费| 乱人视频在线观看| 欧美日韩国产亚洲二区| 51午夜福利影视在线观看| 黄色一级大片看看| 午夜老司机福利剧场| 久久香蕉精品热| 欧美成狂野欧美在线观看| 三级国产精品欧美在线观看| 国产中年淑女户外野战色| 自拍偷自拍亚洲精品老妇| 免费看日本二区| 最新中文字幕久久久久| 最近最新中文字幕大全电影3| 精品乱码久久久久久99久播| 99久久精品一区二区三区| 国产精品爽爽va在线观看网站| 好看av亚洲va欧美ⅴa在| 天美传媒精品一区二区| 黄色配什么色好看| 嫁个100分男人电影在线观看| 亚洲一区二区三区色噜噜| 亚洲成人久久爱视频| 国产av在哪里看| 极品教师在线视频| 高清日韩中文字幕在线| 99精品在免费线老司机午夜| 俄罗斯特黄特色一大片| 国产白丝娇喘喷水9色精品| 99riav亚洲国产免费| 身体一侧抽搐| 国产视频内射| 18禁黄网站禁片免费观看直播| 亚洲国产精品久久男人天堂| 亚洲av成人av| 亚洲成av人片在线播放无| 69人妻影院| 变态另类成人亚洲欧美熟女| 午夜精品久久久久久毛片777| 国产精品影院久久| 一卡2卡三卡四卡精品乱码亚洲| 国产伦精品一区二区三区视频9| 又黄又爽又刺激的免费视频.| 国产一区二区在线av高清观看| 99久久99久久久精品蜜桃| 嫩草影院新地址| 网址你懂的国产日韩在线| 日本一二三区视频观看| 91午夜精品亚洲一区二区三区 | 人妻制服诱惑在线中文字幕| 一边摸一边抽搐一进一小说| 国产高清视频在线播放一区| 久久草成人影院| 国产精品久久久久久亚洲av鲁大| 女人被狂操c到高潮| 亚洲美女视频黄频| 国产色婷婷99| 国产精品精品国产色婷婷| 欧美激情久久久久久爽电影| 一本久久中文字幕| 国产三级中文精品| 校园春色视频在线观看| 精品午夜福利在线看| 亚洲中文日韩欧美视频| 我要搜黄色片| 在线观看免费视频日本深夜| 18禁黄网站禁片午夜丰满| 精品人妻视频免费看| 高清在线国产一区| 麻豆国产97在线/欧美| 国产成+人综合+亚洲专区| 男女视频在线观看网站免费| 亚洲成人久久性| 99久久精品国产亚洲精品| 在线观看一区二区三区| 怎么达到女性高潮| 观看美女的网站| 欧美日韩综合久久久久久 | 久久久国产成人精品二区| 九九久久精品国产亚洲av麻豆| 国产探花在线观看一区二区| 国产伦人伦偷精品视频| 成人特级黄色片久久久久久久| 欧美xxxx性猛交bbbb| 18禁黄网站禁片免费观看直播| 一本精品99久久精品77| 九九热线精品视视频播放| 欧美精品啪啪一区二区三区| 男插女下体视频免费在线播放| 亚洲av成人不卡在线观看播放网| 亚洲美女搞黄在线观看 | 91麻豆精品激情在线观看国产| 成人欧美大片| 日日干狠狠操夜夜爽| 亚洲电影在线观看av| 亚洲欧美激情综合另类| 国产高清激情床上av| 国产探花在线观看一区二区| 一个人免费在线观看电影| 99视频精品全部免费 在线| 熟女人妻精品中文字幕| 午夜福利欧美成人| 欧美最黄视频在线播放免费| 69人妻影院| 永久网站在线| 成年女人毛片免费观看观看9| 美女黄网站色视频| 好看av亚洲va欧美ⅴa在| 欧美成人一区二区免费高清观看| 男人舔奶头视频| 欧洲精品卡2卡3卡4卡5卡区| 国产亚洲欧美在线一区二区| 国产白丝娇喘喷水9色精品| 久久久国产成人免费| 91av网一区二区| av在线老鸭窝| 欧美一区二区亚洲| 国产老妇女一区| 99热6这里只有精品| 国产伦精品一区二区三区视频9| 国产av麻豆久久久久久久| 成人午夜高清在线视频| 一区二区三区高清视频在线| 极品教师在线视频| 国产国拍精品亚洲av在线观看| 欧美+日韩+精品| 久久性视频一级片| 亚洲成人精品中文字幕电影| 亚洲人成网站高清观看| 男女之事视频高清在线观看| 欧美日韩黄片免| 免费看日本二区| 精品午夜福利在线看| 一个人看视频在线观看www免费| 高清毛片免费观看视频网站| 亚洲精品一区av在线观看| 精品久久久久久久久久久久久| 亚洲av电影在线进入| 波野结衣二区三区在线| 在线免费观看不下载黄p国产 | 亚洲欧美清纯卡通| 2021天堂中文幕一二区在线观| 国产在视频线在精品| 欧美性猛交╳xxx乱大交人| 中文亚洲av片在线观看爽| 久久99热6这里只有精品| 国产精品一及| 国产精品99久久久久久久久| 亚洲av免费高清在线观看| 亚洲成人中文字幕在线播放| 亚洲av成人精品一区久久| 51午夜福利影视在线观看| av在线蜜桃| 校园春色视频在线观看| 成人特级av手机在线观看| av天堂中文字幕网| 亚洲欧美激情综合另类| 欧美黄色片欧美黄色片| 两性午夜刺激爽爽歪歪视频在线观看| 99在线视频只有这里精品首页| 18禁黄网站禁片免费观看直播| 一卡2卡三卡四卡精品乱码亚洲| 精品一区二区免费观看| 色av中文字幕| av天堂在线播放| 国产淫片久久久久久久久 | 91av网一区二区| 国内精品久久久久久久电影| 久久久久亚洲av毛片大全| 极品教师在线视频| 亚洲熟妇中文字幕五十中出| 黄色一级大片看看| 三级毛片av免费| 国产精品一及| 国产探花极品一区二区| 能在线免费观看的黄片| 婷婷六月久久综合丁香| 搡老熟女国产l中国老女人| 欧美乱妇无乱码| 在线观看一区二区三区| 国产精品人妻久久久久久| 亚洲五月婷婷丁香| 一级a爱片免费观看的视频| 国产精品野战在线观看| 亚洲精品久久国产高清桃花| 无遮挡黄片免费观看| 人妻丰满熟妇av一区二区三区| 欧美在线一区亚洲| 日韩欧美国产一区二区入口| 国产色爽女视频免费观看| av视频在线观看入口| 97碰自拍视频| 真实男女啪啪啪动态图| 我要搜黄色片| 三级男女做爰猛烈吃奶摸视频| 在现免费观看毛片| 成年女人毛片免费观看观看9| 人人妻人人澡欧美一区二区| 亚洲一区高清亚洲精品| 久久精品久久久久久噜噜老黄 | 精品99又大又爽又粗少妇毛片 | 日本黄色片子视频| 亚洲av美国av| 老熟妇乱子伦视频在线观看| 99热这里只有是精品在线观看 | 老司机深夜福利视频在线观看| 能在线免费观看的黄片| 高清日韩中文字幕在线| 欧美日韩瑟瑟在线播放| 少妇人妻精品综合一区二区 | 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 99热这里只有是精品50| 久久6这里有精品| 午夜免费成人在线视频| 看免费av毛片| 亚洲精品久久国产高清桃花| 美女大奶头视频| 色综合婷婷激情| 日韩亚洲欧美综合|