?

大黃魚野生和養殖群體的COI序列變異分析*

2023-10-24 06:30雷鳳玲陳敏芳蒙揚鑫牛素芳吳仁協
廣西科學 2023年4期
關鍵詞:組群粵西大黃魚

雷鳳玲,陳敏芳,蒙揚鑫,牛素芳,吳仁協,**,潘 瀅

(1.廣東海洋大學水產學院,廣東湛江 524088; 2.寧德市富發水產有限公司,大黃魚育種國家重點實驗室,福建寧德 352103)

大黃魚Larimichthyscrocea是石首魚科Sciaenidae中主要的經濟魚類之一,俗稱黃花魚、黃魚、大鮮等,為暖水性集群洄游的中下層魚類,從山東半島以南至雷州半島以東海域均有分布。大黃魚為我國特有的地方性海水魚類,素有“國魚”之稱,大黃魚漁業曾位居中國“四大海洋漁業”之首[1]。早在20世紀50年代,國內學者就開始對大黃魚的漁業生物學、生態學以及種質資源狀況進行了持續調查。以往的研究報道通常將中國沿海的大黃魚分為3個地理族群:岱衢族(黃東海種群),分布在黃海南部至東海中部;閩粵東族(閩粵東種群),從東海南部至珠江口;硇洲族(粵西種群),限于珠江口以西至瓊州海峽以東[2-4]。經過1955-1980年3個時期的過度捕撈,中國近海大黃魚資源從興旺繁盛時代逐漸衰退到嚴重枯竭狀態,早已不再形成漁汛[5]。大黃魚生活史類型從之前的K選擇演化為當前的r選擇,種群結構趨向簡單,資源穩定性變差。

為此,我國政府大力發展東海區的大黃魚人工養殖和增殖放流,以減輕對野生大黃魚的捕撈壓力[3]。然而,養殖實踐和人工放流容易影響魚類種質資源狀況和改變野生種群的遺傳背景?;诙喾N分子標記的大量研究表明,東海區大黃魚養殖群體遺傳多樣性處于較低水平[6-8],野生群體的種質資源整體退化嚴重,可供種質改良所需的遺傳基礎急劇減少[9]。目前有關大黃魚粵西種群遺傳背景的調查研究較少,僅見Wang等[10]、Han等[11]、Liu等[12]應用線粒體D-loop序列或SNP標記分析了粵西大黃魚的群體遺傳變異,且該種群與閩粵東種群的遺傳學關系仍存有疑惑和爭議。此外,東海區大黃魚養殖個體是否對粵西種群的遺傳多樣性產生影響尚未知曉。因此,有必要使用其他分子標記來進一步開展東海區和粵西海域大黃魚不同群體間的遺傳學研究,以了解養殖和野生群體間的遺傳特征差異以及養殖對野生群體的遺傳學影響。

線粒體COI作為線粒體DNA中功能穩定而又具有較多變異的基因序列,已成為動物物種DNA條形碼與系統進化、譜系地理學和遺傳分化等方面研究的重要分子標記。此外,也有研究指出COI基因序列在魚類群體遺傳學研究中具有較好的解析能力,已廣泛應用于鯔Mugilcephalus[13]、棘頭梅童魚Collichthyslucidus[14]、彈涂魚科Periophthalmidae[15]等魚類的群體遺傳多樣性研究。本研究通過分析大黃魚粵西野生群體、閩粵東野生群體、東海區養殖群體的線粒體COI基因序列變異,研究不同群體間的遺傳多樣性、遺傳分化和歷史動態,以期為粵西海域和東海區大黃魚種質資源保護提供基礎資料和科學依據,亦可為養殖群體種質資源改良提供參考。

1 材料與方法

1.1 樣品采集

本研究所用的兩個大黃魚野生群體分別于2017年12月和2019年9月采自粵西的徐聞東岸(XWD,34尾)和硇洲島近岸(NZD,23尾),兩個養殖群體于2011年11月分別采自浙江象山岱衢族大黃魚養殖實驗基地(XSYZ,21尾)和福建寧德三都灣大黃魚養殖漁排(NDYZ,20尾)。實驗魚經形態學鑒定后,取背部肌肉裝入95%酒精溶液中,于-20 ℃保存備用。

1.2 序列擴增和測序

