经济贝类是重要的食品来源，也是海洋生态系统中不可或缺的部分。我国科研团队已经对一些重要的经济贝类开展了大量的遗传研究，包括基因组、转录组、表观遗传和分子标记等，如栉孔扇贝(Chlamys farreri)^[1]、近江牡蛎(Crassostrea ariakensis)^[2]、硬壳蛤(Mercenaria mercenaria)^[3]、缢蛏(Sinonovacula constricta)^[4]等全基因组精细参考图谱的陆续发布。贝类转录组测序的研究为功能基因、分子标记开发提供了数据支撑，加速了贝类遗传育种研究的进程：通过转录组测序得到了温度胁迫下马氏珠母贝(Pinctada fucata martensii)和缢蛏的差异表达基因，为双壳贝类对高低温胁迫的响应机制提供了理论基础^{[5, 6]}；通过转录组测序获得了泥蚶(Tegillarca granosa)耐低氧的候选基因^[7]；利用转录组学分析可以研究贝类在不同毒性胁迫下的免疫应答机制，了解其机体抗氧化酶系统及解毒代谢机制^{[8, 9]}。对贝类不同组织和不同发育阶段进行转录组解读，可了解其器官分化的相关基因^{[10, 11]}和变态发育过程中的关键基因^[12]。

弓獭蛤(Lutraria rhynchaena Jonas, 1844，同物异名Lutraria arcuata Reeve, 1854)隶属于软体动物门(Mollusca)瓣鳃纲(Lamellibranchia)帘蛤目(Veneroida)蛤蜊科(Mactridae)獭蛤属(Lutraria)^{[13, 14]}。弓獭蛤广泛分布在我国、日本、菲律宾和越南等沿海国家，在我国主要分布在广东、海南、广西和台湾沿海^{[13, 14]}，栖息于潮间带至水深80 m的浅海泥沙质海底，与獭蛤属的大獭蛤(L.maxima)、施氏獭蛤(L.sieboldii)和澳洲獭蛤(L.australis)等形态相似，在我国沿海居民皆俗称它们为象鼻螺。象鼻螺个体较大，每只成年个体质量90-120 g，我国自1997年开始通过浅海开放式底播方式人工养殖象鼻螺，养殖范围主要集中在南方，2020年北海的象鼻螺产量达6万多吨^[15]。2022年，北海象鼻螺获得了中华人民共和国农产品地理标志登记证书(中华人民共和国农业农村部公告第532号)。

因弓獭蛤具有可观的经济价值，国内外学者对其开展了多方面的研究，国外关于弓獭蛤的研究主要由越南学者发表。越南学者发布了基于第三代纳米孔和第二代Illumina NovaSeq6000平台测序组装的弓獭蛤基因组草图，基因组大小为545-547 Mbp^[16]，但此后未有基于该基因组的拓展研究发表。研究者对越南弓獭蛤的地理分布、繁育技术、线粒体基因组、微卫星等也做了一些报道^[17-19]。目前国内关于弓獭蛤研究的报道主要局限于苗种繁育技术^[20]、受精卵的细胞学观察^[21]、重金属响应^[22]等。有研究者对弓獭蛤使用了高浓度的镉胁迫，发现其不同组织对重金属的胁迫响应存在显著差异，但是并未进行分子生物学层面的深入研究。总体来说，国内目前对弓獭蛤的遗传育种、功能基因表达等方面研究还处于起步阶段。

本研究以弓獭蛤的肌肉(斧足)、外套膜、鳃、肝胰腺、虹吸管、雌性性腺及雄性性腺7种组织为研究对象，采用RNA-seq测序技术进行转录组测序，通过数据组装、各组织表达量的比较和对差异表达基因注释获得弓獭蛤各组织的基因表达信息，为进一步开展弓獭蛤遗传育种、分子标记、功能基因的挖掘等方面研究提供数据参考。

1 材料与方法 1.1 材料

