宿主转录组揭示共感染期间防御-繁殖权衡的减轻及其机制

本文通过宿主转录组学分析，揭示了在寄生虫共感染过程中，由于寄生虫间的拮抗作用降低疾病严重性，宿主得以减轻防御与繁殖之间的代价权衡。研究利用线虫-细菌模型，结合基因组学与基因共表达网络（WGCNA）分析，发现保护性共感染可显著减弱宿主免疫投资而不改变核心应答基因，同时通过细菌基因组注释提示其保护机制可能与公共资源（如铁载体介导的铁摄取）竞争相关。该研究为理解共生微生物调控宿主免疫trade-off提供了新视角。

摘要

本研究利用转录组学技术，探究了在线虫宿主^{Caenorhabditis elegans}中，单一感染与共感染条件下宿主基因表达的变化，以及这些变化如何影响宿主在防御和繁殖之间的资源分配权衡。通过比较强毒寄生虫^{Leucobacter musarum}（LM）和保护性寄生虫^{Leucobacter celer}（LC）的感染效应，研究发现共感染（LM与LC同时存在）可显著降低宿主的疾病严重程度，并减轻宿主的免疫-繁殖权衡。

1 引言

感染可驱动宿主采取多种防御策略以最小化对健康和适应性的代价。宿主可能通过免疫、行为或共生介导的防御来避免、抵抗或耐受寄生虫。然而，防御是有成本的，跨物种的证据表明，在分配有限资源到防御、繁殖和体质之间存在张力。例如，昆虫的交配努力可能导致免疫功能和寿命减少，而植物的防御化合物生产可能以生长为代价。这些适应度权衡在感染毒力减弱时可能减轻。几乎所有生物都被多种生物定植，涵盖从共生到寄生的整个 continuum。共感染下的 symbionts 可能通过三种主要模式相互竞争：剥削性竞争（如争夺资源）、干扰竞争（如分子战争）或表观竞争（如修饰或诱导宿主免疫）。竞争性共感染的结果可能导致“保护性共感染”，即毒力降低。例如，保护性共感染可通过微生物竞争公共 goods（如铁）、产生抗菌超氧化物或免疫 priming 其宿主来实现。微生物共生体的保护可能导致宿主可塑性减少对驱动强毒感染的寄生虫的特定防御投资，并最终减少对防御机制的选择。本研究旨在通过测量宿主在动物-微生物寄生虫系统中的转录反应，来研究这些权衡。

2 方法

2.1 宿主-微生物系统与感染测定

研究使用^{C. elegans}线虫作为宿主，在涂布^{Escherichia coli}（OP50）作为食物源的Nematode Growth Medium（NGM）平板上维持。强毒由^{L. musarum}驱动，保护由^{L. celer}赋予。感染方法大致遵循Bates等人（2021）的方法。通过漂洗分离线虫卵，在M9缓冲液中孵化L1幼虫，然后在^{E. coli}平板上生长至L4期。所有三种细菌从冷冻甘油 stock 在摇动lysogeny broth（LB）中 overnight 培养，在30°C（^{E. coli}）或25°C（^Leucobacter）培养至光密度（600 nm）约1.0（^{E. coli}）或0.3（^Leucobacter）。感染暴露平板通过在9 cm NGM平板上接种80%^{E. coli}培养物和20%^Leucobacter培养物（总体积200 μL）来制备细菌 lawn。四种处理类型：对照（80%^{E. coli}，20%磷酸盐缓冲盐水）、'LC'（80%^{E. coli}，20%^{L. celer}）、'LM'（80%^{E. coli}，20%^{L. musarum}）和共感染（80%^{E. coli}，10%^{L. celer}，10%^{L. musarum}）。所有平板（生物学重复）在25°C孵化24小时。然后将1500只L4线虫转移到每个处理的八个重复平板，并在25°C孵化所有平板。

