?

基于UPLC-Q-Exactive MS 構建麻痹性貝類毒素的特征指紋溯源技術*

2024-02-24 08:45鄭關超王瀟瀟翟毓秀譚志軍吳海燕
海洋與湖沼 2024年1期
關鍵詞:貽貝貝類攝食

張 帆 鄭關超 王瀟瀟 翟毓秀 譚志軍 吳海燕①

(1. 農業部水產品質量安全檢測與評價重點實驗室 中國水產科學研究院黃海水產研究所 山東青島 266071; 2. 中國海洋大學食品科學與工程學院 山東青島 266003; 3. 海水養殖生物育種與可持續產出全國重點實驗室 中國水產科學研究院黃海水產研究所 山東青島 266071)

麻痹性貝類毒素(paralytic shellfish toxins, PSTs)主要由Alexandrium、Gymnodinium和Pyrodinium等海洋甲藻產生, 通過貝類濾食進入食物鏈, 在世界范圍內造成嚴重生態和食品安全風險, 是國際社會重點關注目標和管控對象(Viscianoetal, 2016; Nicolasetal, 2017)。研究表明PSTs 分布具有地域性特點(Villalobosetal, 2019), 而其風險形成過程主要受兩個方面的影響: 首先是產毒藻的種類及來源差異, 鏈狀亞歷山大藻(Alexandriumcatenella)主要分布在渤海、黃海和東海海域(顧海峰等, 2011; 唐瑩瑩等,2018), 而太平洋亞歷山大藻主要分布在東海和南海,鏈狀裸甲藻在福建近海曾多次導致中毒事件(于仁成等, 2020)。其次, 不同地區貝類如紫貽貝攝食產毒藻后毒素成分也有顯著差異, 在2019 年大連大窯灣海域紫貽貝中PSTs 成分以C1&2 和GTX2&3 為主(許道艷等, 2014), 而在浙南海域的紫貽貝中PSTs 成分則以C2、GTX5 和NEO 為主(張樹剛等, 2011), 這種區域差異在其他研究中也得到了充分的證明(Yaoetal,2019; Liangetal, 2022)。因此, 貝類PSTs 組分的區域差異對其風險具有一定的指示作用, 可有效指向風險源地并用于產地溯源, 對于提升貝類中PSTs 風險的監管效果具有重要作用。

近年來, 指紋溯源技術廣泛應用于生物地球化學、醫學領域、環境監測、食品安全等多個領域, 通過對生物體指紋生物標志物的篩選, 結合數據庫比對可快速實現源頭識別(劉靜等, 2022; 冷雪等,2023)。該技術是基于代謝組學技術, 以高通量、高靈敏度、高分辨率的現代儀器分析方法為手段, 對細胞、體液和組織中全部代謝物進行定性與定量分析,隨后結合多元統計分析來構建指紋溯源技術模型(孫曉珊等, 2021; 高淑芳等, 2022)。其中液相色譜-質譜聯用(LC-MS)(徐天潤等, 2020)能夠在任何單一分析平臺中實現最高的代謝組覆蓋率, 并且需要最少的樣品前處理, 包括蛋白質沉淀和代謝物提取(Bujaket al, 2015), 適宜的提取溶劑、提取時間和處理溫度等條件是影響前處理的重要參考指標。有文獻報道, 基于LC-MS 的代謝指紋識別鯉魚和虹鱒魚精漿中的代謝物, 篩選出精子質量的新型生物標志物, 為優化人工繁殖奠定了基礎(Dietrichetal, 2019)。同時也有研究針對澳大利亞翡翠貽貝新西蘭綠唇貽貝和進行代謝指紋檢測, 發現兩種貽貝之間存在明顯的代謝差異, 根據差異代謝物可以實現貽貝產地溯源(Rochfortetal, 2013)。此外, 也有研究采用同位素標記(Zhaoet al, 2019)、脂質(Shinetal, 2008)、微量元素(Morrisonetal, 2019)等技術識別雙殼貝類腸道內容物的方式,實現貝類產地溯源(Gaoetal, 2006)。