测序用的弓獭蛤贝(其贝壳见图 1)采集于北海市金不换水产有限公司位于广西北海市侨港镇(109°9′30.85″N，21°18′2.92″E)的天然海域养殖场，均为成年活贝，平均湿重为(109.4±30.6) g、壳长为(92.5±17.8) mm、壳高为(47.1±3.8) mm、壳宽为(31.9±3.9) mm。取斧足、外套膜、鳃、肝胰腺、虹吸管、雌性性腺及雄性性腺7种组织，每种组织取9个个体样本，设置3个生物学重复，每3种组织样本混合后提取RNA。实时荧光定量PCR(Real Time qPCR)验证实验用的弓獭蛤贝购自北海市南珠市场。

1.2 方法 1.2.1 RNA的提取及检测

使用MiniBEST通用RNA提取试剂盒(TAKARA公司)提取各组织总RNA，提取步骤按照试剂盒说明书进行。使用NanoDrop分光光度计(Thermo Fisher Scientific公司)检测总RNA纯度，使用Qubit 2.0 (Thermo Fisher Scientific公司)检测总RNA浓度，琼脂糖凝胶电泳检测总RNA完整性。

1.2.2 文库构建及转录组测序

检测合格的总RNA送上海凌恩生物科技有限公司，通过带有Oligo (dT)的磁珠富集具有polyA尾巴的真核mRNA后，采用Illumina Truseq RNA sample prep Kit方法完成7种组织共21个样品的文库构建，经质检合格后开展基于Illumina HiSeq 2500的RNA-seq测序。

1.2.3 生物信息学分析

使用FastQC对测序得到的原始数据进行质量评估，使用trimmomatic软件(v0.36)去除reads接头序列、将序列末端(3′端)质量较低(质量值小于20)的碱基修剪掉，去除含N比率超过10%的reads，舍弃长度小于75 bp的序列，过滤后得到高质量的测序数据(Clean reads)。利用Trinity软件(v2.8.6)进行无参考转录组从头组装得到非冗余基因(Unigene)并预测开放阅读框(ORF)。用Salmon(v1.1.0)软件将质控后得到的高质量序列和拼接结果进行比对(Mapping)。将组装的序列分别在NR(Non-Redundant protein sequence database)、GO(Gene Ontology)、eggNOG(evolutionary genealogy of genes: Non-supervised Orthologous Groups)、KEGG(Kyoto Encyclopedia of Genes and Genomes)及Swiss-Prot等数据库中进行基因功能注释。

使用RESM v1.3.3软件计算各转录本在样本中的表达丰度，根据RPKM(Reads Per Kilobase per Million mapped reads)判断基因表达量，以P＜0.05和|log₂(FoldChange)|≥1为标准筛选差异表达基因(Differentially Expressed Genes，DEGs)。通过对不同组织之间DEGs进行GO和KEGG通路富集分析，确定DEGs参与的生化过程和信号转导途径。

1.2.4 实时荧光定量PCR验证分析

随机选取9个差异基因进行实时荧光定量PCR实验，以判断RNA-seq数据结果是否可信。经前期实验筛选内参后，以β-actin作为内参基因。利用Primer premier 5.0软件设计引物如表 1所示。使用PrimeScrip^TM RT reagent Kit with gDNA Eraser试剂盒(TAKARA公司)将总RNA反转录合成cDNA。以反转录产物为模板，使用TB Green Premix Ex Taq Ⅱ(Tli RNaseH Plus)试剂盒(TAKARA公司)在QuantStudio^TM3实时荧光定量PCR仪(Thermo Fisher Scientific公司)上进行基因表达相对定量分析。PCR反应体系为20 μL，PCR程序参照产品说明书进行。样品进行3次生物学重复和3次技术重复，用2^-ΔΔCt方法计算基因的相对表达量。

2 结果与分析 2.1 转录组测序结果 2.1.1 测序数据及转录本拼接