在早期感染（10小时）和晚期感染（20小时）两个时间点收集实验样品。选择这些时间是因为在这些^{L. musarum}感染方法中，24小时内观察到显著的^{C. elegans}死亡率（Bates等人，2021）。每个采样时间，从每个处理的四个重复平板收集线虫，用5 mL M9缓冲液 flood 平板。样品通过离心用M9缓冲液洗涤四次。将洗涤后的线虫悬浮在800 μL冰 cold RNA Shield（Zymo Research）中，并立即通过 bead beating 在2800 rpm（Disruptor Genie，Scientific Industries）下进行5个循环的2分钟 disruption—2分钟冰 incubation 来 homogenize。然后将样品在−80°C冷冻直至RNA提取。

使用Quick-RNA Minprep Plus Kit（Zymo Research）按照制造商协议从样品中提取总RNA。总RNA送至CGR进行RNAseq文库制备和Illumina测序。在CGR，使用RiboZero Plus kit（Illumina）去除rRNA以富集mRNA。使用Ultra II Directional kit（NEB）制备 paired-end、indexed、strand-specific RNA文库。文库在NovaSeq（Illumina）上测序，产生2×150 bp reads。每个样品产生约50–70 million reads。CGR使用Cutadapt（v. 1.2.1）进行初步read trimming，使用3 bp匹配阈值开始trimming（Martin，2011）。CGR还使用Sickle（v. 1.200）进行质量trimming，使用Q20窗口和最小read长度15 bp（Joshi和Fass，2011）。

2.2 RNAseq分析以识别感兴趣基因

将reads（仅paired）分配给^{C. elegans}、^{L. celer}、^{L. musarum}或^{E. coli}进行量化。宿主RNA参考转录本可从NCBI under accession GCF 000002985.6获取。使用GffRead（v. 0.12.7）（Pertea和Pertea，2020）为所有三种细菌物种生成转录本序列，使用本文 presented 的^Leucobacter组装和^{E. coli} assembly GCF 009496595.1。使用Salmon（v 1.4.0），将reads pseudo-map 到 concatenated all-organisms transcriptome 并量化RNA转录本丰度（Patro等人，2017）。从这些 mapping 中，然后分离宿主和三种细菌中每一种的转录本丰度数据。每个样品约有20–35 million read pairs map 到宿主转录组，并使用tximport（v. 1.30.0）在RStudio（v. 2023.09.1）与R（v. 4.3.2）（R Core Team，2022；RStudio Team，2020；Soneson等人，2016）中将转录本与基因名称关联。对于LC和共感染处理，每个样品发现7000–350,000^{L. celer} reads。对于LM和共感染处理，每个样品发现34,000–336,000^{L. musarum} reads。跨所有处理，每个样品发现3000–900,000^{E. coli} reads。认为两种^Leucobacter物种的测序深度太低，无法进行稳健分析，因此仅将reads map 到细菌转录组以防止非特异性map 到宿主转录组。丢弃来自非mRNA转录源（如残留rRNA或非编码piRNA）的计数。

使用DeSeq2（v. 1.42.0）标准化计数数据并计算实验处理之间的 pairwise 差异基因表达（Love等人，2014）。从Wald测试中选择上调和下调的差异表达基因（DEGs），零假设为处理间基因表达变化< |1.5|倍，Benjamini–Hochberg调整错误发现率^p ≤ 0.05。成功测试的^{C. elegans}基因数量通过DeSeq2独立过滤，范围在16,785–18,292之间，取决于特定的 pairwise 处理比较。使用DeSeq2标准化和rlog转换（选项blind = False）的基因计数在R中进行主成分分析检查样品聚类按处理，并使用ggplot2（v. 3.4.4）绘图（Wickham，2016）。