采用高通量篩查技術獲得的數據規模大且來源復雜, 因此需要結合化學計量學方法, 降低數據維度篩選出具有識別能力的變量(徐天潤, 2020)。多變量分析是通過降低維度來簡化數據的復雜程度, 并通過相關軟件進行結果可視化, 主要包括主成分分析(principal component analysis, PCA)、正交偏最小二乘判別分析(orthogonal signal correction partial least squares-discriminant analysis, OPLS-DA)、聚類分析等。其中PCA 主要通過質量控制樣本(QC)在PCA 圖中分布情況對數據進行初步考察, 譚芷晴等(2023)為篩選炔諾酮暴露斑馬魚后的生物標志物, 運用PCA分析發現處理組與代謝組有明顯分離。劉天亮等(2022)運用PCA 分析成功將來自山東、河南和河北產區的金銀花樣品分成了3 類。而OPLS-DA 可以使分類信息主要集中在主成分上, 并且通過變量投影重要度(Variable importance in projection, VIP)來確定差異化合物, 判別效果更具有說服力(李思源等, 2021)。然而, 基于多變量分析的代謝指紋從未用于研究攝食麻痹性貝類毒素的紫貽貝。

本研究選擇我國近海主要麻痹性貝類毒素產毒藻種, 通過室內暴露貽貝模擬指紋信息傳遞過程, 采用UPLC-Q-Exactive MS 對五種前處理方法進行比較分析, 結合PCA 和OPLS-DA 化學計量學方法, 挖掘攝食不同產毒藻后貽貝體內的化學差異, 篩選區分攝食不同產毒藻貽貝中的特征指紋物質, 構建特征指紋溯源技術模型, 快速精準的識別潛在風險以便提高麻痹性貝類毒素風險源頭識別。

1 材料與方法

1.1 實驗材料

實驗所用藻株包括鏈狀亞歷山大藻(A.catenella,GY-H25)、微小亞歷山大藻(A.minutum, GY-H46)和鏈狀裸甲藻(G.catenatum, GY-H65)購自上海光語生物科技有限公司。實驗室內以f/2 培養液單種培養, 溫度為(20±1) °C, 光照為6 000 lx; 光暗比12 h∶12 h。餌料藻選用小球藻(Chlorellavulgaris), 置于-20 °C冰箱保存。

樣品前處理及色譜質譜條件優化采用取自秦皇島貽貝重點養殖區的紫貽貝(Mytilusgalloprovincialis),經檢測為麻痹性貝類毒素陽性貝, 毒素含量為1 008 μg STXeq/kg。

指紋溯源技術建立采用購于山東青島碼頭紫貽貝(Mytilusgalloprovincialis), 平均規格為: 殼長(40.0±1.9) mm, 殼寬(22.6±0.8) mm, 總重(6.0±0.6) g,購入后先在實驗室進行馴化, 每天持續通氣并更換海水。實驗共設置3 組, 每200 只紫貽貝為一組, 培養于24 h 通氣的天然海水養殖容器環境下。

1.2 實驗方法

1.2.1 紫貽貝室內暴露實驗 各取1 mL 混勻后的鏈狀亞歷山大藻、微小亞歷山大藻和鏈狀裸甲藻的藻液, 加入20 μL 魯哥試劑(Lugol’s agent)固定, 在顯微鏡下計數藻細胞密度。按照A.catenella組和A.minutum組投喂量約為1×105cells/(inds.·d),G.catenatum組投喂量約為 2×105cells/(inds.·d), 早晚各喂一次,輕微曝氣。實驗共持續14 d, 投喂產毒藻7 d, 投喂餌料藻7d 。分別取0 (對照)、1、3、7 和14 d 的紫貽貝樣品, 每組樣品設置3 個生物學重復, 毒素分析樣品為3 只貝混樣; 特征物質分析樣品為6 只貝混樣, 均解剖紫貽貝并采集全部軟組織, 于-80 °C 下冷凍保存。

1.2.2 產毒藻和貽貝中麻痹性貝類毒素組成分析 將指數生長期的藻液取20 mL 于離心管中, 2 760 ×g離心5 min, 收集藻細胞置于15 mL 離心管中, 加入5 mL 1%乙酸水溶液。將其置于冰浴中, 使用超聲波細胞破碎儀(JY92- ⅡN, 寧波新芝生物科技股份有限公司)對藻類樣品進行破碎(設置全程時間5 min、間隔時間2 s、功率比50%)。鏡檢確認藻細胞完全破碎后, 將離心管置于高速冷凍離心機(Himac CR 22G,日本 Hitachi 公司)中, 2 760 ×g離心5 min, 供測試。將貽貝開殼用生理鹽水清洗后取軟組織, 瀝水均質, 參考Boundy 等(2015)等方法進行麻痹性貝類毒素提取及測試。

