董忠典 黎學友 廖 健 張 寧 郭昱嵩 王中鐸

(廣東海洋大學水產(chǎn)學院 南海水產(chǎn)經(jīng)濟動物增養(yǎng)殖重點實驗室 湛江 524088)

肝臟是脊椎動物最重要的器官之一， 在物質(zhì)能量代謝、激素合成及免疫反應(yīng)等生理過程中具有重要意義(van der Oostet al， 2003； Tomet al， 2004； Laiet al， 2015； Zhanget al， 2017)。在魚類中， 雌性肝臟還是卵黃蛋白原(vitellogenins，vtgs)、透明帶蛋白前體(choriogenins，chgs)、雌激素受體(estrogen receptor，er)和細胞色素P450(cytochrome P450s，cyp450s)等基因的主要表達器官， 這些基因?qū)Νh(huán)境中雌激素物質(zhì)較為敏感， 當水體環(huán)境中存在雌激素類物質(zhì)時， 可以誘導(dǎo)雄魚肝臟表達這些基因(Leeet al， 2002； Chenet al， 2008； Laiet al， 2015)。

弓背青鳉(Oryzias curvinotus)是一種廣泛分布于我國南海沿岸水域的小型魚類， 為紅樹林水域的常見物種(Hayakawaet al， 2015； Wanget al， 2017)。弓背青鳉具有體型小、生長快、性成熟時間短和胚胎發(fā)育透明的特點， 可以用于實驗室飼養(yǎng)繁殖； 其對鹽度有較強的適應(yīng)能力， 具備開發(fā)成海洋環(huán)境毒理及鹽度適應(yīng)研究新模式種的潛力(Hayakawaet al， 2015)。同時， 弓背青鳉的性染色體為XX/XY 型， 具有性別標記， 是研究海洋魚類性別決定和分化機制的潛在模型(Matsudaet al， 2003； 董忠典等， 2018)。同屬的日本青鳉(Oryziaslatipes)和海水青鳉(Oryziasmelastigma)已被廣泛應(yīng)用于淡水和海水環(huán)境檢測研究(Metcalfeet al， 2000； Chenet al， 2008； Honget al， 2015； Conget al， 2017； Abdel-moneimet al， 2018)。然而， 缺乏基因組和轉(zhuǎn)錄組信息， 在很大程度上限制了弓背青鳉資源的開發(fā)和利用。本研究采用RNA-Seq 技術(shù)分別對雌、雄弓背青鳉肝臟進行了轉(zhuǎn)錄組測序、基因功能注釋和表達譜分析， 以期為弓背青鳉基因資源的開發(fā)提供基礎(chǔ)數(shù)據(jù)。

1 材料與方法

1.1 弓背青鳉的采集及總肝臟RNA 提取

本研究所用性成熟弓背青鳉(平均體長29.25mm， 平均體質(zhì)量0.27g)采集自湛江雷州半島近岸的高橋國家級紅樹林保護區(qū)(21°36′24"N， 109°47′8"E)， 實驗室 暫 養(yǎng) 30d， 溫 度 25±1°C， 鹽 度 15， 光 暗 周 期14h：10h。解剖通過性腺類型確定弓背青鳉生理性別， 雌雄各 3 尾取肝臟分別進行總 RNA 提取， 使用DNase I(寶生物， 大連)去除基因組污染。使用1%瓊脂糖凝膠電泳分析RNA 完整性， Nanodrop 檢測RNA的純度(OD260/OD280比值)， Qubit 檢測總RNA 濃度， RNA 滿足進行下一步的文庫構(gòu)建的要求。對3 尾雌性肝臟總RNA， 分別取等量(1.0μg)混合為一個RNA樣品池(FL)， 用同樣的方法獲得雄性肝臟RNA 樣品池(ML)， 對FL 和ML 各取1.5μg 總RNA， 分別進行cDNA 文庫構(gòu)建和轉(zhuǎn)錄組測序。

1.2 文庫構(gòu)建及測序