使用WGCNA（v. 1.72–1）创建与实验处理类型相关的宿主基因共表达网络模块（Langfelder和Horvath，2008）（补充数据S1）。通过去除在32个样品中超过28个零计数的低表达基因来减少虚假测试，留下18,782个分析基因。使用rlog转换表达数据的 biweight midcorrelations（'bicor'）构建 signed-hybrid 基因网络，使用 binary categorical trait data 的推荐选项。选择软功率阈值5，其中拓扑拟合> 0.9且平均连接性< 1000，符合WGCNA包的标准使用（Langfelder和Horvath，2008）。为识别感兴趣模块，保守要求（i）模块 eigengene 与LC、LM或共感染处理类型有显著（^p ≤ 0.05）且强（相关性≥ |0.5|）的关联，且（ii）该模块与对比组样品没有相同符号（正或负相关）的显著相关（^p ≤ 0.05）。LC或LM模块的对比组是对照处理。对于共感染，对比是LM。还使用WGCNA包计算了每个模块 eigengene 解释的基因表达方差量。

为利用RNAseq数据进行功能推断，假设基因转录与翻译和蛋白质功能相关。mRNA和蛋白质之间的这种联系已被证明用于响应实验处理的DEGs（Koussounadis等人，2015）。尽管在某些情况下，由于转录后/翻译调控，基因转录可能与翻译关系不佳（Vogel和Marcotte，2012）。选择对线虫群体（重复平板）而非个体进行测序是基于实际约束， favoring 更高功率检测DEGs（Wang等人，2022），并最好地复制了启发本研究的实验设计（Bates等人，2021）。

2.3 感兴趣基因的GO term富集分析

使用基因集富集分析（GSEA）与基因表达排名或超几何富集分析与基因分类（Benjamini-Hochberg多重测试调整FDR ^p ≤ 0.05）进行多个GO term富集分析（The Gene Ontology Consortium等人，2023）。GSEAs使用由DEseq2计算的 fold-change 表达排名的基因。基准测试研究表明GSEA是一种高灵敏度方法，适用于探索性分析（Mathur等人，2018），尽管其他测试建议谨慎对待假阳性（Tarca等人，2013）。使用WGCNA模块成员资格在超几何富集分析中相对于背景转录组分类感兴趣基因。对于所有富集分析，使用clusterProfiler（v. 4.12.0）（Yu等人，2012）。^{C. elegans}的基因注释从R包org.Ce.eg.db（v. 3.19.1）获取（Carlson，2019）。这里报告了富集的生物过程GO terms，并在数据S2中包括了分子功能和细胞组件本体的额外结果。将富集的GO terms可视化为通过相关性连接的注释网络。使用Revigo（v. 1.8.1）通过仅显示冗余、密切相关 terms 的代表性 term（0.7相似性阈值）来减少视觉杂乱（Supek等人，2011）。GO term网络由Revigo产生，并在Cytoscape（v. 3.10.2）中最终确定（Shannon等人，2003）。网络布局通过边权重进行 prefuse force directed，然后手动调整以提高可读性和分组孤立节点。所有图最终在Inkscape（v. 1.3.2）中完成。

2.4 细菌基因组组装

从HiFi PacBio序列为^{L. celer}和^{L. musarum}细菌（革兰氏阳性，Phylum Actinomycetota）产生^{de novo}基因组组装。对于两种物种，从LB琼脂平板选择单菌落，在15 mL LB中在25°C轻摇培养3天至饱和培养。然后通过离心收集细胞，并使用NEB Monarch Genomic DNA提取试剂盒（T3010）和协议提取DNA，附加 overnight 化学裂解步骤使用溶菌酶和 mutanolysin 消化细胞壁。样品质量由利物浦大学基因组研究中心（CGR）验证和测序。两种^Leucobacter物种样品在Sequel II仪器上 multiplexed，为每种^{L. celer}和^{L. musarum}产生超过650,000 Q20+ reads，中位长度 above 18 kbp。然后使用seqtk（v. 1.3）将每组测序数据随机 down-sampled 到估计100倍覆盖度的read数据。组装前的基因组大小估计基于两种^{L. celer}（RefSeq GCF_001273835.1）和^{L. musarum}（RefSeq GCF_001273845.1）的已发表短读组装（Clark和Hodgkin，2015），并通过我们的read数据使用jellyfish（v. 2.3.0）和GenomeScope（v. 1.0）的 k-mer计数数据分析证实。