测序原始数据及质控数据统计结果如表 2所示。原始数据、有效数据、有效碱基为3个生物学重复样品之和，质控后Q20和Q30碱基比例为3个生物学重复样品的平均值。共测得1 039 382 252条raw reads，质控过滤后得到939 197 518条clean reads，各组织有效碱基数从15.39至29.26 G不等。各组织Phred数值分别大于20、30的碱基占总体碱基的百分比均分别大于等于96.93%和91.96%，数据质量达到后续分析条件。

2.1.2 基因功能注释

拼接得到56 773条Unigene，N50为1 373，平均长度为1 033.12 bp。将预测基因的蛋白序列在NR、GO、eggNOG、KEGG和SWISS-Prot等数据库中比对，分别得到24 557、13 094、17 524、10 352和13 857条注释信息(各数据库注释韦恩图见图 2)。根据GO注释将基因功能分为3大类别：生物过程(BP)、细胞成分(CC)和分子功能(MF)，细分为45个分支(图 3)。通过eggNOG数据库进行基因对应的COG(Clusters of Orthologous Groups of proteins)注释，按注释结果进行功能分类得到了25个功能类别(图 4)，其中注释基因数量最多的功能为信号转导机制，其次是功能未知的基因，然后是翻译后修饰、蛋白转运及伴侣蛋白。通过KEGG数据库得到10 352条注释信息，涵盖了5个功能大类：新陈代谢、遗传信息处理、环境信息处理、细胞过程和有机系统的生物过程通路(图 5)。通过SWISS-Prot数据库得到13 857条注释信息，注释结果可用作基因功能描述的参考。经5个数据库比对共有24 593条Unigene获得注释，注释率为43.32%。

2.2 弓獭蛤不同组织间的差异分析 2.2.1 不同组织基因表达分析

RPKM能够消除基因长度和测序数据分析中的基因表达水平估算产生的影响。基于RPKM估算基因转录本丰度，7种组织中基因表达数量由多到少依次为肝胰腺＞鳃＞外套膜＞雄性性腺＞肌肉＞虹吸管＞雌性性腺(图 6)。在此基础之上，对雌性性腺与雄性性腺、鳃与肝胰腺及各组织与肌肉组织之间的表达基因进行差异比对，肝胰腺vs肌肉组织获得差异基因数量最多，有11 074条，其次是鳃vs肌肉组织有8 938条，鳃vs肝胰腺有6 896条DEGs，虹吸管vs肌肉组织有5 174条DEGs，雄性性腺vs肌肉组织有5 128条DEGs，外套膜vs肌肉组织有2 398条DEGs，雌性性腺vs雄性性腺有703条DEGs，雌性性腺vs肌肉组织有658条DEGs，其各自上调与下调的差异基因数量如图 7所示。

2.2.2 组织间差异基因的GO富集分析

对差异基因进行GO注释，各组织间DEGs的GO条目富集在生物过程、细胞成分和分子功能3大类上的数量分布情况如表 3所示。

雌性性腺vs雄性性腺的DEGs富集程度较显著的GO条目主要包括纤毛移动(Cilium movement)、轴动蛋白复合体组装(Axonemal dynein complex assembly)、基因丝装配(Axoneme assembly)、9+2运动纤毛(9+2 motile cilium)和基因丝(Axoneme)。鳃vs肝胰腺的DEGs富集程度较显著的GO条目主要为甘露糖结合(Mannose binding)、磺基转移酶活性(Sulfotransferase activity)、角质细胞增殖的调控(Regulation of keratinocyte proliferation)和碳水化合物结合(Carbohydrate binding)。雌性性腺vs肌肉组织的DEGs富集程度较显著的GO条目主要为内质网靶向蛋白(Protein targeting to ER)、内质网蛋白定位的建立(Establishment of protein localization to endoplasmic retic)及游离胞浆核糖体(Cytosolic ribosome)。鳃vs肌肉组织的DEGs富集程度显著的GO条目为共质体(Symplast)和盐胁迫响应(Response to salt stress)。雄性性腺vs肌肉组织DEGs富集的GO条目主要有多聚核糖体(Polysomal ribosome)、SRP依赖协同翻译蛋白靶向膜(SRP-dependent cotranslational protein targeting to membrane)、色素体(Plastid)及协同翻译蛋白靶向膜(Cotranslational protein targeting to membrane)。外套膜vs肌肉组织的DEGs富集程度显著的GO条目主要有凋亡细胞的识别(Recognition of apoptotic cell)、伪足膜(Pseudopodium membrane)和吞噬作用、识别(Phagocytosis, recognition)。肝胰腺vs肌肉组织的DEGs富集程度较显著的GO条目包括共质体、色素体和外封结构(External encapsulating structure)。虹吸管vs肌肉组织的DEGs富集程度较显著的GO条目有质体间质(Plastid stroma)和多聚核糖体。