采用傳統的苯酚-氯仿法提取樣品的基因組DNA。利用魚類COI條碼通用引物(FishF1、FishR1)來擴增大黃魚線粒體COI基因序列。PCR的反應體系和擴增程序參照Wu等[16]報道的實驗方法。PCR擴增產物經1.2% 瓊脂糖凝膠電泳檢測合格后送至蘇州金唯智生物科技有限公司測序。

1.3 數據分析

除上述4個群體測定的98條COI序列外,本研究還從GenBank上下載了22條閩粵東野生群體(MYD)的大黃魚COI條碼序列,其中9條序列來自福建中南部沿海(MG574434-MG574437、KX777989、KX777990)和臺灣西岸(KX777985、KX777986、KX777988),13條序列來自珠江口近海(EU595167-EU595177、FJ237998、HQ564537、HQ564538)。本研究分析了包括上述5個群體共120條大黃魚的COI條碼序列。

群體COI條碼序列的堿基組成、序列變異、單倍型確定、單倍型網絡關系以及群體間的單倍型分布頻率Exact檢驗和Fst值等分析參照吳仁協等[17]報道的方法。群體基因交流Nm由DnaSP 5.10軟件基于Gammast值進行估算。在Arlequin Ver 3.5.1.2軟件中,通過3種設定方式來進行群體分子變異分析(AMOVA):一是將5個群體歸為一個組群;二是將5個群體劃分為野生組群(MYD、XWD、NZD)和養殖組群(XSYZ、NDYZ);三是根據群體遺傳分化Fst值將5個群體劃分為3個組群(MYD、XWD、NZD,XSYZ,NDYZ)。上述分析所涉及的核苷酸進化模型采用由jModelTest V2.1.10軟件運算出的DNA序列最佳替換模型(K-2P+G,G=0.876),其余參數為分析軟件設置的默認參數。根據樣品來源情況,將3個野生群體合并為一個野生組群(YSQ)、2個養殖群體合并為一個養殖組群(YZQ)進行歷史動態分析,其中,中性檢驗、核苷酸不配對分布、擴張時間分析等參見吳仁協等[17]報道的方法。由于線粒體COI基因分歧速率與Cytb基因相近。因此,本研究采用其他硬骨魚類Cytb基因2%每百萬年的突變速率作為大黃魚COI條碼序列的突變速率,以推算群體擴張時間。

2 結果與分析

2.1 遺傳多樣性

所分析的120條大黃魚COI條碼序列經比對后,保留共有序列長度有591 bp,共編碼197個氨基酸,序列無插入和缺失或終止密碼子。序列T、C、A、G平均堿基含量分別為27.3%、30.0%、22.9%、19.8%,A+T含量(50.2%)稍微高于C+G含量(49.8%)。其C+G含量在3個密碼子中的分布差異較大,第1密碼子的C+G含量(56.6%)明顯高于第2、3密碼子(分別為43.4%、49.6%)。在591 bp長的序列中共檢測到35個變異位點,其中簡約信息位點有13個,單堿基突變位點有22個。這些變異位點共定義了35處核苷酸替換,包括30處替換和5處顛換。

表1顯示,5個群體的單倍型多樣性以NZD群體最高(h=0.877 5)、NDYZ群體最低(h=0.194 7),核苷酸多樣性以XWD群體最高(π=0.003 1)、XSYZ群體最低(π=0.000 8),表現出明顯的群體間差異?;浳?個野生群體(XWD、NZD)的平均單倍型多樣性(h=0.841 6)稍微高于閩粵東野生群體(MYD,h=0.757 6),但二者的核苷酸多樣性相等(π值均為0.002 8),表明二者的遺傳多樣性水平大致相當。野生組群(YSQ)的單倍型多樣性和核苷酸多樣性(h=0.810 8、π=0.002 9)均遠高于養殖組群(YZQ,h=0.376 8、π=0.000 9),前者是后者的2-3倍。所有個體的單倍型多樣性、核苷酸多樣性和平均核苷酸差異數分別為0.690 6、0.002 2、1.327 3。

表1 大黃魚5個群體的遺傳多樣性指數

2.2 單倍型分布