为组装reads，使用Flye（v. 2.9.3）默认设置并抑制这些单倍体细菌组装的替代 contigs。对每种物种进行两次额外的随机 down-sampling 和组装迭代，并验证组装 regardless 使用的随机read子集一致。使用BUSCO（v. 5.6.1）（参考数据库 micrococcales_odb10）、PGAP（v. 2023-10-03.build7061）、eggNOG-Mapper（v. 2.1.12）、InterProScan（v. 5.66–98.0）、geNomad（v. 1.8.0）、PathoFact（v. 1.0，ORF分支）、MMseqs2（v. 13）搜索蛋白质匹配到毒力因子数据库（VFDB，v. July 2024）、SOCfinder（v. 1.2）和FeGenie（v. 1.2）注释和评估组装（表S1和数据S3）（Belcher等人，2023；Camargo等人，2023；Cantalapiedra等人，2021；de Nies等人，2021；Garber等人，2020；Jones等人，2014；Li等人，2021；Liu等人，2022；Manni等人，2021；Steinegger和Söding，2017）。BUSCO使用“通用”直系同源物作为基因组组装完整性的基准（Manni等人，2021）。PGAP、eggNOG-Mapper和InterProScan是注释基因模型、描述、蛋白质结构域和基因本体（GO）terms的成熟方法（Cantalapiedra等人，2021；Jones等人，2014；Li等人，2021）。使用染色体、质粒和病毒标记基因的策划数据库，geNomad检测可能的原病毒和质粒序列，产生总结置信度量“病毒分数”和“质粒分数”（标度0–1）（Camargo等人，2023）。PathoFact利用毒力蛋白质模型（主要来自VFDB）的自定义数据库和随机森林分类器识别可能的毒力因子，如粘附、表面标记、分泌 machinery 和毒素蛋白质（de Nies等人，2021）。考虑由蛋白质模型匹配和随机森林子分析两者选择的PathoFact命中具有“完全”支持，而“部分”命中仅由一种支持。对于我们用MMseqs2直接查询VFDB，使用实验支持的核心毒力因子列表（Liu等人，2022）。为注释可能的“社交基因”以提供这些细菌物种如何相互作用的线索，用SOCfinder注释基因。继SOCfinder的 siderophore相关注释之后，用FeGenie进一步识别可能的 siderophore合成基因，并通过MMseqs2将FeGenie输出与其他基因注释关联。仅使用FeGenie注释推定的 siderophore合成基因，因为SOCfinder未检测到任何 siderophore合成基因簇。

3 结果

3.1 宿主对共感染的反应由更有害的寄生虫驱动

多种分析显示宿主基因表达基于采样时间和^{L. musarum}存在的强烈变化，由^{L. celer}引起的宿主变化较少。宿主基因表达数据的PCA显示时间（由PC1聚类，63%）和感染（由PC2聚类，18%）是处理变异的主要驱动因素（图1A）。在早期（10小时）和晚期（20小时）时间点，PC2将包含^{L. musarum}的处理（即LM和共感染）聚类在一起， largely 与对照或LC（即^{L. celer}单独）处理分开（图1A）。广泛地，这种模式表明^{L. musarum}的存在是宿主对感染反应的主要驱动因素， compared to^{L. celer}。还使用共表达基因网络方法（WGCNA）选择与感染处理强烈相关的基因模块。十三个基因表达网络模块中的十个区分感染宿主与对照（图1B）。所有10个清楚地将^{L. musarum}（共）感染宿主与未感染对照对比，仅一个还区分LC与对照。作为在基因组水平表征变化的第三种方法，使用GSEAs。这些分析将功能GO term注释与基于处理 pairwise 比较中表达水平排名的基因组关联。再次发现LM和共感染比LC处理引起更多变化相对于对照线虫，在早期或晚期感染时间约有200更多富集GO terms（图1C和表S2）。在个体基因水平，早期和晚期对照线虫之间的DEG分析显示宿主基因表达随时间的大变化（2895 DEGs），似乎几乎 entirely 与发育相关（数据S2）。因此，仅直接比较相同采样时间点内的感染处理以避免发育时间和感染信号混淆。DEG计数突出了^{L. musarum}感染比^{L. celer}的更大影响。相对于对照，LM和共感染处理在早期（909或886 DEGs）和晚期时间（2559或2332 DEGs）产生许多DEGs；而LC处理影响 far 较小（早期和晚期时间分别为34和146 DEGs）（图1C和表S3）。