1.2.3 特征指紋物質的前處理方法 取(1.00±0.02) g麻痹性貝類毒素陽性紫貽貝樣品轉移至5 mL 離心管(分為5 組, 每組6 個平行), 按照表1 中五種方法進行前處理, 提取液經12 470 ×g離心5 min 后, 供測試。

1.2.4 液相條件 Q-Exactive 液相色譜質譜聯用儀(Q-Exactive, 美國 Thermo 公司)C8 色譜柱體系:Kinetex C8 (150 mm×2.1 mm, 2.6 μm); 梯度洗脫程序:0~0.5 min, 20% A; 0.5~5.0 min, 20%~90% A;5.0~14.0 min, 90% A; 14.0~15.0 min, 80% A。Amide色譜柱體系: TSK-Amide-80 柱(2 mm×15 cm, 3 μm);梯度洗脫程序: 0~1.0 min, 80% A; 1.0~4.0 min, 80%~40% A; 4.0~6.5 min, 40% A; 6.5~8.0 min, 40%~80%A;8.1~10.0 min, 80% A。

柱溫: 40 °C; 進樣量: 10 μL; 流速: 0.35 mL/min;流動相為: A 相為95%乙腈水(含2 mmol/L 甲酸銨,50 mmol/L 甲酸), B 相為水(含2 mmol/L 甲酸銨,50 mmol/L 甲酸)。

1.2.5 質譜條件 Q-Exactive 液相色譜質譜聯用儀使用加熱電噴霧離子源(HESI), 掃描模式: 全掃描/數據依賴二級子離子掃描(Full MS/dd MS2); 噴霧電壓: 3 500 V; 毛細管和加熱器溫度: 320 和50 °C; 鞘氣和輔助氣: 40 和10 arb; 分辨率用半峰高處的全峰寬(FWHM)表示, 一級全掃描分辨率: 70 000 FWHM;trap 最大容量(AGC target): 1×106; C-trap 最大注入時間: 150 ms, 數據依賴二級離子掃描分辨率: 17 500 FWHM; C-trap 最大容量(AGC target): 1×105; C-trap最大注入時間: 50 ms; 質譜掃描范圍: 質荷比(m/z)100~600 和600~1 500。

1.3 數據處理

1.3.1 樣品前處理及色譜質譜分析 按照時間順序在色譜圖中隨機選取峰強度大于1×105的10 個離子峰, 分別計算每種方法6 個樣品的峰強度均值, 計算每種方法貽貝樣品中所測的10 個離子峰峰強度及其變異系數(Coefficient of variation, CV=標準差/均值)。5 種前處理方法所測峰強度以均值表示, 多組間比較采用隨機區組的單因素方差分析, 若組間差異有統計學意義, 進一步的兩兩比較采用Bonferroni 檢驗。以P<0.05 為差異有統計學意義。

1.3.2 特征指紋物質分析 原始數據使用Compound Discoverer 3.1 軟件進行峰檢測、保留時間對齊等預處理。采用數據庫(ChemSpider, mzCloud 等)從精確質量數、同位素組成及分布與預測分子式的吻合度、一級質譜/二級質譜碎片和數據庫匹配等方面進行化合物初步鑒定, 其中設置參數如下: 保留時間最大漂移值0.2 min; 質量偏差<5×10-6; 信號強度偏差<30%;峰強度最小值度50 000, 同時設置基質空白去除, 不進行統計學分析。由于特征物質沒有標準品, 采用STX 標準品曲線進行相對定量分析計算其濃度。

1.4 統計分析

采用Excel 2021、Origin 2023、SIMCA 14.1 和IBM SPSS Statistics 22 軟件對數據進行統計分析和繪圖。在Excel 2021 中將數據整理為二維矩陣形式, 保留P值、m/z值、保留時間、分子式、峰面積等信息,根據P值篩選出差異顯著的化合物(P<0.01); 以化合物的峰面積為變量,m/z值作為觀測ID 將數據導入SIMCA 14.1 進行PCA 與OPLS-DA 分析; 根據VIP>1作為特征指紋物質。采用逐步判別法對3 組樣品的13 個復合指紋物質特征建立3 個群體的Fisher 線性判別函數。

