《Frontiers in Microbiology》：Elaborating the molecular characteristics of corals’ different tolerance to environmental stress in Sanya Luhuitou based on multi-omics analysis
1 Introduction
珊瑚礁是地球上最重要的生态系统之一，因其高初级生产力和生物多样性常被视为海洋中的热带雨林。它们为人类提供食物供应、渔业、旅游以及抵御自然灾害的海岸保护等多种生态系统服务。尽管仅占海洋表面的0.2%，珊瑚礁却是至少25%已知海洋生物的家园。然而，海洋热浪导致水温异常升高，严重破坏了海洋生态系统和全球生物多样性。全球范围内的海洋热浪严重影响了珊瑚礁生态系统，导致珊瑚白化和死亡。此外，过度捕捞、污染排放和沿海开发等人类干扰正驱使珊瑚礁迅速衰退。因此，珊瑚适应环境扰动的能力对于珊瑚礁在全球气候变化压力下的生存至关重要。
实验研究和野外调查均发现，不同珊瑚物种对环境胁迫的敏感性存在显著差异。例如，在南中国海三亚岸礁，枝状珊瑚Pocillopora比块状珊瑚Porites更容易发生热白化。珊瑚对环境胁迫的生理响应非常复杂，因为珊瑚是与甲藻和微生物组共生的全生物。珊瑚对环境胁迫的耐受性受多种因素影响，如环境记忆、珊瑚群体形态和大小、珊瑚生长形式以及遗传多样性。
珊瑚全生物通过涉及基因表达、蛋白质激活和代谢重编程的多层次调控网络来应对环境胁迫（如热胁迫、酸化、富营养化）。在热胁迫期间，珊瑚宿主中的热休克蛋白（HSPs）过度表达，与珊瑚免疫反应相关的基因上调。在共生甲藻中，与光合作用、代谢、抗氧化活性和免疫反应相关的基因表达发生变化，触发珊瑚全生物中的热胁迫级联反应。这导致活性氧（ROS）的过量产生，并释放到珊瑚宿主细胞中。ROS的过量产生和毒性积累对珊瑚宿主和共生藻类造成损害。同时，在环境胁迫下，共生甲藻保留更多有机碳供自身利用，导致珊瑚处于“碳限制”状态，珊瑚与共生甲藻之间的共生关系被破坏。此外，珊瑚细菌群落的动态与其环境耐受性相关。然而，其潜在的分子特征尚未完全阐明。通常使用转录组学、蛋白质组学或代谢组学等单组学方法来表征环境胁迫下珊瑚的基因表达谱、蛋白质表达模式和代谢物特征。相比之下，多组学整合分析对于解析珊瑚物种间环境耐受性的种间差异至关重要。鉴于环境胁迫导致的大规模珊瑚白化和死亡，揭示与珊瑚环境耐受性相关的分子特征并理解其生理适应背后的分子机制，对于在未来气候变化情景下建立珊瑚礁模型至关重要。
三亚鹿回头在2010年、2015年和2017年观测到珊瑚热白化。然而，在白化事件后，优势珊瑚物种包括Pocillopora damicornis（PD）、Porites lutea（PL）和Galaxea fascicularis（GF）在正常生长条件下（未经人工胁迫处理）的细菌群落、基因表达、蛋白质和代谢物差异的内在分子特征尚不清楚。本研究假设不同珊瑚物种表现出与环境胁迫相关的独特分子特征。本研究评估了所选优势珊瑚物种对环境胁迫的抵抗力和恢复力，从而为未来三亚湾珊瑚礁的环境耐受性提供参考。
2 Materials and methods
2.1 Sample collection
本研究于2021年6月7日从中国三亚鹿回头岸礁（18°12′19′′N, 109°28′27′′E）采集了珊瑚Pocillopora damicornis（PD）、Porites lutea（PL）和Galaxea fascicularis（GF）。每种珊瑚物种包括来自不同个体（彼此相距约10米）的三个生物学重复。珊瑚采集深度为1.5–2米。海水温度为27–28°C，盐度约为34‰。采集的珊瑚样品立即在液氮中速冻，并转移至-80°C保存直至DNA提取。
2.2 DNA extraction, Illumina MiSeq sequencing and analysis
使用E. Z. N. A.? soil DNA Kit（Omega, United States）分离总基因组DNA。16S rRNA基因（V3–V4）的PCR扩增按照我们先前研究中的详细方法进行。使用等摩尔浓度的纯化PCR产物构建测序文库。通过Gene-Denovo（广州）在Illumina HiSeq 2500平台上进行双末端测序。原始测序读数随后存入NCBI序列读长档案（SRA）。
解复用后，合并序列并通过fastp（v0.19.6）进行质量过滤，去除较长序列（>275 bp）和含有模糊碱基对的序列。为了进一步降低噪音，使用Qiime2（版本2020.2）中的DADA2插件对高质量序列进行降噪处理。然后生成扩增子序列变体（ASVs）。将每个样本的序列数稀释至4000条。使用Na?ve Bayes一致性分类器以及SILVA 16S rRNA数据库（v138）对ASVs进行物种分类。基于每个样本的ASV丰度进行主坐标分析（PCoA）。使用相似性分析（ANOSIM）、置换多元方差分析（Adonis）和多重响应置换程序（MRPP）来检验三种珊瑚物种微生物组的差异。
2.3 RNA extraction, sequencing and transcriptomic data processing
按照先前描述的方法，我们使用Trizol试剂（Invitrogen, Carlsbad, CA, United States）提取总RNA。使用Agilent 2100 Bioanalyzer（Agilent Technologies, Palo Alto, CA, United States）评估分离RNA的质量。使用Hieff NGS Ultima双模式mRNA文库制备试剂盒（上海翊圣生物科技有限公司，中国上海）进行文库制备。制备好的文库随后在Illumina HiSeq? 4000系统上进行测序，由基因迪诺生物技术有限公司（中国广州）提供。
使用默认参数的SeqPrep和Sickle进行接头修剪和原始测序读数的质量过滤。使用Bowtie2（版本2.2.8）将读数映射到核糖体RNA（rRNA）数据库，随后去除rRNA映射的读数。使用HISAT2在链特异性模式下通过参考基因组引导比对将读数与参考基因组进行比对。然后使用StringTie（v2.2.1）应用基于参考的转录本组装方法，并使用每千碱基每百万映射读数的读数（RPKM）来估计基因表达。使用DESeq2进行差异表达分析。我们识别了显著差异表达基因（DEGs），标准为|log2FC| > 1且错误发现率（FDR）< 0.05。
共生甲藻的转录组分析使用de novo组装进行，遵循我们先前的方法。通过HISAT2将组装的unigenes映射到上述索引的共生甲藻参考基因组；最后，计算映射的unigenes并标准化为RPKM。
使用KEGG功能富集分析来识别在代谢通路中显著富集的DEGs（与全转录组背景相比，Bonferroni校正p值≤ 0.05）。使用Goatools进行KEGG通路分析。
2.4 Protein extraction and analysis
使用含有7 M尿素、2 M硫脲和1% SDS的缓冲液从珊瑚样品中提取总蛋白。使用BCA蛋白测定法测定蛋白浓度。将50微克蛋白转移至新的Eppendorf管中，并用1 M DTT 1微升将最终体积调整至50微升，并在55°C下孵育1小时。在室温下将样品在5微升1 M IAA中孵育1小时。加入五倍体积的-20°C预冷丙酮沉淀蛋白质过夜。用预冷的90%丙酮水溶液洗涤两次后，将沉淀重悬于碳酸氢铵（50 mM）中。使用测序级修饰胰蛋白酶（Promega, Madison, WI）以1:50（酶:蛋白，重量:重量）的比例在37°C下将蛋白质消化过夜。
使用反相分离在高pH值下分级肽混合物。分级分离完成后，使用含有0.1%甲酸的水的溶剂A重新溶解肽段。然后使用与EASY-nLC 1200系统（Thermo Fisher Scientific, MA, United States）耦合的Orbitrap Fusion Lumos进行在线纳喷雾LC-MS/MS分析。将肽段样品加载到分析柱（Acclaim PepMap C18, 75 mm × 25 cm）上。使用120分钟的梯度洗脱色谱柱，其中溶剂B的比例从5%增加到35%。溶剂B含有0.1%甲酸的乙腈。分析柱的流速为200 nL/min，电喷雾电离相对于质谱仪入口的电压为2 kV。使用数据非依赖采集（DIA）模式获取质谱数据，允许全面且无偏地分析肽离子。
使用Spectronaut X（Biognosys AG, Switzerland）处理DIA数据，使用默认设置针对从RNA-seq结果获得的珊瑚和共生甲藻蛋白质数据库进行分析。对所有通过过滤器的前体进行定量。通过平均通过1% Q值阈值的前三个过滤肽段来计算主要组数量。使用Student t检验识别差异表达蛋白（DEPs）。使用两倍变化阈值和Q值 < 0.05来识别DEPs。使用KEGG数据库注释所有鉴定的蛋白质，然后使用DEPs分析KEGG通路。
2.5 Metabolite extraction and analysis
对于代谢物提取，将100微升样品与300微升甲醇和20微升内标混合，涡旋约30秒，并在冰水浴中进行超声提取5分钟。然后将混合物在-20°C下储存2小时，随后在4°C下以13,000 rpm离心15分钟。将200微升上清液等分转移到2毫升样品瓶中用于后续LC-MS分析。通过汇集所有实验样品提取物的等体积来制备质量控制（QC）样品，以评估样品处理的重现性。在仪器分析过程中，每6个测试样品插入一个QC样品，以监控整个分析过程的稳定性。
使用UHPLC系统与Q Exactive Orbitrap高分辨率质谱仪联用进行代谢物检测，获取一级和二级质谱数据。针对离子模式优化流动相：正离子模式下，0.1%甲酸水溶液（A）和乙腈（B）；负离子模式下，5 mM乙酸铵水溶液（用氨水调节pH至9.0，A）和乙腈（B）。使用Full MS-ddMS2方法分析QC样品，具有四个分段扫描范围（70–200 m/z, 190–400 m/z, 390–600 m/z, 和590–1,000 m/z），每个扫描一次以增强代谢物鉴定的二级数据采集。质谱分析后获得的原始数据经过基线过滤、峰识别、积分、保留时间校正、峰对齐和归一化处理，产生包含保留时间、质荷比和峰强度的数据矩阵。通过比对HMDB和Metlin等公共数据库对代谢物进行定性鉴定。
使用Bioconductor的ropls R包进行正交偏最小二乘判别分析（OPLS-DA），以确定可比组之间的全局代谢变化。通过累积确定系数（R2）和交叉验证参数（Q2）验证模型有效性。尽管存在夸大组分离的潜在局限性，但使用双阈值筛选策略识别显著代谢物：来自OPLS-DA模型的变量重要性投影（VIP）得分 > 1.0 结合配对t检验p值 < 0.05。识别两组之间不同的代谢物。然后绘制这些代谢物的生化通路。使用KEGG进行代谢富集和通路分析。
2.6 mRNA and protein correlation analyses
对基因和蛋白质之间的关联进行了定量分析。分别分析显示差异表达的基因和蛋白质的转录组和蛋白质组。使用九象限图分析来识别在GF和PL中一致变化的DEGs/DEPs。使用R（版本3.5.1）进行分析。为了进一步分析珊瑚全生物是否在功能或分子类别上也存在差异，我们基于具有显著一致变化的DEGs/DEPs，使用阈值（Q < 0.05）进行了KEGG通路富集分析。
3 Results
3.1 Diversity and community assembly of bacterial community in three coral species
根据稀疏曲线，所有珊瑚样品均显示出足够的覆盖度。三种珊瑚物种相关细菌群落的Alpha多样性不同，Pocillopora damicornis（PD）的Shannon指数显著低于其他两种珊瑚物种（p < 0.05）。然而，Porites lutea（PL）和Galaxea fascicularis（GF）的Shannon指数之间没有显著差异（p > 0.05）。珊瑚相关细菌群落的其他Alpha多样性指数与Shannon指数模式相似。细菌群落的主坐标分析（PCoA）结果表明，三种珊瑚物种的珊瑚相关细菌群落存在显著差异。这些差异得到了三种差异检验的证实，包括ANOSIM、Adonis和MRPP（p < 0.01），它们都证实了PCoA的结果。
物种分类表明，变形菌门（Proteobacteria）是所有珊瑚样品相关细菌群落中的优势群，分别占PD的70.6%、PL的49.9%和GF的57.0%。我们发现，包括变形菌门、浮霉菌门（Planctomycetes）、芽单胞菌门（Gemmatimonadetes）、硝化螺旋菌门（Nitrospirae）、螺旋体门（Spirochaetes）和古菌门（Nanoarchaeaeota）在内的六个门在PD和PL之间的丰度存在显著差异，而11个门在PD和GF之间存在显著差异（p < 0.05）。此外，PD中浮霉菌门、芽单胞菌门和硝化螺旋菌门的丰度显著低于PL和GF。在PL和GF之间，只有两个门（厚壁菌门（Firmicutes）和螺旋体门）显示出显著差异（p < 0.05）。在属水平上，PD中不动杆菌属（Acinetobacter）、Endozoicomonas属、Ralstonia属和Sphingomonas属的丰度在统计学上显著高于PL和GF（p < 0.05）。相比之下，PD中Ruegeria属、Chlorobium属和Pelagibius属的丰度低于PL和GF（p < 0.05）。与PD相比，PL和GF的细菌群落组成更相似，前11个属的丰度在PL和GF中均高于或低于PD，除了Fulvivirga属。
3.2 Differentially expressed genes among coral host and Symbiodiniaceae
三种珊瑚物种和相关共生甲藻的De novo组装转录组分别包含112,655和36,938个推定基因。初步分析显示，PDGF（PD versus GF）是一个异常值，在珊瑚宿主转录组中检测到64,740个DEGs，而在PDPL（PD versus PL）和PLGF（PL versus GF）中分别为82,227和83,555个。与PL相比，PD中共有37,185个基因表现出较高的基础表达，而45,042个基因表现出较低的基础表达。与GF相比，PD也表现出32,838个基因具有显著较高的基础表达，31,902个基因具有显著较低的基础表达。与珊瑚宿主类似，差异基因表达表明PDGF是一个异常值，在共生甲藻转录组中检测到1,535个DEGs，而在PDPL和PLGF中分别为36,503和36,353个。具体而言，大多数DEGs（975或64%）在PD中相对于GF表现出显著较低的基础表达。与PL相比，PD中有18,543个基因表现出显著较高的基础表达，而17,960个基因表现出显著较低的基础表达。当比较PL和GF时，PL中有17,820个基因具有显著较高的基础表达，而18,533个基因具有显著较低的基础表达。
发现24个环境胁迫相关的DEGs在PDGF和PDPL组之间共享。八个基因共享相同的表达模式（即均表现出较高的基础表达）：Maf、基质金属蛋白酶（MMP）、G蛋白偶联受体（GPCR）、肌球蛋白VIIA（MYO7A）、核糖体蛋白L37（RPL37）、肿瘤坏死因子受体相关因子3（TRAF3）、过氧化物酶（PXDN）和延伸因子1α（EF1α）。此外，我们发现PL和GF之间热胁迫相关基因的表达模式相似。
3.3 Differences in the proteomes among three coral host and Symbiodiniaceae
蛋白质组学分析显示，大多数表达的蛋白质在10至25 kDa范围内，分别为珊瑚宿主和共生甲藻鉴定出16,485和626个蛋白质。基于Q值 < 0.05时≥2倍变化的选择标准识别差异表达蛋白（DEPs）。然后，在珊瑚宿主的PDPL、PDGF和PLGF组中分别鉴定出492、574和449个蛋白质作为DEPs。PDPL组中绝大多数DEPs（331或67%）和PDGF组中绝大多数DEPs（326或57%）在PD中表现出较低的基础表达。特别是，与GF相比，PL中显示出较低基础表达的DEPs较少（207或46%）。在共生甲藻的蛋白质组中，在PDPL、PDGF和PLGF组中分别发现总共30、31和17个蛋白质是DEPs。与GF和PL相比，更多的DEPs在PD中表现出较高的基础表达（在PDPL和PDGF中分别为16和18个）。在PLGF组中，PL中有11个蛋白质相对于GF具有较高的基础表达。PDGF和PDPL组共享五个热胁迫相关的DEPs。其中，三个DEPs具有相同的表达谱：核糖体蛋白S9、钙调蛋白（CaM）和FK506结合蛋白12（FKBP12）。
3.4 Metabolic profiling analysis of three coral holobionts
基于三种珊瑚全生物的代谢谱分析，我们发现PDPL、PDGF和PLGF组中分别有778、764和607个代谢物显著差异表达。超过一半的差异表达代谢物在PD中相对于PL和GF表现出较高的基础表达（分别为457和401个）。然而，更多的差异表达代谢物在PL中相对于GF表现出较高的基础表达（385或63%）。包括硫酸胆固醇和睾内酯在内的几种代谢物在PD中富集。然而，与PD相比，脱氧胞苷和L-谷氨酰胺在GF中显著更高，在PL中也发现了脱氧胞苷的相似模式。尽管不显著，LysoPC在PD中的丰度高于GF和PL。
3.5 Correlation between the transcriptome and the proteome
我们分析了珊瑚宿主和共生甲藻的基因和蛋白质表达之间的相关性，结果显示所有珊瑚宿主之间存在类似弱但高度显著的相关性（p值 < 0.01）。在PDPL和PLGF组中，共生甲藻的基因和蛋白质之间存在弱且不显著的相关性。因此，珊瑚物种不直接影响基因和蛋白质之间的相关性。此外，我们评估了珊瑚宿主和共生甲藻的DEGs和DEPs之间的相关性。有趣的发现表明，先前识别的弱相关性在珊瑚宿主和共生甲藻中均得到增强。我们假设在珊瑚物种之间无显著差异的蛋白质可能驱动了弱相关性。很可能DEPs会对应在相同方向上差异富集的基因。因此，去除非显著表达的蛋白质导致了更高的相关性。
此外，我们通过九象限关联分析研究了基因和蛋白质的关联，发现对于珊瑚宿主的PLGF组，更多的基因（11或73%）在GF中表现出较低的基础表达。相反，对于PDGF（11或16%）和PDPL（18或16%）组，只有少量基因（或蛋白质）在PD中表现出较高的基础表达。在共生甲藻样品中，很少有基因（或蛋白质）具有一致的显著变化。对于PDGF组，只有1个基因在PD中表现出显著较低的基础表达，而对于PDPL组，有1个基因在PD中表现出显著较高的基础表达。
3.6 Function investigation based on KEGG enrichment analysis
为了评估功能富集，对具有一致变化的DEGs/DEPs进行了KEGG通路分析。对于珊瑚宿主，与PL相比，GF中有8条通路显示基础水平升高，33条通路基础水平降低。与PD相比，GF有110条通路基础水平升高，26条基础水平降低，而PL有239条通路基础水平升高，20条基础水平降低。与PD相比，与肌动蛋白细胞骨架调节相关的通路在GF和PL中富集程度最高。此外，与溶酶体相关的通路在GF中相对于PD也富集，而名为紧密连接、致病性大肠杆菌感染、Rap1信号通路、沙门氏菌感染和志贺氏菌病的通路同时在PL中相对于PD富集。特别是，与精氨酸和脯氨酸代谢相关的通路在PL和GF中相对于PD均富集，并且一些与蛋氨酸、色氨酸、丙氨酸、丝氨酸、缬氨酸和异亮氨酸代谢相关的通路也在PL中相对于PD富集。
4 Discussion
珊瑚对环境胁迫的敏感性与环境记忆、生长形式和珊瑚群体形态有关。具体来说，预先经历胁迫的珊瑚表现出更高的环境胁迫抵抗力。具有较厚组织层的块状珊瑚通常对环境胁迫具有更高的耐受性。相比之下，枝状珊瑚通常表现出对环境胁迫更高的敏感性。在此，我们使用多组学方法结合微生物群落来探索PD、PL和GF之间不同环境胁迫敏感性的分子特征。
4.1 The microbial patterns across PD, PL, and GF
细菌微生物组成员通过一系列过程促进珊瑚健康，包括硫循环、固氮和抵御病原体，珊瑚相关微生物群落与其对环境胁迫因子的敏感性和抵抗力相关。珊瑚礁生态系统正在经历日益衰退，这凸显了研究珊瑚宿主与其相关细菌群落之间关系的重要性。在本研究中，我们发现与PD相关的细菌群落的Alpha多样性显著低于PL和GF。多样化和健康的微生物组对于宿主在不断变化的环境条件下生存至关重要。虽然先前的研究表明PD在白化期间表现出增加的微生物多样性，但健康PD中细菌群落相对较低的多样性可能导致其对环境胁迫的脆弱性更大。进一步的研究表明，珊瑚微生物组的多样性与疾病易感性无关。相反，微生物组中少数优势微生物类群是珊瑚疾病易感性的预测指标。不动杆菌属（Acinetobacter）被认为既是有效的第一道防线，也是潜在的病原体；此外，高温使其能够在珊瑚微生物组中承担核心作用。Ralstonia被发现共定位于珊瑚宿主的内共生甲藻和胃层细胞内，它也被发现是一种机会性病原体，其丰度在弧菌感染后显著增加。Sphingomonas是珊瑚礁内共生甲藻细胞中的细胞内核心细菌。此外，该属中的一些物种也是潜在的珊瑚病原体，并且已发现未鉴定的Sphingomonas物种可能参与“白瘟疫”，导致珊瑚组织以异常速率降解。值得注意的是，不动杆菌属、Ralstonia属和Sphingomonas属的丰度在PD中均高于PL和GF。先前的研究表明，几种Ruegeria菌株对V