3.2 强毒单一感染促进宿主抗菌免疫 over 繁殖投资

暴露于^{L. musarum}寄生虫的宿主中最强和最一致的转录变化表明免疫和应激反应基因上调 coupled with 繁殖投资基因下调。总体，^{L. musarum}感染宿主（无论LM或共感染）在早期感染期间与对照线虫区分有 slightly under 1000 DEGs。到晚期时间点，这些DEGs增加到 over 2000（图1C）。许多这些DEGs是相同基因，甚至跨时间点。在这些重叠DEGs中，那些具有最大转录变化的包括编码具有免疫（例如上调的C50F7.5、^clec-60、^spp-21和^ugt-18）（Irazoqui等人，2010）和繁殖投资功能（例如下调的^vit-1、^vit-3和^vit-4）（Perez和Lehner，2019）的蛋白质的基因，（表S4）。防御和繁殖之间的这种张力在基因组水平进一步证明（图2）。指示强免疫反应，所有^{L. musarum}包含处理有“对革兰氏阳性细菌的防御反应”作为上调基因中最富集的生物过程GO term（数据S2）。同时，繁殖和发育GO terms如“卵子发生”或“染色体分离”在^{L. musarum}感染期间更 lowly 表达的基因中富集（数据S2）。

基因表达模块还突出了与强毒感染相关的基因组，这些基因组将LM处理宿主与对照动物区分。模块2、7、9和10都与晚期LM感染正相关（相关性= 0.5、0.5、0.5和0.6，^p = 0.002、0.004、0.003和3e-4， respectively）。下面讨论模块2和7，因为它们进一步对比LM和共感染。然而，模块9和10类似地与晚期LM或共感染处理相关，并覆盖一系列过程如刺激反应、神经元发育和碳水化合物代谢。模块1与晚期LM处理负相关，主要与代谢和线粒体功能相关（相关性= −0.6，^p = 6e-4）。

3.3 保护性寄生虫施加亚致死防御-繁殖适应度权衡

被保护性^{L. celer}单独定植的线虫宿主适度激活防御反应。早期LC感染相对于对照上调33个DEGs，包括转运蛋白/载体蛋白、转移酶、细胞色素、化学感应受体、抗菌^ilys溶菌酶和感染响应肽（Gravato-Nobre等人，2016；Omi和Pujol，2019）（数据S4）。在早期时间点仅一个DEG下调，一个推定的胶原蛋白基因（数据S4）。晚期LC感染相对于对照增加宿主免疫和应激反应信号的指示由上调DEGs（ⁿ = 128）包括那些如^clec、^spp、^tsp、DAF-16响应和CUB结构域包含蛋白基因（数据S4）（Di等人，2018；Dierking等人，2016；Ermolaeva和Schumacher，2014；Irazoqui等人，2010；Martineau等人，2021；Shivers等人，2008；Visvikis等人，2014；Yang等人，2023）。一些免疫基因在晚期LC处理期间也下调；下调DEGs（ⁿ = 18）包括^pud基因（响应DAF-2信号）、一个^clec基因和一个编码有限范围抗菌溶菌酶的基因（数据S4）（Berndt等人，2024；Ding等人，2013）。总体，由^{L. celer}诱导的宿主DEGs far 较少且 less 指示强反应 than DEGs来自LM处理（图1C和表S3）。