2.2.3 组织间差异基因的KEGG富集分析

对弓獭蛤各组织转录组的DEGs经过KEGG通路注释分析显示(图 8)，雌性性腺vs肌肉组织差异基因所注释的显著富集(P＜0.05)的通路仅6条，富集的DEGs最多的通路是核糖体(Ribosome)，其次是钙信号通路(Calcium signaling pathway)。此外，在心肌收缩(Cardiac muscle contraction)、间隙连接(Gap junction)、蛋白质的消化和吸收(Protein digestion and absorption)、氧化磷酸化通路(Oxidative phosphorylation)均有差异基因显著富集。雄性性腺vs肌肉组织差异基因所注释的显著富集(P＜0.05)的通路有30条，其中富集DEGs较多的通路主要为次生代谢物的生物合成(Biosynthesis of secondary metabolites)、核糖体(Ribosome)、不同环境下微生物的代谢(Microbial metabolism in diverse environments)及马达蛋白(Motor proteins)等通路，核糖体、马达蛋白、氧化磷酸化和碳代谢等通路的富集程度较高。雌性性腺vs雄性性腺富集通路中，显著富集(P＜0.05)的通路只有8条，富集DEGs数量最多的通路是生热作用(Thermogenesis)，其次是与能量代谢相关的氧化磷酸化通路。其中逆行内源性大麻素信号通路(Retrograde endocannabinoid signaling)仅在雌雄性腺间的差异基因中存在显著富集，因此该通路可能参与了弓獭蛤性腺的分化过程。

作为参与运动的肌肉类组织，外套膜、虹吸管与肌肉组织间的差异基因主要在生物活动、生物合成及代谢相关的通路存在显著富集。外套膜vs肌肉组织差异基因所注释的显著富集(P＜0.05)的通路有26条，各通路富集的DEGs数量较均匀，其中富集程度最强的有光传导(Phototransduction)、粘蛋白O型-聚糖生物合成(Mucin type O-glycan biosynthesis)、间隙连接、嗅觉转导(Olfactory transduction)及唾液分泌(Salivary secretion)等通路。虹吸管vs肌肉组织差异基因所注释的显著富集(P＜0.05)的通路有33条，其中富集DEGs最多的主要有代谢通路(Metabolic pathways)、次生代谢物的生物合成和不同环境下微生物的代谢。富集程度最高的通路主要为碳代谢(Carbon metabolism)、氨基酸的生物合成(Biosynthesis of amino acids)及核糖体。

鳃vs肌肉组织差异基因所注释的显著富集(P＜0.05)的通路有40条，其中富集DEGs最多的通路主要有次生代谢物的生物合成、不同环境下微生物的代谢、碳代谢等通路。肝胰腺vs肌肉组织差异基因所注释的显著富集(P＜0.05)的通路有35条，其中富集DEGs最多且富集程度强的通路主要为代谢通路、次生代谢物的生物合成、不同环境下微生物的代谢和吞噬体(Phagosome)。鳃vs肝胰腺差异基因所注释的显著富集(P＜0.05)的通路有63条，其中富集DEGs最多的通路主要是代谢通路，其次是吞噬体，在唾液分泌、吞噬体及光传导等通路的富集程度较其他通路更显著。根据差异基因富集的通路分析，弓獭蛤的鳃和肝胰腺主要负责机体的消化吸收、免疫应答等重要生理活动。