所分析的120條COI條碼序列共定義34個單倍型,單倍型間的核苷酸差異數為1-9。表2的單倍型分布表明,單倍型H1為5個群體所共享,H8、H15為3個群體所共享,H4、H7、H9、H10、H17為2個群體所共享,其余26個單倍型均為某一個群體所特有。其中單倍型H1出現頻率最高,占5個群體總個數的55%,說明該單倍型是大黃魚種群在長期進化過程中形成的較為穩定的優勢基因型。2個養殖群體均只有3個單倍型,且以單倍型H1的個體數為主,分別占XSYZ、NDYZ群體個數的67%、90%;而3個野生群體的單倍型數量較多,有11-17個。單倍型中介鄰接網絡關系圖呈現出以一個主體單倍型H1為中心、其他單倍型呈輻射狀分布的星狀結構(圖1),未顯示出與序列來源地相對應的譜系結構。這種結構與群體發生擴張的譜系特征相一致。

The sizes of circles are proportional to haplotype frequency

表2 大黃魚5個群體的單倍型分布情況

2.3 遺傳分化和歷史動態

遺傳分化分析結果顯示(表3),XSYZ與NDYZ、MYD群體間的Fst值具有顯著性差異(Fst=0.128 5、0.086 0,P<0.05),表明象山養殖群體與寧德養殖群體、閩粵東野生群體均存在顯著的遺傳分化。3個野生群體間的Fst值均為負值,說明閩粵東野生群體與粵西野生群體內個體間的遺傳差異大于群體間的遺傳差異水平。3個野生群體與NDYZ間的Fst值不具有顯著性差異(Fst=-0.010 6-0.000 2,P>0.05),基因流分析也顯示這4個群體間的Nm值較大(19.55-38.04)(表3),表明3個野生群體與寧德養殖群體間存在遺傳同質性,群體間的基因交流明顯。

表3 大黃魚5個群體間的遺傳分化指數Fst (對角線以下)及基因交流Nm值(對角線以上)

3種AMOVA結果均表明大黃魚5個群體的遺傳變異絕大部分來自群體內部(占總變異的98.03%-99.05%),來自組群間或群體間的遺傳變異僅占極少量(0.04%-1.25%),且統計檢驗差異均不顯著(P=0.153-0.723 >0.05)(表4)。同時,AMOVA結果不支持5個群體劃分為野生組群和養殖組群(Fct=0.000 4,P=0.496),也不支持根據Fst值分為3個組群的假設(Fct=0.012 5,P=0.382)。但是,Exact檢驗分析表明,單倍型在群體間的分布頻率差異具有顯著性(P<0.001),不支持5個群體是一個隨機交配種群的零假設。

表4 大黃魚5個群體的分子變異分析(AMOVA)

表5顯示野生組群(YSQ)的Tajima′sD檢驗(D=-2.274,P=0.002)和Fu′sFS檢驗(FS=-27.438,P=<0.001)均為顯著的負值,提示野生大黃魚可能經歷了群體擴張事件。野生組群的核苷酸不配對分布呈單峰分布(圖2),其擬合優度檢驗的SSD值(0.000 3)和Hri指數(0.022)均統計不顯著(P=0.94、0.83)(表5),表明所觀測的野生組群核苷酸不配對分布與群體擴張模型相吻合。根據群體數量擴張參數τ值(1.30),推算出野生群體約在10.9萬年前即晚更新世時期經歷了群體數量擴張事件。群體擴張后、前的θ1/θ0比值為19.4,揭示了野生群體發生擴張后其有效群體數量有實質性的增大。在養殖組群(YZQ)中,Tajima′sD檢驗(D=-1.445,P=0.062)和Fu′sFS檢驗(FS=-1.506,P=0.123)均為統計檢驗不顯著的負值(表5),該結果與群體發生擴張的分子信號不相符。雖然養殖組群的核苷酸不配對分布的擬合優度檢驗結果(SSD=0.009 3、P=0.50,Hri=0.196、P=0.71)并未顯著偏離群體擴張模型(表5、圖2),但其群體擴張參數θ1/θ0比值僅為1.2,提示養殖群體的有效群體數量在其歷史動態過程中變化并不明顯??梢?歷史動態分析結果支持野生大黃魚經歷了晚更新世的群體數量擴張,但不支持養殖大黃魚發生群體擴張事件。

The histogram is the observation distribution,and the curve is the expected distribution under the population expansion model.

表5 大黃魚野生組群和養殖組群的擴張參數以及擬合優度檢驗