使用帶Oligo(dT)的磁珠富集真核生物mRNA， 加入fragmentation buffer 將mRNA 打斷成短片段， 以mRNA 為模板， 用六堿基隨機引物合成第一鏈cDNA， 然后加入緩沖液、dNTPs 和DNA 聚合酶和RNase H合成第二鏈cDNA， 再用AMPure XP beads 純化雙鏈cDNA。純化的雙鏈cDNA 先進行末端修復(fù)、加A 尾并連接測序接頭， 再用 AMPure XP beads 選擇150—200bp 片段。隨后進行PCR 擴增， 并用AMPure XP beads 純化PCR 產(chǎn)物， 得到最終的文庫。文庫構(gòu)建完成后， 用Qubit2.0 進行初步定量并稀釋至1.5ng/μL， 用 Agilent 2100 確認文庫插入片段大小， 最后用Q-PCR 方法對文庫的有效濃度進行準確定量(文庫有效濃度＞2nmol/L)， 以保證文庫質(zhì)量。庫檢合格后， 進行Illumina HiSeq 測序， cDNA 文庫的制備和測序由諾禾致源生物信息科技有限公司(天津)完成。

1.3 轉(zhuǎn)錄組測序質(zhì)量控制和組裝

為了保證轉(zhuǎn)錄組分析質(zhì)量， 對原始測序數(shù)據(jù)(raw reads)過濾， 去除帶接頭的reads、N(N 表示無法確定堿基信息)的比例大于0.1%的reads、低質(zhì)量(質(zhì)量值Qphred＜=20 的堿基數(shù)占整個reads 的50%以上)的reads， 得到高質(zhì)量測序數(shù)據(jù)(clean reads)。采用Trinity軟件對clean reads 進行拼接(Grabherret al， 2011)。

1.4 基因功能注釋