2 結果與分析

2.1 樣品前處理及色譜質譜條件選擇

通過對比不同色譜和質譜條件結果如圖1 所示,兩種色譜柱均為方法D 檢測到的化合物數量最多。采用C8 色譜柱時(圖1a), 正離子模式m/z100~600 掃描范圍下檢測到的化合物數量最多為3 374 個, 負離子模式m/z100~600 掃描范圍次之; 采用Amide 色譜柱時(圖1b)負離子模式m/z100~600 的掃描范圍下檢測到的化合物最多為1 785 個, 正離子模式m/z600~1 500 掃描范圍次之。因此方法D 檢測到的化合物達到了6 811 個, 更加適合對特征物質的提取。

圖1 不同前處理方法在每種色譜質譜條件下的特征物質數量Fig.1 The amount of compounds in each mode for different pretreatment method

此外, 根據不同方法的峰強度變異系數(CV 值)來分析不同的質譜條件。采用C8 色譜柱時, 離子峰強均值在正離子模式m/z100~600 和m/z600~1 500掃描范圍下大于1×107且差異無統計學意義(P>0.05);比較其CV 值, 方法D 在正離子模式m/z600~1 500掃描范圍下CV 值最小。采用Amide 色譜柱時, 方法D 在正離子模式m/z600~1 500 掃描范圍下的峰強最高為5.49×107且差異具有統計學意義(P<0.05), CV 值<25% (表2)。因此, 選擇方法D 在C8 色譜柱正離子模式m/z100~600 的掃描范圍和Amide 色譜柱正離子模式m/z600~1 500 的掃描范圍下更適用于貽貝的特征物質分析, 覆蓋的化合物數量占總檢測量的40.4%。在方法D 兩種色譜質譜條件下, 分別在QC 樣品色譜圖的前、中、后隨機共選取10 個離子峰進行分析, 其CV 值均≤1.00%, 各峰面積和峰強度的CV 值均<25.00%。以上結果表明本研究樣品分析過程中儀器的精密度、穩定性及重復性均較好, 滿足分析要求。

表2 不同前處理方法測得的峰強及CV 值Tab.2 Peak intensity and CV values of different pretreatment methods

2.2 紫貽貝攝食不同產毒藻后體內麻痹性貝類毒素分析

紫貽貝攝食不同產毒藻后麻痹性貝類毒素的變化情況如圖2 所示, 在前7 天為毒素蓄積階段, 紫貽貝中的麻痹性貝類毒素含量不斷上升, 由暴露1 天毒素含量最低為54.1 μg STXeq/kg, 到7 天時毒素含量最高為5 538 μg STXeq/kg; 7 天到14 天為代謝階段(停止投喂產毒藻), 至 14 天時毒素含量最低降至267 μg STXeq/kg。

圖2 紫貽貝攝食不同產毒藻后體內麻痹性貝類毒素含量及組成Fig.2 Variation of toxin content and composition after ingestion of different toxin producing algae by M. galloprovincialis

在投喂產毒藻的紫貽貝體內共檢出11 種麻痹性貝類毒素成分, 不同組樣品毒素成分差異顯著。投喂A.catenella的紫貽貝中麻痹性貝類毒素主要成分為C1&2 和GTX5, 在整個實驗中C1&2 占比從第1 天68.6%下降到第14 天的49.8%; 投喂A.minutum的紫貽貝中麻痹性貝類毒素主要成分為GTX1&4 和GTX2,且在整個實驗期間各成分占比穩定; 投喂G.catenatum的紫貽貝中麻痹性貝類毒素主要成分為dcGTX2&3、dcSTX 和C1&2, 其中dcSTX 占比從第1 天的10.9%增高到第14 天的21.5%; 而C2 的占比從21.5%降至4.02%。因此, 攝食不同產毒藻后紫貽貝體內的毒素成分有明顯差異且隨時間發生代謝轉化, 還需對紫貽貝體內的代謝物進行檢測。

2.3 多元統計分析