3 討論

3.1 大黃魚野生和養殖群體的遺傳多樣性比較

本研究的大黃魚野生群體和養殖群體顯示出不同的遺傳多樣性模式,前者為高h(>0.5)、低π(<0.005)類型,后者為低h、低π類型。這種差異在以往的大黃魚群體線粒體COI序列變異研究中也有相同的報道[18,19]。歷史動態分析結果支持野生大黃魚經歷了晚更新世的群體擴張事件,但不支持養殖大黃魚發生群體擴張事件。因此,野生大黃魚的這種遺傳多樣性模式很可能是因其在種群數量動態演化史中發生了群體擴張,即經歷了一個由較小的有效群體快速擴張成一個較大群體的過程,該過程產生的許多新突變可增加單倍型多樣性,但因缺乏足夠的進化時間來積累核苷酸序列變異,最終形成高h、低π模式[17]。養殖大黃魚呈現低h、低π模式可能是因為其群體經歷過近期瓶頸效應或由單一、少數群體所產生的奠基者效應引起的[19]。與野生群體相比,本研究的養殖大黃魚群體遺傳多樣性下降較為明顯(前者是后者的2-3倍)。這種下降與基于其他分子標記研究東海區大黃魚養殖群體遺傳變異的分析結果相一致[6,9,20-22]。以往研究普遍認為用于養殖繁育的親本數量少、遺傳背景單一,由此引起的遺傳瓶頸效應及其伴隨發生的遺傳漂變和近交衰退等作用,是導致東海區養殖大黃魚種質資源退化的主要原因[7,20-22]。此外,養殖過程中的人為選擇也可能是造成養殖大黃魚群體遺傳多樣性下降的另一個因素。人工選育通常篩選目標群體的優良性狀,隨著世代培育及代際遺傳,其定向選擇勢必增強,容易引起群體中某些稀有基因喪失和雜合子缺失,從而造成子代群體遺傳多樣性的丟失[23]??梢?不同進化力量在大黃魚野生和養殖群體中的作用機制導致了二者遺傳多樣性特征的顯著差異。

與其他研究報道的大黃魚線粒體COI序列遺傳多樣性相比,如呂泗野生群體(h=0.915、π=0.004 12)和浙閩2個養殖群體(h=0.503、π=0.001 53)[18],以及黃東海8個野生群體(h=0.842 5、π=0.003 454)、福建2個養殖群體(h=0.463 5、π=0.001 1)、江浙4個養殖群體(h、π均為0)[19],本研究的野生和養殖群體的遺傳多樣性均低于上述研究結果(因共同親本繁育而來的江浙4個養殖群體除外),提示所檢測的大黃魚種質資源狀況不容樂觀。然而,Han等[11]報道了黃東海的大黃魚野生群體線粒體D-loop序列具有較高的遺傳多樣性(h=0.989 5-1.000 0、π=0.021 2-0.024 6)。除了分子標記、樣品數量、采樣時間的差異引起不同研究的差別之外,一個不爭的事實是過度捕撈導致中國近海的大黃魚野生群體資源早已處于嚴重枯竭狀態。而養殖群體通常來源于少數野生大黃魚繁育的后代個體,有效親本數量較小導致種群出現遺傳瓶頸效應和近交衰退[20,21]。因此,本研究未能在野生群體和養殖群體中檢測到較為豐富的遺傳多樣性屬于情理之中。

與東海區大黃魚不同的是,粵西大黃魚尚無長期的人工養殖或種苗放流增殖歷史。因此,影響遺傳多樣性的瓶頸效應、遺傳漂變和近親雜交等不利因素在大黃魚粵西種群遺傳中的作用應當不明顯,粵西種群可能保存相對較多的稀有等位基因。然而,本研究的粵西野生群體遺傳多樣性水平與閩粵東野生群體相當,明顯低于上述報道的黃海野生群體的遺傳多樣性(h=0.915-0.946、π=0.004 12-0.004 70)[18,19]。究其原因可能有2個:一是粵西種群是大黃魚3個地理種群中產卵群體數量、資源量及捕撈產量、棲息地面積等均最少的一個種群[2,4],其有效群體數量相對較少,所蘊含的遺傳多樣性豐富度不如黃東海種群;二是浙閩不僅是野生大黃魚的主要產區,也是主要的養殖區和人工放流區,其養殖逃逸及放流苗種數量較多[19,22],東海區的這些養殖個體可能通過中國沿岸流進入南海,與南海的野生大黃魚發生基因交流,從而影響了粵西大黃魚野生群體的遺傳多樣性。研究表明養殖魚類可對野生群體產生遺傳稀釋或基因漸滲,導致野生群體遺傳多樣性的降低[23]。

