盧梓雅，陳曉琳，黃思敏，段茜茜，彭博言，肖子健，項(xiàng)載盈，郭曉梅，劉雅琪，林思晴，TAN Karsoon，張洪寬，鄭懷平

（汕頭大學(xué)海洋科學(xué)研究院/廣東省海洋生物技術(shù)重點(diǎn)實(shí)驗(yàn)室/廣東省亞熱帶海水養(yǎng)殖工程技術(shù)研究中心，廣東汕頭 515063）

0 引言

牡蠣養(yǎng)殖在中國貝類養(yǎng)殖業(yè)中占有重要地位，是中國傳統(tǒng)四大養(yǎng)殖貝類之一，沿海各地都有養(yǎng)殖[1]。太平洋牡蠣(Crassostrea gigas)又稱長牡蠣、褶牡蠣，隸屬于軟體動物門(Mollusca)、瓣鰓綱(Lamellibranchia)、翼形亞綱(Pteriomorphia)、珍珠貝目(Pterioida)、牡蠣科(Ostreidae)、巨蠣屬(Crassostrea)。太平洋牡蠣原產(chǎn)于日本，20 世紀(jì)80 年代被中國引進(jìn)，因其生長快、產(chǎn)量高、口感細(xì)嫩、味鮮清香，很快成為中國主要海水養(yǎng)殖品種，帶動了中國牡蠣養(yǎng)殖業(yè)的發(fā)展[2]。2020年，中國牡蠣養(yǎng)殖產(chǎn)量達(dá)542.5萬t，位居世界首位[3]。牡蠣口感鮮嫩爽滑，富含優(yōu)質(zhì)蛋白且營養(yǎng)成分易被機(jī)體吸收，深受廣大消費(fèi)者喜愛[4]。因此，研究牡蠣響應(yīng)環(huán)境脅迫的分子機(jī)制具有重大的應(yīng)用和理論價值。

近年來，由于人類活動加劇，近海海洋生態(tài)系統(tǒng)的穩(wěn)定性被破壞[5]，其中最主要的就是海洋酸化。工業(yè)革命以來，人類活動釋放了包括CO2在內(nèi)的大量溫室氣體，導(dǎo)致全球變暖，海洋作為地球上最大的CO2存儲體，吸收了大量的CO2（超過100萬t/h），極大地緩解了全球變暖，但卻使海水pH 下降，最終導(dǎo)致海洋酸化[6]。據(jù)統(tǒng)計(jì)，自工業(yè)革命以來，海洋水體的pH 已下降了0.1。根據(jù)IPCC 預(yù)測模型(A1F1)推測，至2100 年表層海洋水體pH 將會下降0.3～0.4，氫離子濃度則增加100%～150%。海洋酸化正在逐漸改變海水的理化性質(zhì)、海洋生物的生理機(jī)能和海洋生物群落，進(jìn)而威脅到了海洋生態(tài)平衡。最近的研究表明，人為氣候變化的影響從根本上改變了海洋生態(tài)系統(tǒng)，包括海洋生產(chǎn)力下降、食物網(wǎng)動態(tài)變化、棲息地物種豐度減少以及物種分布變化等[7-8]。

海洋酸化對太平洋牡蠣的幼體和成體都有很大的危害。生活在海洋酸化環(huán)境下的牡蠣，雖然幼體牡蠣的生長速率要高于成體牡蠣，但是幼體牡蠣為了保障機(jī)體功能的正常運(yùn)行，往往會消耗更多的能量，并且長期處于該環(huán)境下，幼體的死亡率會更高[9]。與此同時，幼體牡蠣外殼和軟體部的生長情況在酸化條件下會比生活在正常環(huán)境下的差很多，且貝殼十分脆弱；對于成體牡蠣來說，一方面本身長速變緩，另一方面由于新陳代謝的速率比幼體慢，因此生長速率并沒有受到太大的影響[9]。但是在這種弱酸性的環(huán)境條件下，成體牡蠣的外套膜增加了碳酸酐酶的表達(dá)，碳酸酐酶是一種專門調(diào)節(jié)pH的酶，有助于殼體內(nèi)碳酸氫鹽的生成，這種酶含量的增加，表明成體牡蠣能夠在一定程度上平衡酸化條件造成的機(jī)體損傷[9]。