2.3.1 主成分分析 本研究使用上述乙腈/甲醇/丙酮(體積比為1∶1∶1)作為提取劑的提取方法結合篩選出的兩個離子通道, 經過軟件處理后, 提取了3 886 個化合物, 利用P<0.01 篩選差異化合物并結合毒素成分進行PCA 分析, 如圖3 所示, 前兩個主成分分別解釋了總變量的 36.9%和 18.3%, 共解釋了55.2%的變量。三組樣品分別沿PC1 或PC2 軸分離,表明攝食不同產毒藻貽貝中化合物存在差異。

圖3 投喂三種不同產毒藻后貽貝中差異化合物PCA 圖Fig.3 PCA of different compounds and toxins in M.galloprovincialis fed with three different toxic microalgae

2.3.2 正交偏最小二乘判別分析 對投喂A.catenella、A.minutum和G.catenatum的貽貝經篩選后結合毒素進行OPLS-DA 分析, 如圖4a~4c 所示。三組樣品中的化合物顯著分離, 自變量擬合指數分別為0.724、0.860、0.846, 自變量擬合指數分別為1.00、1.000、0.969, 模型預測指數(Q2)為0.985、0.982、0.988, 表明當前OPLS-DA 模型穩定可靠, 具有良好的預測性。經OPLS-DA 模型分析分別生成S-plots 圖, 如圖4d~4f 所示。OPLS- DA 模型中變量VIP>1.0 說明該變量對整體模型的貢獻度高于平均水平, 作為本研究的毒素和特征指紋物質。此外, 毒素成分NEO 和dcNEO 作為攝食A.catenella和G.catenatum后紫貽貝的特有成分也被作為指紋物質,共計13 種復合指紋物質。根據精確分子質量、保留時間以及二級質譜碎片等與數據庫進行比對參數如表3 所示。

表3 復合指紋物質Tab.3 Composite fingerprint substances

2.4 判別分析

根據篩選的復合指紋物質, 建立基于Fisher 判別函數的一般判別方法。以篩選出來的13 種復合指紋物質作為判別分析的自變量, 對貽貝樣品進行多變量判別分析, 排除了7 種對分類作用較小的指紋物質,僅保留6 種指紋物質用于判別分析, 建立Fisher 線性判別函數如下:

式中:Y分別表示暴露于不同產毒藻貽貝的判別得分;C329.22、C381.14、C461.31、C374.21、CNEO、CdcNEO分別表示分子質量為329.22、381.14、461.31、374.21 的特征指紋物質和毒素NEO、dcNEO?;卮鷻z驗正確判別率為100%, 交叉驗證正確率為88.9%, 表明6 種特征指紋成分對貽貝的判別效果較好。

3 討論

麻痹性貝類毒素的風險形成過程, 亦是從產毒藻信息向貝類傳遞的過程。通過傳統的監管手段不僅耗費大量的人力物力財力, 對于近年來常見的小規模有害赤潮爆發作用有限(Fernandes-Salvadoretal,2021)。研究利用UPLC-Q-Exactive MS 結合PCA 分析與OPLS-DA 分析, 從攝食3 種麻痹性貝類毒素產毒藻的紫貽貝樣品里共篩選出13 種復合指紋物質,對攝食不同藻的紫貽貝進行成功溯源。由于有毒藻種類和生態學特征不同, 其所產PSTs 的種類和含量也有很大差異(Guetal, 2013), 因此造成了貝類中PSTs風險表征差異明顯, 這為識別PSTs 風險源頭提供了重要依據。Lewis 等(2022)用蘇格蘭的鏈狀亞歷山大藻和英格蘭南部的鏈狀亞歷山大藻分別投喂紫貽貝測定其體內的毒素譜, 根據毒素譜差異成功區分了攝食兩種產毒藻的樣品。同時有研究(Yanetal, 2022)表明不同產毒藻的毒素組成及毒性大小存在著較大差別, 并且我國的毒素分布也具有區域性特征。由于雙殼貝類攝食產毒藻后毒素會在體內發生代謝轉化,如攝食鏈狀亞歷山大藻后的紫貽貝中毒素成分隨時間發生轉化, C2 占比減小而GTX5 占比增加(張海濤等, 2023), 蛤仔在攝食微小亞歷山大藻后毒素的積累階段, GTX1&4 的占比逐漸增加, 在排出階段占比下降( 焦玥 等, 2010), 這與本研究結果一致。在之前研究中發現這種代謝轉化會受到投喂產毒藻量(Chenetal,2001)、環境的酸堿性(Cheetal, 2020)等因素影響, 因此僅根據毒素信息無法進行準確分類。