2.2.4 实时荧光定量PCR验证

为验证所送组织转录组测序结果的可靠性，从获得的全部DEGs中随机挑选了9个进行实时荧光定量PCR验证，通过2^-ΔΔCt计算出相对表达量，结果显示TRINITY_DN32671_c2_g1仅在鳃中具有极显著表达，TRINITY_DN24180_c2_g1、TRINITY_DN39293_c3_g10、TRINITY_DN41770_c1_g1等基因在弓獭蛤肌肉组织中具有极显著表达(图 9)，与测序结果(表 4)基本趋于一致，表明弓獭蛤各组织转录组测序结果可信。

3 讨论 3.1 转录组技术在双壳贝类中的应用

多组学技术和分子生物学技术已经广泛应用于水产品的育种、繁殖、病害鉴定、杂交种鉴定等方面。在我国贝类养殖中具有传统优势的物种，如牡蛎、扇贝等，均已构建有丰富的全基因组、转录组、蛋白组、小RNA组及表观遗传组资源库。对经济贝类进行多组学资源库的构建和研究对推动养殖业的优良种选育和产业的可持续发展有着重要的作用。

转录组可揭示特定组织在特定状态下的基因表达情况，多应用于组织间差异、环境胁迫或毒理条件下的差异表达和功能基因筛选。有研究者对虾夷扇贝(Mizuhopecten yessoensis)的外套膜进行转录组测序，获得16 816个Unigenes^[23]；对海湾扇贝(Argopecten irradians)发育阶段个体进行转录测序，获得70 929个Unigenes^[24]；从牛角江珧(Atrina pectinata)的外套膜、鳃、性腺、后闭壳肌和肝胰腺5个组织中提取RNA，转录测序获得127 263个Unigenes^[25]。

弓獭蛤作为我国南方重要的养殖经济贝类，是北海市的地理标志农产品，其多组学研究还比较匮乏。本研究利用RNA-seq技术对成体弓獭蛤的7种组织进行转录组测序，旨在丰富弓獭蛤的基因数据库资源。虽然NCBI上已有越南学者上传的基于肌肉组织测序的弓獭蛤基因组，但因獭蛤属的几个种形态学极为相似较难区分，因此本文对测序数据采用无参考转录组从头拼接的方法进行组装。对7种不同组织的转录组进行测序分析，比研究单一组织或较少组织在胁迫条件下的基因变化获得的数据更加全面和可信，也为后期进一步开展弓獭蛤遗传、生理、养殖等研究建立数据库资源。

本研究通过RNA-seq技术共获得939 197 518条clean reads。组装及ORF查找后获得了56 773个Unigenes，N50为1 373，平均长度1 033.12 bp，将其放入5个公共数据库中进行对比，有24 593条基因得到注释信息。通过各组织间差异基因的比较发现，肝胰腺vs肌肉组织和鳃vs肌肉组织的差异基因数量较多。

3.2 KEGG通路

雌性性腺vs肌肉组织富集的胰岛素信号通路(Insulin signaling pathway)和卵母细胞成熟(Oocyte meiosis)均与弓獭蛤性腺的发育相关，这两条通路在纵肋织纹螺(Nassarius variciferus)^[26]和钝缀锦蛤(Tapes conspersus) ^[27]中也被提出与性腺发育有关。鳃vs肝胰腺和鳃vs肌肉组织之间在唾液分泌、胃酸分泌通路、碳水化合物、蛋白质的消化吸收等通路有显著富集，这表明鳃是弓獭蛤的消化器官之一。肝胰腺vs鳃及肝胰腺vs肌肉组织间存在吞噬体、病毒体-黄病毒、TRP通道的炎症介质调节和药物代谢-其他酶等通路的显著富集，可证明肝胰腺是负责抵抗疾病、免疫和解毒的关键组织，且经对比在肝胰腺中存在甘氨酸、丝氨酸、苏氨酸、丙氨酸、天冬氨酸和谷氨酸等多种必需氨基酸的代谢通路。外套膜与肌肉组织存在光传导通路的富集，因此推测外套膜会对外界光照强度刺激产生应答，此外，贝类外套膜中表达的基因参与色素沉着、金属转录和免疫^[28]，在本研究中，参与黑色素生成的关键限速酶的酪氨酸酶(Tyrosinase)^{[29, 30]}在外套膜中显著高表达。