目前研究太平洋牡蠣的學(xué)者們雖然從不同方面展示了海洋酸化對牡蠣生長發(fā)育的影響，但關(guān)于牡蠣響應(yīng)海洋酸化脅迫下的分子機(jī)制尚未深入研究。因此，本研究利用公共數(shù)據(jù)庫NCBI 中的RNA-seq 數(shù)據(jù)（包括4 個組：pH 7.8、7.4、7.0、6.6）進(jìn)行分析，進(jìn)而在分子層面探究太平洋牡蠣對海洋酸化響應(yīng)的機(jī)制。

1 材料與方法

1.1 實(shí)驗(yàn)材料

從NCBI 網(wǎng)站的SRA 數(shù)據(jù)庫(Sequence Read Archive)上下載了太平洋牡蠣在pH 6.6、7.0、7.4、7.8的轉(zhuǎn)錄組數(shù)據(jù)(https://www.ncbi.nlm.nih.gov/bioproject/PRJNA735889)以及Genome 數(shù)據(jù)庫中太平洋牡蠣的參考基因組(https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/902/806/645/GCF_902806645.1_cgigas_uk_roslin_v1/GCF_902806645.1_cgigas_uk_roslin_v1_genomic.fna.gz)，本研究于2022 年7—8 月在汕頭大學(xué)海洋生物研究所開展。

1.2 方法

1.2.1 數(shù)據(jù)格式轉(zhuǎn)換 利用SRAtoolkit 2.11.3 軟件的fastp-dump 將SRA 格式的文件轉(zhuǎn)化為fastq.gz 格式的文件。

1.2.2 數(shù)據(jù)質(zhì)控利用Fastp 軟件對轉(zhuǎn)換格式后的數(shù)據(jù)進(jìn)行質(zhì)控，過濾低質(zhì)量的reads，去除adapter，數(shù)據(jù)截取，最后生成HTML和JSON格式的質(zhì)控報告，統(tǒng)計(jì)過濾前和過濾后Q20、Q30及GC含量數(shù)據(jù)。

1.2.3 數(shù)據(jù)的比對首先建立太平洋牡蠣參考基因組的索引，然后通過Hisat2 將質(zhì)控得到的數(shù)據(jù)比對到參考基因組，得到sam格式的輸出文件和報告，統(tǒng)計(jì)數(shù)據(jù)的整體比對率(overall alignment rate)，唯一比對率(Unique mapping)，多重比對率(Multiple mapping)等信息。

1.2.4 數(shù)據(jù)排序由于從Hisat2得到的sam格式的文件占用存儲較大，利用Samtools 將sam 文件轉(zhuǎn)換為bam文件，并且對bam文件進(jìn)行排序。

1.2.5 數(shù)據(jù)的定量分析與標(biāo)準(zhǔn)化對于處理后得到的bam 文件，通過使用R 語言中的FeatureCounts 來計(jì)算readscounts，并轉(zhuǎn)化為TPM值(transcripts per million)。

1.2.6 差異表達(dá)分析通過R語言中的tidyverse包進(jìn)行多表關(guān)聯(lián)，將每組的差異表達(dá)基因篩選出來（基因差異表達(dá)的參數(shù)設(shè)置為|log2FoldChange|≥1和padj＞0.05），然后利用EnhancedVolcano 和pheatmap 包生成火山圖和熱圖，將差異表達(dá)基因可視化。

1.2.7 富集分析（GO 和KEGG）通過基迪奧生物信息云平臺(https://www.omicshare.com）進(jìn)行GO 分析(Gene Ontology)和KEGG 分析(Kyoto Encyclopedia of Genes and Genomes)。利用GO數(shù)據(jù)庫，通過計(jì)算差異表達(dá)基因參與生物學(xué)過程、細(xì)胞組分、分子功能的超幾何分布，進(jìn)行功能分類注釋和富集分析。同時，用KEGG數(shù)據(jù)庫獲得代謝通路注釋信息。

2 結(jié)果

2.1 測序數(shù)據(jù)的統(tǒng)計(jì)

對4 個組別pH 6.6、pH 7.0、pH 7.4、pH 7.8 的數(shù)據(jù)分別進(jìn)行質(zhì)控，獲得高質(zhì)量的測序數(shù)據(jù)。在未過濾前，這些數(shù)據(jù)的Q20 值最低為90.65%，最高為98.22%；Q30 值最低為84.59%，最高為98.17%，而GC 含量在42.04%～94.97%；過濾后Q20最低值達(dá)到96.82%，而最高值達(dá)到98.54%；Q30 值最低達(dá)到92.11%，最高達(dá)到95.34%，而GC 含量在42.33%～94.67%，表明測序數(shù)據(jù)經(jīng)過fastp質(zhì)控后質(zhì)量得到提升。利用將hisat2軟件將質(zhì)控后的數(shù)據(jù)與建好的牡蠣基因組索引文件進(jìn)行比對。測序reads 的總比對率值最低為76.26%，最高為79.36%（見表1），數(shù)據(jù)的唯一對照率最低為71.31%，最高為74.60%；多重對照率最低為4.23%，最高為6.12%，說明質(zhì)控后的數(shù)據(jù)質(zhì)量較高，可用于研究分析。

表1 測序數(shù)據(jù)比對到參考基因組的統(tǒng)計(jì)結(jié)果

2.2 差異表達(dá)基因的統(tǒng)計(jì)與分析

通過使用R語言中的tydiverse包對3組樣本進(jìn)行了差異表達(dá)基因分析，將log2FoldChange ＜1 且padj ＞0.05 的基因篩選出來，生成火山圖和熱圖（圖1 和圖2）。分別對pH 6.6 vs pH 7.8 組、pH 7.0 vs pH 7.8 組和pH 7.4 vs pH 7.8 組的顯著差異表達(dá)基因進(jìn)行分析，發(fā)現(xiàn)pH 6.6 vs pH 7.8 組總共有943 個顯著差異表達(dá)基因，其中上調(diào)基因有295個，下調(diào)基因有648個；pH 7.0 vs pH 7.8組總共有93個顯著差異表達(dá)基因，其中上調(diào)基因有46 個，下調(diào)基因有47 個；pH 7.4 vs pH 7.8 組總共有61個顯著差異表達(dá)基因，其中上調(diào)基因有31個，下調(diào)基因有30 個。通過火山圖將差異表達(dá)基因可視化（圖1）。

圖1 牡蠣各組差異表達(dá)基因火山圖

圖2 牡蠣各組差異表達(dá)基因熱圖

利用韋恩圖顯示pH 6.6 vs pH 7.8 組、pH 7.0 vs pH 7.8 組、pH 7.4 vs pH 7.8 組間共有及特有的基因數(shù)（圖3），響應(yīng)酸脅迫的基因總共有1059 個。當(dāng)pH 6.6 vs pH 7.8 時，基因個數(shù)已有943 個，pH 7.0 vs pH 7.8時，基因個數(shù)為93個；pH 7.4 vs pH 7.8時，基因個數(shù)達(dá)到61個；其中有4個基因在pH 6.6、pH 7.0和pH 7.4均差異表達(dá)。由此可知，隨著pH的下降，牡蠣響應(yīng)酸化脅迫的基因數(shù)也在增加。

圖3 太平洋牡蠣酸脅迫差異表達(dá)基因韋恩圖

利用韋恩圖可視化pH 6.6 vs pH 7.8 組、pH 7.0 vs pH 7.8組、pH 7.4 vs pH 7.8組間的差異基因。

2.3 GO功能分析

由于pH 7.0 vs pH 7.8 組和pH 7.4 vs pH 7.8 組的差異基因較少，因此，本研究僅對pH 6.6 vs pH 7.8 組的差異表達(dá)基因進(jìn)行了GO(Go Ontology)富集分析。將該組的943個差異表達(dá)基因在基迪奧生物信息云平臺(https://www.omicshare.com)上進(jìn)行GO 富集分析，對該組中差異表達(dá)基因的前53個GO term進(jìn)行列舉，發(fā)現(xiàn)其中8個GO term差異富集明顯。

在pH 6.6 vs pH 7.8組中的差異表達(dá)基因，在生物過程方面，大部分富集到了細(xì)胞過程、代謝過程、生物調(diào)節(jié)等過程中，如糖代謝、脂類代謝；在細(xì)胞組分方面，細(xì)胞器成分和細(xì)胞成分等組分被集中富集；在分子功能方面，以催化和分子結(jié)合為主（圖4）。

圖4 差異表達(dá)基因的GO統(tǒng)計(jì)柱狀圖

3 討論

（1）代謝方面

在海洋酸化脅迫下，很多海洋貝類和珊瑚等海洋生物的新陳代謝活動受到了影響[10]。紫貽貝(Mytilus edulis)在海水酸化的影響下能量代謝出現(xiàn)抑制現(xiàn)象[11]；淺珊瑚(Mesophotic corals)的代謝率在海洋酸化脅迫下顯著降低[12]。通過研究發(fā)現(xiàn)（見圖5），太平洋牡蠣在海洋酸化條件下的糖代謝、氨基酸代謝、脂類物質(zhì)代謝等都受到抑制。糖類和脂肪都是太平洋牡蠣的主要能源物質(zhì)，當(dāng)這些主要能源物質(zhì)的代謝被抑制時，會減少對太平洋牡蠣的能量供應(yīng)從而對太平洋牡蠣的生長、發(fā)育、繁殖等正?；顒釉斐刹焕绊慬13]。在研究中發(fā)現(xiàn)，糖代謝方面的下調(diào)差異基因數(shù)目(8)多于上調(diào)差異基因數(shù)目(3)，表明太平洋牡蠣的糖代謝由于受到海水酸化的脅迫而被抑制。在海洋酸化條件下(pH 6.6)，太平洋牡蠣的糖代謝差異表達(dá)基因gene-LOC105328227（與糖原磷酸化酶有關(guān)）的下調(diào)倍數(shù)最大，對太平洋牡蠣的糖代謝抑制起到主導(dǎo)作用。糖原磷酸化酶是關(guān)于糖原分解的酶。在糖原分解過程中，糖原磷酸化酶催化糖原的磷酸解作用，從而實(shí)現(xiàn)給太平洋牡蠣提供能量。但是在海水pH 逐漸降低的情況下，gene-LOC105328227（與糖原磷酸化酶有關(guān)）的表達(dá)量逐漸減少，甚至趨向于0（見圖6），對太平洋牡蠣的糖代謝極其不利。海水酸化通過抑制糖原磷酸化酶對糖原的分解能力來抑制太平洋牡蠣的糖代謝，進(jìn)而減少由糖原提供的能量。這將導(dǎo)致太平洋牡蠣在海洋酸化脅迫下將供給能量的方式由糖原代謝轉(zhuǎn)化為蛋白質(zhì)消耗，這將嚴(yán)重?fù)p耗它們自身的肌肉組織等，從而影響到正常的生命活動[14]。厚殼貽貝和紫貽貝在海水酸化下出現(xiàn)類似現(xiàn)象，它們在海水pH 逐漸下降的情況下，代謝底物由碳水化合物轉(zhuǎn)化為蛋白質(zhì)，出現(xiàn)糖代謝抑制現(xiàn)象[8]。在研究中發(fā)現(xiàn)，在海洋酸化脅迫下，太平洋牡蠣的新陳代謝快慢受到氨基酸代謝的影響，氨基酸代謝的下調(diào)基因數(shù)目遠(yuǎn)多于上調(diào)基因數(shù)目，表明太平洋牡蠣的氨基酸代謝受到抑制。其下調(diào)基因gene-LOC117691793和gene-LOC10532916的下調(diào)倍數(shù)顯著，最高達(dá)7 倍，并且這2 個基因都與精氨酸激酶有關(guān)。精氨酸激酶是無脊椎動物調(diào)節(jié)能量代謝的最重要的酶之一，在細(xì)胞活動活躍時，精氨酸激酶會促進(jìn)高能磷酸鍵轉(zhuǎn)移并形成ATP[15]。在海洋酸化下，gene-LOC117691793和gene-LOC10532916的表達(dá)量相較于在中性海水嚴(yán)重減少（見圖6），很大程度上抑制形成精氨酸激酶，以至于太平洋牡蠣形成的ATP 減少，最終導(dǎo)致呼吸代謝等減弱[16]。

圖6 各個差異表達(dá)基因在不同pH下的表達(dá)量

在研究中發(fā)現(xiàn)，脂肪代謝的下調(diào)基因數(shù)目(6)遠(yuǎn)高于上調(diào)基因數(shù)目(1)，其中以基因gene-LOC105318687差異倍數(shù)為最高。gene-LOC105318687與脂肪酸合成酶(Fatty Acid Synthase,FAS,E.C.)相關(guān)，屬于脂質(zhì)代謝基因，通過影響脂質(zhì)代謝過程的效率，從而影響組織脂肪的含量[17]。脂肪酸合成酶又稱脂肪酸合酶，是內(nèi)源性脂肪酸合成的關(guān)鍵酶。脂肪酸合成酶能夠催化丙二酰輔酶A等酶，在體內(nèi)合成擁有16碳的長鏈飽和軟脂酸鹽。海洋酸化脅迫下，其表達(dá)水平顯著下降（見圖6），表明其酶的含量也顯著下降，這將影響太平洋牡蠣脂肪組織的發(fā)育、脂肪的沉積，進(jìn)而打破脂肪合成代謝與分解代謝之間的平衡狀態(tài)，影響牡蠣供能物質(zhì)的攝入量與其體內(nèi)能量消耗的平衡，進(jìn)而導(dǎo)致能量代謝的紊亂[18]。在其他水產(chǎn)動物中也有相似的研究，比如通過在飼料中添加膽汁酸，齊口裂腹魚(Schizothorax prenanti)降低其體內(nèi)的脂肪酸合成酶的表達(dá)量，從而控制翻譯水平下降，減少體脂合成，減少脂肪沉積[19]；在短期酸化脅迫10 d，紫貽貝(Mytilus edulis)選擇降低攝食，減少吸收的能量，增加呼吸代謝能量消耗，從而減少貝類生長的能量，使其生長下降[20]。

（2）免疫方面

在海洋酸化的影響下，不僅導(dǎo)致太平洋牡蠣代謝功能的減弱，其免疫功能也被顯著抑制。通過比較分析發(fā)現(xiàn)，在pH 6.4組和pH 7.8組相比時，免疫代謝方面下調(diào)表達(dá)數(shù)目(12)多于上調(diào)表達(dá)基因(9)，表明太平洋牡蠣在酸化脅迫下免疫作用受到抑制。貝類作為無脊椎動物，具有不同于脊椎動物的免疫機(jī)制來抵抗外界病原體的入侵，其宿主防御包括細(xì)胞和體液系統(tǒng)，血細(xì)胞可以清除病原微生物和死亡細(xì)胞，參與損傷修復(fù)，是貝類主要的免疫防御機(jī)制[21]。當(dāng)太平洋牡蠣的免疫功能受到抑制時，太平洋牡蠣的生長發(fā)育等生理活動都會受到影響。在研究中發(fā)現(xiàn)，在海洋酸化脅迫下(pH 6.6)，太平洋牡蠣起免疫防御功能的差異表達(dá)基因gene-LOC105345760（CD151 抗原有關(guān)）下調(diào)倍數(shù)最大，表明該基因可能對其在海洋酸化條件下的免疫功能起著關(guān)鍵作用。CD151 蛋白是一種細(xì)胞表面糖蛋白，參與多種病理生理過程，包括細(xì)胞的黏附、遷移、半橋粒結(jié)構(gòu)的形成以及通過囊泡進(jìn)行整合素的運(yùn)輸?shù)萚22-23]，因此，CD151基因的表達(dá)量減少可能導(dǎo)致太平洋牡蠣免疫系統(tǒng)的發(fā)生一定程度的受損。此外，從圖6 中可以看出，在海水pH 逐漸降低時，CD151的表達(dá)量不斷減少，甚至趨向于0，表明海水酸化顯著抑制了CD151抗原的表達(dá)，這可能會導(dǎo)致太平洋牡蠣體內(nèi)滲透壓改變，免疫系統(tǒng)紊亂，降低太平洋牡蠣抵擋外界環(huán)境變化的能力。研究表明，酸化(p CO2101～203 Pa)可以引起太平洋牡蠣(C.gigas)血細(xì)胞凋亡和活性氧生成的增加，并且升高的p CO2對一些抗氧化酶活性具有抑制作用，導(dǎo)致消化腺中的谷胱甘肽水平降低，過氧化氫酶顯著下降[24-25]。泥蚶(Tegillarca granosa)在酸化（pH 7.4 和7.8）脅迫下的總血淋巴細(xì)胞數(shù)和吞噬率顯著減少，嗜酸性和嗜堿性粒細(xì)胞的百分比分別顯著減少和增加，信號傳導(dǎo)也受到了一定的負(fù)面影響[26]。由此可知，海洋酸化對貝類的免疫系統(tǒng)有顯著的影響，能夠顯著地抑制雙殼貝類的免疫應(yīng)答作用[27]。