更广泛地，基因组水平的表达变化表明LC宿主面临权衡通过将投资从繁殖和发育（下调基因）转移到检测、防御和免疫（上调基因）相对于对照线虫。由GSEA选择的富集GO terms显示在早期和晚期LC处理中与宿主繁殖、细胞周期、生长和发育以及RNA/蛋白质处理相关的基因表达广泛减少（图3和数据S2）。在晚期LC感染期间，宿主繁殖和后代发育的减少尤其明显。相对于它们的时间匹配对照，晚期LC在下调基因中富集的繁殖相关GO terms比早期LC处理多约4倍（图3和数据S2）。LC处理还以相对于对照线虫增加的免疫和应激反应基因表达为特征（图3，数据S2）。与感觉感知和刺激检测相关的上调过程最清楚与早期LC处理相关 rather than 晚期（图3和数据S2）。总共，变化程度（GO terms数量）从对照到LC substantially less than 对照到 either LM或共感染（表S2）。类似地，仅一个WGCNA基因网络模块清楚区分对照线虫和LC宿主，模块2。该模块与早期感染期间的LC负相关，并与晚期LM正相关（图1）。模块2基因富集了与RNA处理相关的GO terms（数据S2），尽管可能在表达塑造免疫和适应度的关键基因中相关，但单独并不高度生物信息丰富。

3.4 共感染和内源性宿主防御的组合减少适应度权衡

宿主对强毒LM感染和减弱共感染的反应具有 largely 重叠的基因表达谱。仅三个DEGs被发现区分LM和共感染样品（^far-3、^lact-6和^spds-1，所有早期感染时间点并在共感染中上调）（图1C）和（数据S4）。相反，在早期和晚期时间点分别检测到304和1411个DEGs between LC和共感染处理（数据S4）。这种模式揭示了对^{L. musarum}的 largely 保守宿主转录反应在有或没有^{L. celer}保护的情况下被激活。

尽管通过直接DEG比较 essentially 不可区分，但基因组水平的变化表明LM和共感染宿主在防御和繁殖投资方面确实不同。GSEA和WGCNA方法检测到共感染和LM处理之间基因转录的微妙（ below DEG阈值）但广泛（基因组）变化（图1C）。与这种对比一致，在对照和LM处理之间识别的DEGs在共感染期间具有类似但 marginally less 极端表达值（图S1）。GSEAs显示共感染宿主比LM处理宿主更强烈表达注释有与感觉检测相关的GO生物过程 terms 的基因，在早期和，尤其，晚期感染期间（图4）。感觉过程也被识别为早期LC期间的响应，相对于对照线虫（见上文，图3）（数据S2）。这些LC相关的化学感应注释包括“化学刺激的感觉感知”和“嗅觉行为”，由许多在共感染期间激活的相同化学感应基因支持（图S2），表明^{L. celer}在促进保护性宿主刺激检测中的可能作用。

通过晚期共感染，宿主出现保护迹象，减少免疫基因表达和增加发育基因表达（与单独被^{L. musarum}寄生虫感染的线虫相比）（图4）。例如，在共感染期间，GO term“先天免疫反应”在下调基因中富集，包括那些编码推定的免疫效应物如溶菌酶、C-type lectin-like结构域蛋白、saposins和caenopores以及neuropeptide-like蛋白的基因。同时，这些共感染宿主更高表达驱动繁殖GO terms如“胚胎发育”和“繁殖”富集的基因。许多这些基因注释有功能如精子激活、真皮细胞命运、位置发育（Hox）和减数分裂（数据S2）。晚期共感染基因表达的差异将^{L. celer}介导的保护与宿主从免疫反应分配更多资源到持续发育和繁殖的能力联系起来。

基因网络模块还区分了更有害寄生虫单独（LM处理）和较低毒力共感染的影响。尽管LM和共感染共享与许多模块的相关性，到晚期感染时间点，四个模块区分LM和