3.3 功能注释

在不同物种的同一组织、同一物种的不同组织之间的转录组表达往往存在相似之处也存在差异，本次对所获得的基因和两两组织间的DEGs进行GO功能注释，结果表明均在细胞过程(Cellular process)、生物调节(Biological regulation)、细胞结构体(Cellular anatomical entity)、结合(Binding)等功能上有较多基因富集，这与泥蚶^[7]、钝缀锦蛤^[27]和缢蛏^[31]等贝类具有相似性。本研究弓獭蛤注释率为43.32%，与大部分双壳贝类注释率相当，但与其他研究较深的常见经济动物相比，注释率较低，这表明在贝类中仍有许多基因或序列的功能等待发掘。

双壳贝类对重金属镉暴露有显著的物种、时间和剂量效应的差异^[32]，本研究中，参与贝类重金属镉胁迫响应的抗氧化酶系的相关基因如谷胱甘肽过氧化物酶(Glutathione peroxidase)、硫氧还蛋白(Thioredoxin)、超氧化物歧化酶(Superoxide dismutase)、金属硫蛋白(Metallothionein)、谷胱甘肽硫转移酶(Glutathiones S-transferase)、过氧化氢酶(Catalase)、谷氧还蛋白(Glutaredoxin)等基因在鳃、肝胰腺、外套膜和虹吸管有较高表达量。贝类免疫主要以体液免疫和血细胞免疫的方式进行^[33]，在本研究的结果中，参与体液免疫的溶酶体(Biogenesis of lysosome-related organelles complex)、toll样受体(toll-like receptor)、髓样分化因子(Myeloid differentiation primary response 88)、同种异体移植炎症因子(Allograft inflammatory factor)等相关基因以及凝集素(Lectin)等多集中于鳃、肝胰腺和外套膜表达，而血细胞免疫相关的基因如巨噬细胞相关基因和受体等多集中分布于鳃和外套膜中。

通过对转录组数据的分析，本研究在弓獭蛤的雌雄性腺、肌肉组织、鳃和肝胰腺中分别得到了大量与弓獭蛤性状、生长发育、生理代谢及免疫调节相关的通路及功能基因，但这些通路和基因在弓獭蛤各组织内具体的表达和调控机制仍待进一步研究。

4 结论

多组学研究在经济贝类的遗传育种等方面发挥着重要作用。本研究基于第二代高通量RNA-seq测序技术获得我国南方一种重要的养殖贝类弓獭蛤的转录组序列数据﹐通过组装拼接得到了56 773条预测基因。利用BLAST软件比对NR、GO、eggNOG、KEGG和SWISS-Prot数据库﹐进行了功能基因注释﹐共有24 593条基因得到注释。通过组织间的两两对比获得了雌性性腺vs雄性性腺(703条)、鳃vs肝胰腺(6 896条)、雌性性腺vs肌肉(658条)、鳃vs肌肉(8 938条)、雄性性腺vs肌肉(5 128条)、外套膜vs肌肉(2 398条)、肝胰腺vs肌肉(11 074条)和虹吸管vs肌肉组织(5 174条)的差异表达基因。对各组之间DEGs进行GO注释和KEGG富集分析，讨论了各组之间可能存在的功能差异。本研究为开展弓獭蛤分子遗传学研究提供了基因资源，也为进一步开展弓獭蛤功能基因挖掘、遗传育种研究和分子标记的开发提供了基础参照数据。