與彈涂魚科(h= 0.913-0.987、π=0.002-0.005)、鱭屬Coilia(h= 0.556-0.933、π=0.002-0.005)、石斑魚屬Epinephelus(h= 0.884-0.955、π=0.002)、銀鯧Pampusargenteus(h= 0.619、π=0.002)等魚類[17]相比,本研究的粵西野生大黃魚群體遺傳多樣性(h=0.757 6-0.805 7、π=0.002 8-0.003 1)處于中等水平。這與Wang等[10]和Han等[11]基于線粒體D-loop序列的研究結果相一致??梢?粵西大黃魚野生群體的遺傳多樣性水平雖比不上南黃海野生群體,但該群體仍保持了一定程度的遺傳變異,具有較強的環境適應能力、生存能力和進化潛力。此外,考慮到粵西大黃魚的群體資源開發和利用程度,以及受人工養殖影響程度均是大黃魚3個地理種群中相對較輕的一個,粵西種群可能是中國近海大黃魚非常寶貴的一個天然種質資源庫,可視為大黃魚增養殖繁育群體的一個潛在的適合補充群體。因此,必須對大黃魚粵西種群進行有效的資源保護和科學管理,以維護其良好的種質資源狀況。

3.2 大黃魚不同來源群體間的遺傳分化

Fst值和Exact檢驗分析表明,大黃魚5個群體間存在一定程度的遺傳差異,主要是象山養殖群體與寧德養殖群體、閩粵東野生群體之間具有顯著的遺傳分化。這種遺傳分化在其他研究報道的大黃魚群體間也有發現。例如,諶微等[19]的線粒體COI序列分析顯示黃東海大黃魚的養殖和野生兩大組群間具有顯著的遺傳分化,但組群內群體間無遺傳差異;Wang等[7]的SSR標記研究也表明大黃魚養殖和野生群體之間存在明顯的遺傳結構。研究表明,養殖過程中發生的奠基者效應、隨機遺傳漂變以及人工選育是造成養殖和野生群體間遺傳分化的主要原因,這在其他經濟魚類的群體遺傳變異研究中已有廣泛的報道,如大西洋鮭Salmosalar[24]、草魚Ctenopharyngodonidella[25]、斜帶石斑魚Epinepheluscoioides[26]等。由于遺傳多樣性較低的養殖魚苗大規模放流到野生種群的生活海區中,一旦產生基因交流就會改變野生種群的遺傳組成和適合度,從而降低野生種群的生存能力和適應環境能力[26],這對大黃魚野生群體資源的保護是非常不利的。因此,對于目前東海區大規模開展的大黃魚苗種增殖放流工程需謹慎實施[27],尤其是在當前養殖群體遺傳多樣性已出現明顯下降之時。此外,本研究發現象山養殖與寧德養殖群體間存在中等偏高程度的遺傳分化水平。這兩個養殖群體均來源于東海區大黃魚的主養區,也都經歷過長期的人工選育,該結果表明人工選育的強度和作用在東海區大黃魚的主要養殖群體中已呈現明顯的遺傳效應。這提示在大黃魚養殖過程中,除了注重親本來源外,對養殖實踐中的人為選擇問題也應當給予足夠的重視。