綜上所述，本研究通過轉(zhuǎn)錄組分析獲得太平洋牡蠣在不同pH(7.8、7.4、7.0、6.6)環(huán)境下的差異表達(dá)基因，發(fā)現(xiàn)在海水酸化(pH 6.6)時，太平洋牡蠣的代謝功能和免疫功能都受到顯著抑制。在代謝方面中，糖代謝相關(guān)的gene-LOC105328227（糖原磷酸化酶）、氨基酸代謝相關(guān)的gene- LOC117691793和gene-LOC10532916（精氨酸激酶）和脂類物質(zhì)代謝相關(guān)的gene-LOC105318687（脂肪酸合成酶）表達(dá)量都隨著pH 的下降而不斷減少，有的甚至趨于零，表明在海洋酸化條件下太平洋牡蠣的代謝功能受到嚴(yán)重的抑制；在免疫方面中，gene-LOC105345760（CD151 抗原）的表達(dá)量顯著下降，表明海洋酸化脅迫很有可能會導(dǎo)致太平洋牡蠣體內(nèi)滲透壓改變，免疫系統(tǒng)紊亂，降低其抵擋外界環(huán)境變化的能力。面對全球海洋酸化的大趨勢下，希望本研究成果能夠彌補(bǔ)海洋貝類響應(yīng)海洋酸化分子機(jī)制研究中的不足，能夠?yàn)楹Ｑ筘愵惖倪z傳改良提供一定的理論基礎(chǔ)。

4 結(jié)論

本研究利用生物信息學(xué)對海洋酸化脅迫下的太平洋牡蠣轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行了分析，結(jié)果顯示pH 6.6 vs pH 7.8 組總共有943 個顯著差異表達(dá)基因，pH 7.0 vs pH 7.8組總共有93個顯著差異表達(dá)基因pH 7.4 vs pH 7.8 組總共有61 個顯著差異表達(dá)基因。對pH 6.6 vs pH 7.8組的差異表達(dá)基因進(jìn)行GO和KEGG分析發(fā)現(xiàn)太平洋牡蠣在海洋酸化脅迫下其代謝和免疫功能受到顯著抑制，表明海洋酸化極可能導(dǎo)致牡蠣機(jī)體代謝和免疫系統(tǒng)的紊亂，降低太平洋牡蠣抵擋外界環(huán)境變化的能力。