研究表明, 含水的有機提取劑如甲醇水溶液(Sidwicketal, 2017)可以有效地提取極性代謝物, 而非極性代謝物常選擇氯仿、二氯甲烷等有機試劑(Caoetal,2020)進行提取。通過比較了5 種高效前處理方法, 經過扣除空白貽貝基質干擾后, 最終選擇了乙腈/甲醇/丙酮作為提取劑的提取方法, 可以有效地移除大分子蛋白質(Lietal, 2017), 提高弱極性和非極性物質(Mushtaqetal, 2014)的提取效率。通過比較, 特征物質的提取方法和色譜質譜方法差異系數較小, 表明儀器穩定較好(李倩倩等, 2021), 可為后續特異性物質的篩選提供了方法。采用PCA 分析發現, 三組樣品有明顯分離; 進一步采用OPLS-DA 分析, 剔除無關數據使得模型更加易于解釋, 可視化效果更加明顯(鐘若梅等, 2018)。劉鴿等(2021)采用PCA 與OPLSDA 分析, 通過VIP 值>1 以及P<0.05 篩選出暴露無機砷后三疣梭子蟹的100 種差異代謝產物。Vera 等(2019)等采用主成分分析(PCA)并比較獨立建模(SIMCA)和偏最小二乘判別分析(PLS-DA)作為化學計量學分類方法, 實現了對品種特級初榨橄欖油的地理來源高靈敏和無偏差的鑒定。進一步采用Fisher線性判別模型構建判別函數, 可在變量服從正態分布、沒有顯著相關或變量的平均值和方差不相關使效果較好(李洪成等, 2017)。歐陽建等(2022)通過對黃金茶綠茶風味的測定將52 個樣品分為醇厚型、鮮爽型和其他型3 種滋味類型。開建榮等(2022)構建枸杞原產地鑒別的判別模型, 模型的正確判別率達到了97.3%。目前, 未見Fisher 線性判別模型用于對攝食麻痹性貝類毒素產毒藻的溯源報道。本研究引入6 種物質構建的判別模型交叉驗證準確率達到了88.9%,預測效果較為理想。指紋識別技術需要結合識別和鑒定(Cuadros-Rodríguezetal, 2021)兩個過程, 為了達到更好地溯源目的還需要進一步的現場驗證。同時,由于在同一海域往往會存在多種產毒藻且毒素成分會更加復雜, 如遼東灣曾檢測出安氏亞歷山大藻、偽亞歷山大藻和李氏亞歷山大藻等(宋倫等, 2018), 長江口海域則為塔瑪亞歷山大藻和鏈狀亞歷山大藻(高巖等, 2016)。因此, 需持續采集典型藻株特征物質識別以增加溯源技術的準確度和適用性。

4 結論

通過優化前處理方法, 確定了甲醇/乙腈/丙酮(體積比為 1∶1∶1)作為提取劑; 采用 Q-Exactive Orbitrap 篩選 C8 色譜柱中正離子掃描模式m/z100~600 的掃描范圍和Amide 色譜柱正離子掃描模式m/z600~1 500 的掃描范圍, 可獲得40.4%的指紋信息。通過多元統計分析技術篩選出攝食不同產毒藻貽貝的13 種復合指紋物質, 構建了判別模型交叉驗證準確率達到了88.9%。研究成果有助于提高麻痹性貝類毒素區域風險源頭的準確識別能力, 為防范食用安全風險提供了科學依據和技術支撐。

猜你喜歡
貽貝貝類攝食
我國海水貝類養殖低碳效應評價
兩種不同投喂策略對加州鱸攝食量和生長的影響
貽貝、海虹、青口、淡菜……到底有何區別
QuEChERS-液相色譜-高分辨質譜法測定貝類中6種親脂性貝類毒素
鮮美貝類可能暗藏毒素
輕則攝食減慢,重則大量死魚!加州鱸養殖亞硝酸鹽超標,預防處理如何做好?
噪音太吵,貽貝受不了
籠養灰胸竹雞入冬前日攝食量規律的研究
全球貽貝產量減少導致貿易受挫
基于Ⅲumina平臺的厚殼貽貝外套膜轉錄組從頭測序
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合