本研究檢測到寧德養殖群體與3個野生群體之間具有高度的遺傳同質性,即大黃魚閩粵東種群與粵西種群間存在明顯的基因交流。這在以往的大黃魚群體線粒體DNA序列研究中也有發現。Han等[11]研究表明黃海、東海、南海的4個大黃魚野生群體間D-loop序列分歧很小,可視為同一個種群。Wang等[8]聯合COI和Cytb序列分析顯示黃海、東海5個大黃魚野生群體的2個譜系遺傳分化并不顯著。大黃魚在分布區缺乏顯著的地理譜系結構可能反映了冰期后大黃魚在中國近海陸架區發生重新殖化的歷史進程。本研究的歷史動態分析表明大黃魚野生群體經歷了晚更新世的群體擴張,且群體間只檢測到一個單倍型譜系??紤]到目前野生大黃魚遺傳多樣性最豐富的群體位于黃海南部,其分布中心位于東海區??赏茰y中國近海大黃魚種群的冰期避難所應該只位于東海的深水區,在晚更新世冰期之后,來自冰期避難所的大黃魚孑留群體向黃、東、南海陸架區擴張,擴散后新建立的群體尚未在遷移和遺傳漂變之間取得平衡[17,28],這可能是大黃魚閩粵東種群與粵西種群間遺傳差異不顯著的主要歷史原因。這種歷史動態進程在中國近海魚類種群遺傳研究中具有普遍性[17,28,29]。

大黃魚是一種近岸洄游性的海洋魚類,幼體浮游期較長(3-4周)、游泳能力較強[1],具有長距離擴散能力[11]。對該物種而言,東海和南海北部的海洋環境缺乏明顯的隔離因素[15,17]。因此,東海區的養殖逃逸和人工放流大黃魚有可能借助中國沿岸流進入南海北部,從而促進兩海域間大黃魚群體的基因交流。通常群體間Nm>1就能發揮遺傳均質化作用,可有效抑制由遺傳漂變引起的遺傳分化[30]。只要群體的每一個世代有一個個體能有效遷移到另一個群體的繁育群體中,就可阻止群體間因為遺傳漂變或選擇作用而產生的遺傳分化[29]。因此,養殖逃逸、人工放流以及大黃魚較強的擴散能力不但使東海區大黃魚養殖和野生群體間具有很高的遺傳一致性,還可能促進了大黃魚閩粵東種群和粵西種群間發生遺傳均質化。

考慮到本研究所分析的群體數量、樣品涵蓋范圍以及分子標記存在的不足之處,本研究還難以對中國近海大黃魚的種群劃分問題進行深入的討論。有必要在后續研究中收集更多的大黃魚野生和養殖群體樣品,采用分辨力更高的分子標記(如D-loop、SSR、SNP等)[31,32],綜合探討中國近海大黃魚的種群劃分和遺傳資源保護。

4 結論

本研究的線粒體COI序列分析表明,大黃魚野生群體和養殖群體具有完全不同的遺傳多樣性特征,這種差異反映了不同進化力量在大黃魚群體遺傳中的作用機制。晚更新世的群體擴張事件是形成現有野生群體遺傳多樣性的重要因素,而養殖繁育親本數量少和人工選育則可能是導致目前養殖群體遺傳多樣性下降的主要原因。大黃魚粵西野生群體的遺傳多樣性處于中等水平,具有一定程度的遺傳變異,可視為大黃魚增養殖繁育群體的一個潛在的適合補充群體。當前的養殖實踐和人工選育可能導致東海區大黃魚群體間產生顯著的遺傳分化,包括養殖和野生群體間、不同養殖群體間。晚更新世冰期后的大黃魚歷史殖化事件以及東海區大黃魚養殖個體向南海北部的擴散,可能促使大黃魚閩粵東種群(包括野生和養殖)和粵西種群具有遺傳同質性。本研究結果提示,中國近海大黃魚種質資源總體狀況不容樂觀,對于當前東海區大規模開展的大黃魚增殖放流工程需謹慎實施,養殖實踐、增殖放流、人工選育等因素對養殖大黃魚群體遺傳特征產生了重要影響,這應當引起高度重視和需要長期、持續的研究監測,以保護好大黃魚這一我國特有的寶貴魚類資源。

猜你喜歡
組群粵西大黃魚
350人參會!恒興蝦苗在粵西再次火爆,為何深受養戶青睞?
加州鱸在粵西掀起養殖潮!看上上生物如何引領行業轉型升級
73個傳統建筑組群組團出道!帶你活進從前的慢時光
“組群”“妙比”“知人”:小學語文古詩群文閱讀的三個途徑
粵西地區慢性丙型肝炎患者病毒基因型分布特征
28元/斤的輝煌不再!如今大黃魚深陷價格“泥沼”,休漁期或初現曙光
QC新七大工具之五:矩陣圖法
寧德迎來大黃魚豐收季
磁盤組群組及iSCSI Target設置
膳食鋅對飼料引起的大黃魚銅中毒癥的保護作用
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合