為獲得全面的基因功能信息， 使用7 大數(shù)據(jù)庫對Trinity 拼接獲得的unigenes 進行功能注釋。使用在線BLAST 程序?qū)nigenes 在Nr (http：//www.ncbi.nlm. nih.gov)期望值e-value ＜ e-5、Nt (http：//www.ncbi.nlm. nih.gov) e-value ＜ e-5、Swiss-Prot (http：//www.expasy. ch/sprot) e-value ＜ e-5、KOG (http：//www.ncbi.nlm.nih. gov) e-value ＜ e-3、KO (https：//www.kegg.jp/kegg/ko. html) e-value ＜ e-10和GO (http：//www.geneontology. org/) e-value ＜ e-6進行比對注釋。通過HMMER3.0 (e-value ＜ 0.01)程序?qū)nigenes 在Pfam 數(shù)據(jù)庫中進行蛋白家族注釋(Eddy， 2011)。

1.5 基因表達水平和差異分析

將Trinity 拼接得到的轉(zhuǎn)錄組作為參考序列， 通過RSEM 對每個樣本的基因表達水平進行估計： 將Clean data 比對到組裝好的參考序列上， 根據(jù)比對結(jié)果得到每個基因的Readcount 數(shù)目， 采用TMM 對read count 數(shù)據(jù)進行標準化處理， 再用DEGseq 進行差異分析， 篩選閾值為qvalue ＜0.005 且|log2(foldchange)|＞1 (Wanget al， 2010)。

1.6 差異基因GO 和KEGG 富集分析

使用R 軟件中的GOseq 對弓背青鳉雌雄肝臟轉(zhuǎn)錄組差異基因進行顯著性富集分析(Younget al， 2010)。使用KOBAS 軟件檢測差異基因在KEGG 通路中富集情況(Maoet al， 2005)。

1.7 實時熒光定量PCR(RT-qPCR)驗證

利用RT-qPCR 技術(shù)， 弓背青鳉雌雄肝臟轉(zhuǎn)錄組中選擇 12 個基因進行定量驗證(引物見表 1)。RT-qPCR 反應(yīng)體系15μL， 包含7.5μL 2×Power Green qPCR Mix (東盛生物， 廣州)， 0.6μL 上下游引物(10μmol/L)， 1.5μL cDNA， 4.8μL ddH2O。每個樣品技術(shù)重復(fù)3 次， 反應(yīng)置于Roche LightCycler 96(羅氏， 瑞士)上運行。反應(yīng)程序如下： 95°C 3min； 95°C 15s， 60°C 15s， 72°C 30s (采集熒光)， 40 個循環(huán)； 熔解曲線分析檢測PCR 產(chǎn)物的特異性。gapdh和adp作為內(nèi)參基因， 使用SPSS17.0 采用2-ΔΔCt方法對目的基因表達水平進行統(tǒng)計(Livaket al， 2001)。

表1 本研究所用引物及序列 Tab.1 The sequences of primers used in this study

續(xù)表

1.8 微衛(wèi)星標記(SSR)檢測

采 用 MISA (http：//pgrc.ipk-gatersleben.de/misa/ misa.html)對弓背青鳉肝臟轉(zhuǎn)錄組拼接獲得的unigenes進行微衛(wèi)星標記(SSR)檢測， 設(shè)置各重復(fù)單元類型的最少重復(fù)次數(shù)為： 1—10、2—6、3—5、4—5、5—5、6—5。以1—10 為例， 該設(shè)置表示單核苷酸重復(fù)類型至少重復(fù)10 次才被算為微衛(wèi)星。

2 結(jié)果

2.1 測序結(jié)果與組裝

本研究分別構(gòu)建了雌、雄弓背青鳉肝臟cDNA 文庫， 并進行了轉(zhuǎn)錄組測序。測序結(jié)果經(jīng)過質(zhì)控后， 分別獲得80095044 和87043984 條Clean reads。經(jīng)Trinity 組裝， 共得到49912 個unigenes， 其中N50 長度為2394bp， 其中長度大于1000bp 的unigenes 共計25423 個(圖1)。

圖1 unigenes 序列長度分布圖 Fig.1 The distribution of sequence length of the unigenes

2.2 轉(zhuǎn)錄組注釋

通過BLAST 比對， 分別在NR、NT、SwissProt、PFAM、KOG 數(shù)據(jù)庫中注釋到32651、47988、19005、25587 和19005 個unigenes。與NR 數(shù)據(jù)庫比對發(fā)現(xiàn)有32649 個unigenes 與390 個物種基因序列高度同源， 其中77.6%(25349 個)的unigenes 和日本青鳉基因序列同源， 其次有 3.6%的 unigenes 和深裂眶鋸雀鯛(Stegastes partitus)同源(圖2)。為進一步了解弓背青鳉雌、雄肝臟表達基因集的功能， 將拼接的unigenes在GO 和KO 數(shù)據(jù)庫中進行了比對分析。25630 個unigenes 被歸類到GO 三個大類(生物學過程、分子功能和細胞成分)55 個功能分類中(圖3)。生物學過程包 括25 個功能分類， 參與“細胞過程”(15094)和“代謝過程”(12821)的基因最多； 細胞成分包含20 個功能分類，“細胞”(8419)和“細胞部分”(8419)基因最多； 分子功能包含10 個功能分類， “結(jié)合”(16155)和“催化活性”(10114)基因最多。對基因做KO 注釋后， 可根據(jù)它們參與的KEGG 代謝通路進行分類。共有18887個unigenes 在KO 數(shù)據(jù)庫中獲得注釋， 富集到232 個KEGG 信號通路中， 其中“PI3K-Akt”信號通路含有最多(657)的 unigenes， 其次是“內(nèi)吞”(615)、“斑黏連”(521)和“肌動蛋白骨架”的調(diào)節(jié)通路(508)。

圖2 弓背青鳉unigenes 同源物種分布 Fig.2 Distribution of homologous species of O. curvinotus unigenes

2.3 差異基因表達及富集分析

通過DEGseq 對弓背青鳉雌、雄肝臟轉(zhuǎn)錄組基因進行差異表達分析(差異基因篩選條件為： qvalue ＜ 0.005 & |log2(foldchange)| ＞ 1)， 結(jié)果如圖4 所示， 有207 個unigenes 在雌魚肝臟中上調(diào)， 364 個unigenes 在雄魚肝臟上調(diào)。對差異基因(DEGs)進行GO 和KEGG 富集分析， GO 富集結(jié)果顯示雌性肝臟高表達基因主要 參與蛋白質(zhì)合成過程， 雄性高表達基因主要參與免疫及氧化還原反應(yīng)(圖5)。KEGG 分析結(jié)果顯示DEGs 富集到244 個KEGG 通路中， 涉及核糖體、PPAR 信號通路、卵巢類固醇激素生成、類固醇生物合成、不飽和脂肪酸合成和內(nèi)質(zhì)網(wǎng)蛋白質(zhì)加工等通路(圖6)。

2.4 RT-qPCR 驗證

本研究選擇了在雌、雄弓背青鳉肝臟RNA-Seq分析中的12 個差異表達基因， 通過RT-qPCR 技術(shù)對RNA-Seq 結(jié)果進行驗證。所選的驗證基因包括環(huán)腺苷三磷酸依賴的轉(zhuǎn)錄因子camp， 內(nèi)質(zhì)網(wǎng)基因endoplasmin， UDP-葡萄糖醛酸轉(zhuǎn)移酶基因utp， 磺基轉(zhuǎn)移酶家族2B 基因sult2b及部分檢測環(huán)境雌激素物質(zhì)的標記基因如卵黃蛋白原基因vtg1、vtg2、vtglike， 透明帶蛋白前體基因chgl、chgh、chghm， 細胞色素P450 基因p4502k1、p45027c。結(jié)果表明所選基因的表達模式與RNA-Seq 分析相一致， 說明RNA-Seq 分析結(jié)果可信(見圖7)。

圖6 弓背青鳉雌雄肝臟差異基因KEGG 富集 Fig.6 KEGG assignment of differentially expressed genes of O. curvinotus

圖7 RT-qPCR 檢測RNA-Seq 結(jié)果 Fig.7 Validation of RNA-Seq data using RT-qPCR

2.5 簡單重復(fù)序列SSR 分析

利用MISA 軟件對弓背青鳉進行SSR 標記鑒定， 用于遺傳多樣性研究。有18326 條unigenes 序列含有SSR 標記， 共計28195 個， 其中6373 條unigenes 包含一個以上的SSR。各種類型的SSR 標記出現(xiàn)的頻率不同(表2)， 單核苷酸重復(fù)SSR 出現(xiàn)的頻率最高， 占總SSR 的73.17%， 其次是三核苷酸重復(fù)SSR， 占總SSR 的19%。另外， 不同重復(fù)序列的相對豐度差異很大。A/T 在單核苷酸SSR 中最常見， 在雙核苷酸SSR中， AG/CT 最為常見。AGG/CCT 和AAAC/GTTT 是三、四核苷酸SSR 中最常見的元件。

3 討論

3.1 弓背青轉(zhuǎn)錄組測序分析

為了闡明弓背青鳉肝臟表達基因的功能和涉及的生物學過程， 本研究分別構(gòu)建了性成熟雌、雄弓背青鳉肝臟 cDNA 文庫， 并在 Illumina 平臺進行了RNA-Seq 測序， 分別獲得了80095044 和87043984 條Clean reads。N50 長度是評價RNA-Seq 組裝質(zhì)量的重要參數(shù)， 本研究弓背青鳉肝臟轉(zhuǎn)錄組N50長度為2394bp， 高于多個已發(fā)表的水生動物肝臟轉(zhuǎn)錄組， 包括海水青鳉(2162bp) (Laiet al， 2015)， 七鰓鰻(lampetra japonica)(1447bp) (李 慶 偉 等， 2018)， 脂 鯉(Pterygoplichthys anisitsi) (1571bp) (Parenteet al， 2017)， 低 于 海 鱸(Dicentrarchus labrax) (3257bp) (Magnanouet al， 2014)。對NR 數(shù)據(jù)庫比對， 發(fā)現(xiàn)有77.64% (32649 條unigenes 中的25349 個)的可比對序列與日本青鳉基因同源(圖2)， 可能是因為弓背青鳉和日本青鳉均屬于青鳉屬， 親緣關(guān)系密切(Wanget al， 2017)。

表2 弓背青鳉肝臟轉(zhuǎn)錄組SSR 類型及數(shù)量 Tab.2 SSR summary from the O. curvinotus liver transcriptome

肝臟在魚類生長和物質(zhì)能量代謝中起著重要的作用。在弓背青鳉肝臟轉(zhuǎn)錄組中， “細胞過程”、“細胞”和“結(jié)合”分別在GO 三大功能類別中富集(圖3)。這與其他水生動物的肝臟RNA-Seq 結(jié)果類似， 說明屬于這幾類功能的基因在物種間保守(Magnanouet al， 2014； Chenet al， 2016； 張燕萍等， 2018)。KEGG 的注釋結(jié)果表明， 弓背青鳉肝臟轉(zhuǎn)錄組中unigenes 主要參與了信號轉(zhuǎn)導(dǎo)(3130)、內(nèi)分泌系統(tǒng)(1536)及免疫系統(tǒng)(1491)的生物學過程。結(jié)果可能提示了肝臟在弓背青鳉生命過程中的主要功能， 為了解肝臟的功能演變奠定了基礎(chǔ)。

3.2 弓背青 免疫相關(guān)基因

肝臟是魚類的重要免疫器官， 通過KEGG 分析， 從弓背青鳉肝臟轉(zhuǎn)錄組中篩選出了多個可能參與免疫相關(guān)防御途徑的unigenes(表3)： 373 個unigenes 參與了血小板活化， 272 個unigenes 參與了白細胞跨內(nèi)皮遷移， 245 個unigenes 參與了T 細胞受體信號通路， 136 個unigenes 參與了補體和凝血級聯(lián)， 188 個unigenes 參與了自然殺傷細胞介導(dǎo)的細胞毒性， 163 個unigenes 參與了Toll-樣受體信號通路等過程。在七鰓鰻的肝臟轉(zhuǎn)錄組中也有相似的結(jié)果(李慶偉等， 2018)。

表3 弓背青鳉肝臟轉(zhuǎn)錄組中參與免疫反應(yīng)相關(guān)KEGG 通路的unigenes Tab.3 KEGG pathways with immune-related genes enrichment in the liver transcriptomes of O. curvinotus

血小板上具有toll 樣受體， 在機體的先天免疫應(yīng)答中起著關(guān)鍵作用。當病原菌與血小板結(jié)合時， 血小板活化同時分泌抗菌肽， 引起機體免疫反應(yīng)(Coxet al， 2011)。白細胞跨內(nèi)皮遷移在先天性和適應(yīng)性免疫系統(tǒng)應(yīng)答中起重要作用(Muller， 2011)。當機體受損時， 就會發(fā)生炎癥反應(yīng)， 它涉及血液中預(yù)先形成的可溶性元素迅速而短暫地運送到受傷部位， 隨后是較長時間的白細胞運送， 經(jīng)內(nèi)皮細胞遷移可能是炎癥反應(yīng)的不可逆點(Jimenezet al， 2010； Muller， 2011)。補體和凝集通路在機體固有免疫防御中具有重要意義(Bajicet al， 2016)， 在 小 瓜 蟲 感 染 的 大 黃 魚(Larimichthys crocea)早期免疫中發(fā)揮重要作用； 另外， 補體和凝集通路在鰻弧菌感染后半滑舌鰨的免疫反應(yīng)中也發(fā)揮了重要作用(Zhanget al， 2015； Yinet al， 2016)。綜上， 我們推測弓背青鳉的肝臟可能在抵御病原體入侵的免疫防御中發(fā)揮重要作用。

3.3 環(huán)境內(nèi)分泌干擾物的生物標志物

魚類肝臟常被用作檢測環(huán)境中雌激素物質(zhì)的靶器官， 通常用于檢測環(huán)境雌激素物質(zhì)的標記基因在雌魚肝臟中大量表達， 如chgs、vtg、cyp450s和er(Bucheliet al， 1995； Arukweet al， 2001； Leeet al， 2002； Chenet al， 2008)等。Chgs基因編碼放射帶蛋白(ZPs)的前體， ZPs 占據(jù)了魚卵殼的大部分， 由ZP-1、2和ZP-3 兩個主要亞基群組成(Hamazakiet al，1989； Murataet al， 1991)。在性成熟雌魚中， 雌激素誘導(dǎo)Chgs 在肝臟中合成， 隨血液流動整合到放射帶中。當雄魚被環(huán)境雌激素物質(zhì)刺激時， 肝臟也產(chǎn)生Chgs。在日本青鳉中， 雙酚A、壬基酚和乙炔雌二醇均可以誘導(dǎo)Chgl 和Chgh 的表達， 表達量的增加與雌激素劑量呈正相關(guān)(Leeet al， 2002)。Chen 等(2008)研究表明海水青鳉chgl和chgh基因?qū)Νh(huán)境雌激素物質(zhì)也有相似的響應(yīng)。卵黃蛋白原是硬骨魚類卵黃蛋白前體， 雖然Vtgs 在成熟的雌性體內(nèi)自然存在， 而在雄性體內(nèi)則不存在， 但雌性激素或環(huán)境雌激素分泌物可以誘導(dǎo)雄性產(chǎn)生Vtgs。因此， 雄性肝臟中的vtgs 轉(zhuǎn)錄產(chǎn)物， 可用于指示檢測對象過去或當前暴露于雌激素或具有雌激素效應(yīng)的環(huán)境干擾物。本研究在弓背青鳉肝臟轉(zhuǎn)錄產(chǎn)物中發(fā)現(xiàn)了多種環(huán)境雌激素的生物標志物， 這些基因在性成熟雌性肝臟中的表達水平明顯高于雄性肝臟， 可用于將弓背青鳉開發(fā)成廣鹽性環(huán)境監(jiān)測模式物種。

3.4 SSR 標記

SSRs 是一種重要的分子標記， 可作為群體遺傳學研究的資源。本研究， 從18326 個unigenes 中共獲得28195 個SSRs。在已鑒定的SSRs 中， 單核苷酸重復(fù)基序最為豐富， 重復(fù)核苷酸較多的SSRs 較少(表2)， 這與斑石鯛(Oplegnathus punctatus)(Duet al， 2017)、花蟹(Charybdis feriatus)(Zhanget al， 2017)及多鱗鱚(Sillago sihama)(Tianet al， 2019)的SSRs 分析相一致。本研究獲得的弓背青鳉SSR 標記， 為遺傳多樣性研究提供了數(shù)據(jù)， 也有助于進一步繪制弓背青鳉遺傳連鎖圖譜和開發(fā)弓背青鳉遺傳資源。

4 結(jié)論

本研究首次獲得了弓背青鳉雌雄肝臟轉(zhuǎn)錄組。共獲得49912 個unigenes， 48391 個unigenes 被定位到主要數(shù)據(jù)庫， 豐富了弓背青鳉的功能基因資源。通過功能分析推測肝臟可能在 弓背青 鳉抵御病原體入侵的免疫防御中發(fā)揮作用， 同時鑒定到了多個環(huán)境雌激素生物標志基因(vtgs，chgs，er，cyp450等)， 為將其開發(fā)成環(huán)境監(jiān)測物種提供了基礎(chǔ)。此外， 從18326 個含有SSR 序列的unigenes 中獲得了28195 個SSR， 為弓背青鳉遺傳多樣性的研究提供數(shù)據(jù)。