Adaptation to free-living drives loss of beneficial endosymbiosis through metabolic trade-offs

221 Isolation of symbiotic algal strains 222 We isolated four algal symbionts from lab strains of host Paramecium bursaria: two Chlorella 223 variabilis strains (Symb-1660/21-IV from host 1660/21, and Symb-HZ75.5-VI from host 224 HZ75.5) and two Micractinium conductrix strains (Symb-1660/37-III from host 1660/37, and 225 Symb-186b-X from host 186b). Prior to isolation, P. bursaria cultures were maintained in 226 standard Paramecium growth conditions at 20°C under 12 µE m-2s-1 light (10:14 L:D cycle) in 227 modified NCL medium (58) with 0.25 g Protozoan pellet (CBA053, Blades biological LTD) 228 replacing the cereal grass leaves. 229 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 19, 2026. ; https://doi.org/10.64898/2026.04.16.718893doi: bioRxiv preprint 230 Symbiotic algae were isolated by washing Paramecium cells with NCL with ampicillin, lysing 231 cells by sonication (20% power, 8 s), and plati ng the lysate on either Modified Bold Basal 232 Media (MBBM; CCAP, UK), or modified artificial WC media (MWC; (59)) supplemented with 233 amino acids, each with 1.5% agar and ampicillin. Plates were incubated at 25 °C under 50 234 µE m-2s-1 light (10:14 L:D cycle). After 3—4 weeks, individual colonies were moved to stan-235 dard algal growth conditions in liquid MBBM (25 °C, 50 µE m ⁻ ² s⁻ ¹ light, 14:10 L:D cycle, 110 236 rpm shaking). Strain identity was confirmed by sequencing 18s rRNA and ITS2 regions (60). 237 238 Selection in different nitrogen sources 239 In preparation for the selection experiment, fresh isolates were obtained as above. Six colo-240 nies were picked from each as founders for each biological replicate and moved to standard 241 algal growth conditions in liquid MBBM. Each culture was pre-adapted for 3 days for their 242 respective treatments, before being initiated at ~1.2 × 10/i3 cells mL⁻ ¹ in BBM media with one 243 of three nitrogen sources, each providing 0.0088 M nitrogen: BBM-NO ₃ (nitrate), BBM-ARG 244 (arginine), or BBM-GLN (glutamine). 245 246 Populations were serially transferred weekly for 30 weeks by subculturing 25% of the total 247 volume. Cell density was measured before and after each transfer using flow cytometry (Cy-248 toFLEX S, Beckman Coulter). 249 250 Nitrogen use assay 251 To test nitrogen use of the evolved strains, the cultures were spun down (2000 rcf, 5 min), 252 washed in BMM buffer (MBBM lacking sodium nitrate and peptone), and resuspended in 20 253 mL MBBM under standard growth conditions for 9 days to increase biomass. 10 mL of each 254 culture was then pelleted and washed twice in BBM-buffer (final wash at 10,000 xg, 3 min). 255 Pellets from each symbiont strain and selection treatment were then resuspended in either 256 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 19, 2026. ; https://doi.org/10.64898/2026.04.16.718893doi: bioRxiv preprint BBM-NO₃ , BBM-ARG, or BBM-GLN, adjusted to OD ₇₅₀ = 0.1, and transferred into 96-well 257 plates in three technical replicates. OD ₇₅₀ was measured at day 0 and day 7 using a plate 258 reader (CLARIOstar Plus, BMG LABTECH). 259 260 Symbiotic capacity assay 261 A common garden host, P. bursaria lab strain HA1, was used to test the symbiotic capacity 262 of the evolved strains. HA1 cells were cured of their native symbionts by incubating in NCL 263 with 20 µg mL -1 paraquat and 20 µg mL -1 cycloheximide in high light (50 µE m ⁻ ² s ⁻ ¹ light, 264 14:10 L:D cycle, 25 °C) for 12 days followed by 3–5 days in darkness and 5–7 days recovery 265 in standard Paramecium conditions. Cells were inspected by microscopy to confirm the ab-266 sence of symbionts and then reinfected with the evolved algal strains. 267 268 For reinfection, evolved algal strains were pre-cultured for 14 days in MBBM, pelleted, and 269 acclimated in NCL medium for 2 days. Symbiont-free HA1 hosts (10–15 cells) were then co-270 cultured with ~1.5 × 10 /i3 algal cells in 1.5 mL NCL medium under stan-271 dard Paramecium conditions. Four control cultures wi thout algae confirmed that hosts did 272 not re-establish symbiosis with surviving native symbionts. After 14 days, re-infected HA1 273 were transferred to NCL supplemented with bacterial food ( Serratia marcescens) and accli-274 mated for 5 days before experimental assays. 275 276 The fitness effect for the HA1 hosts was quantified by directly competing the re-infected 277 hosts with symbiont-free hosts in darkness (< 3 µE m ⁻ ² s⁻ ¹) and high light (50 µE m⁻ ² s⁻ ¹). 10 278 cells of each were cultured in 1.5 mL of NCL (supplemented with S. marcescens on days 0 279 and 3) for 7 days. The final number of symbiotic and symbiont-free cells were quantified with 280 flow cytometry (CytoFLEX S, Beckman Coulter). 281 282 Sugar release assay 283 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 19, 2026. ; https://doi.org/10.64898/2026.04.16.718893doi: bioRxiv preprint To quantify maltose and glucose secretion of evolved Symb-186b lines were cultured in ei-284 ther BBM-NO₃ , BBM-ARG, or BBM-GLN without sucrose in standard growth conditions on a 285 shaker (120 rpm) for 11 days to ensure cells were in exponential growth. Algal cultures were 286 pelleted (2000 rcf, 15 min), washed with BBM-buffer and pelleted again (2000 rcf, 10 min). 287 Three replicate 2 ml populations of 1.3 × 10 8 cells mL-1 of each strain were resuspended (at 288 equal concentration verified with flow cytometry) in pH 5.5 MBBM medium without N or su-289 crose (pH adjusted with 100 mM MES monohydrate and MES Na salt buffers) and incubated 290 for 6 hours at standard growth conditions. After incubation all samples were filtered (0.22 291 µm) and the supernatant frozen (-20°C) for analysis with ion chromatography. 292 293 The sugar release samples were analysed with DIONEX ICS-6000 SP, AS-AP ion chroma-294 tography at Manchester Institute of Biotechnology (MIB, Mass Spectrometry & Separations 295 Facility). The amount of maltose and glucos e was determined using a CarboPac PA20 296 Guard Column (30 mm) and a CarboPac PA20 Analytical Column (150 mm). Multi-step gra-297 dients of 45 min were created to optimise peak separation of both sugars, using 30 mM 298 NaOH as the mobile phase with a sample injection volume of 10 µl. Sugar concentrations 299 were then quantified using 6 conc entrations (STD1, 5, 10, 50, 100, 250 ppm) of known stan-300 dards. 301 302 Phenotype analysis: phenotypes at transfer, nitrogen use and sugar release 303 All below analyses were carried out using R (v . 4.4.1; (61)), unless otherwise specified. 304 Mixed-effects models were fit with lme4::lmer()(v. 1.1.37; (62)) to characterise each 305 evolved strain’s phenotype at each transfer as we ll as sugar release and nitrogen use at the 306 experiment end point. For cell density, per cell chlorophyll fluorescence, cell size and cell 307 granularity at transfer, raw flow cytometry data was imported into R using flowCore (v. 308 2.16.0; (63)) and gated using the FITC.H channel (> 1000) to ensure only live cells were in-309 cluded. log(cells/ml), PerCP.A (chlorophyll fluor escence), FSC.A (forward scatter, cell size) 310 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 19, 2026. ; https://doi.org/10.64898/2026.04.16.718893doi: bioRxiv preprint and SSC.A (side scatter, cell granularity) were fit as responses with transfer number, symbi-311 ont and treatment as the fixed effects. For s ugar release, maltose and glucose concentra-312 tions in ppm were fit as separ ate response variables with treatment as the fixed effect. For 313 growth rate in different growth media, gr owth rate was calculated as log(OD750 at day 314 7/OD750 at day 0)/7 days and fit as the response variable with treatment, test medium and 315 symbiont as the fixed effects. In all three ca ses, biological replicate was fit as a random ef-316 fect with a variable intercept. Significance testing was done using car::Anova() (v. 3.1.3; 317 (64)), and post hoc testing of marginal means as well as slopes and intercepts was carried 318 out using the emmeans package (v.1.11.2; (65)). 319 320 Phenotype analysis: symbiotic capacity 321 To obtain cell counts, the competition assay flow cytometry data was gated based on for-322 ward side scatter (i.e. size) to identify host cells. Reinfected cells were distinguished from 323 symbiont-free cells based on single cell ch lorophyll fluorescence (excitation 488 nm, emis-324 sion 690 nm). 325 326 As a considerable proportion of the competition assay trials saw re-infected cells completely 327 outcompeted by symbiont-free cells, fixation probability of the symbiont-free phenotype was 328 tested by fitting a Bayesian generalised linear model ( brms package, v. 2.22.0; (66)) with 329 Bernoulli likelihood (logit link). Treatment, symbiont and light level were fit as fixed effects 330 with weakly informative Normal (0,2) priors, while biological replicate was fit as a random 331 effect with a variable intercept. 332 333 DNA extraction and sequencing 334 Samples were collected for DNA sequencing from the 186b selection lines at transfers 1 and 335 30. Prior to extraction, 5 mL culture was grown in 25 mL MBBM for 14 days, pelleted (30,000 336 xg, 10 min), and flash-frozen in liquid nitrogen. 337 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 19, 2026. ; https://doi.org/10.64898/2026.04.16.718893doi: bioRxiv preprint 338 DNA was extracted using a custom CTAB–ethanol precipitation method. Cells were ho-339 mogenised with 0.5 mm zirconium oxide beads in CTAB buffer, incubated with with 50 µl 340 proteinase K (65 °C, 30 min), followed by addition of 3 µl RNase A and further incubation (37 341 °C, 30 min), and then extracted with an 25:1 chloroform:isoamyl alcohol solution. DNA was 342 pelleted (10,000 rpm, 60 min, 4 °C), washed twice with ice-cold 70% ethanol, and eluted in 343 200 µL TE buffer. Sequencing was carried out by CGR using the Illumina NovaSeq 6000 344 platform, using TruSeq PCR-free kit using a single lane of a S4 flow cell. 345 346 Variant Calling and Allele Frequency Estimation 347 Paired-end Illumina reads were adapter and qualit y trimmed using Trim Galore! (v0.6; (67)) 348 using standard settings. Following this, trimmed reads for each biological replicate (and in-349 cluding the ancestral library) were aligned to the Micractinium conductrix 186b nuclear ge-350 nome assembly (19) using minimap2 (v2.28; (68)), with secondary alignments suppressed. 351 Read group information indicating growing media condition and replicate number was added 352 during alignment using SAMtools (v1.21; (69) ). The next step included PCR and optical du-353 plicate removal using GATK MarkDuplicates (v4.6.1; (70)), with duplicates retained in the 354 output but marked to allow downstream filtering. 355 356 Variant discovery was performed using the GATK HaplotypeCaller in Genomic Variant Call 357 Format (GVCF) mode (this allows records for all sites), with ploidy set to 1 (to reflect the 358 haploid nature of the algal genome) and with one GVCF file per biological replicate. Then the 359 individual GVCFs were combined per experiment al condition using GATK CombineGVCFs, 360 continuing with the ancestral sample being treated as an additional condition. 361 362 Joint genotyping was performed individually us ing GATK GenotypeGVCFs using the com-363 bined GCVF per condition and the reference genom e. Finally, variant filtration was applied 364 using GATK VariantFiltration; variant sites were filtered out if they showed low quality by 365 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 19, 2026. ; https://doi.org/10.64898/2026.04.16.718893doi: bioRxiv preprint depth (QD 60.0 or SOR > 4.0), poor mean mapping quality (MQ < 40.0), or insufficient read depth 367 (DP < 10). The final passing variants were extr acted using BCFtools (v1.21; Danecek et al 368 2021), and then per-sample genotype (GT) and allele depth (AD) fields were exported to a 369 TSV file using bcftools query. 370 371 Allele Frequency Calculation and Ancestral Filtering 372 Analyses for this section were performed in R (v4.4; (61)) using the tidyverse packages 373 (71). The per-sample allele depth tables were reshaped from wide to long format, and allele 374 frequencies (AF) were computed for each variant-sample combination as the ratio of ALT-375 supporting reads to total read depth. Sites with fewer than 10 total reads were excluded. To 376 remove potential genetic variation present prior to experimental evolution, any variant posi-377 tion at which the ancestral sample exhibited an ALT allele frequency ≥ 0.1 was identified and 378 then discarded from all evolved condition datasets. 379 380 Functional Annotation 381 Variants passing ancestral filtering were func tionally annotated using snpEff (v 5.2; (72)), 382 with a custom database constructed from the M. conductrix 186b genome assembly, CDS 383 sequences, protein sequences, and GFF3 gene model s. The ANN field from the snpEff VCF 384 output was parsed to recover effect class, predicted impact, gene name, gene identifier, and 385 HGVS nomenclature for each alternate allele. Intergenic variants were excluded from down-386 stream analyses. 387 388 GO term mappings were derived from InterProScan (v 5.73; (73)) output, run against all 389 available member databases with GO term and pathway lookup enabled. GO terms were 390 extracted from the InterProScan TSV output and collapsed to a unique set per gene identi-391 fier. 392 393 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 19, 2026. ; https://doi.org/10.64898/2026.04.16.718893doi: bioRxiv preprint Genomic parallelism and candidate gene identification 394 Initial analyses were focussed on fixed or near-fixed variants in the peak of frequencies be-395 tween 0.95 and 1 (Fig. S7), with parallelism first assessed across the full dataset using these 396 high-frequency variants. Presence–absence matr ices were constructed and Jaccard dis-397 tances were calculated using the vegan package. Genome wide parallelism was tested by 398 comparing within- versus between-treatment distances using permutation tests (n = 9999) 399 and bootstrap resampling. 400 401 To identify candidate loci underlying treatment-specific divergence, the data set was ana-402 lysed separately for each treatment using multiple complementary approaches. First, 403 Fisher’s exact tests were applied to identify individual SNPs or genes disproportionately as-404 sociated with one treatment over the others. Second, constrained or dinations (CAP) were 405 performed using vegan::capscale(), with alleles/genes fitted as explanatory vectors 406 via envfit() (v. 2.7.1; (74)); loci with significant loadings were ranked by effect size (r²). 407 Third, an indicator species anal ysis was conducted using the indicspecies package (v. 408 1.8.0; (75)) to detect sets of loci indicative of a treatment. This was done across a range of 409 allele frequency thresholds (1, 0.99, 0.97, 0.95, 0.90, 0.85, 0.80, 0.75) for determining fixed 410 variants to ensure the results were robust to threshold changes. The output from these mod-411 els was joined and filtered for variants unique to each treatment and parallel in at least 2 rep-412 licate populations and present in at least 75% of runs across the allele frequency thresholds 413 to form our list of genes of interest. 414 415 Candidate gene function identification 416 To summarise functional annotations, shortlisted genes were linked to Gene Ontology (GO) 417 terms and analysed for GO enrichment against the full list of GO terms from the reference 418 annotation using the topGO package (v. 2.56.0; (76)). Semantic similarity of enriched GO 419 terms was calculated using the rrvgo package (v. 1.16.0; (77)) fo r plotting. Candidate pro-420 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 19, 2026. ; https://doi.org/10.64898/2026.04.16.718893doi: bioRxiv preprint tein sequences were further annotated by m apping to KEGG orthologs using GhostKOALA 421 (https://www.kegg.jp/ghostkoala/) and by sequence similarity searches against the SwissProt 422 Viridiplantae database using blastp through the NCBI web tool 423 (https://blast.ncbi.nlm.nih.gov). IDs for the blastp output were mapped using the Uniprot web 424 tool (https://www.uniprot.org/id-mapping). 425 426 Untargeted metabolic fingerprinting using DESI-MSI 427 Culture samples were pelleted by centrifuging (2000 rcf, 5 min), weighed and flash frozen at 428 transfers 1, 5, 10, 20 and 30 and stored in -80 °C. Samples were dissolved in 1 mL of a 429 methanol-water mix (4:1) for 10 minutes before being centrifuged (10 mins, 16 163 xg) and 430 the supernatant removed. The supernatant was put in a SpeedVac (Thermo Fisher, UK) for 431 3 hours 30 minutes, until samples were completely dry. Samples were then reconstituted 432 with 80:20 methanol water solvent at a ratio of 10 µL solvent per 1 mg wet mass. 433 434 2 µl of each sample was spotted onto Waters PTFE coated microscope glass slides (Waters 435 Corporation, Wilmslow, UK) in a randomised block design (3 reps) and were analysed on a 436 SYNAPT XS mass spectrometer (Waters Corpor ation, Wilmslow, UK) coupled with a modi-437 fied DESI XS source in negative ionisation mode. DESI solvent spray was composed of a 438 methanol-water mixture (98:2) at a flow rate of 1.5 μ L/min with nitrogen gas flow set to 0.5—439 1 bar, a capillary voltage of 0.57 kV with 25 V sampling cone, heated transfer line at 450 °C 440 and source temperature of 100 °C; the trap and transfer cell voltages were set at 4 and 2 V, 441 respectively. We used a resolution of 150 µm moving at 1500 µm/s. The mass spectrometer 442 was operated using MassLynx v4.2 (Waters Corporation, Wilmslow, UK). Deposited spot 443 regions were extracted and pre-processed (0.2 Da window, m/z between 50-2000 were se-444 lected, top 2000 peaks) using HDI software (v1.5) . Data were further processed using R for 445 lock mass correction, clustering of peaks with a 10 ppm tolerance and calculating median 446 weighted centroids, requiring presence in at least 2/3 technical replicates and keeping only 447 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 19, 2026. ; https://doi.org/10.64898/2026.04.16.718893doi: bioRxiv preprint features above a minimum intensity threshold (sample mean > median of all means in sam-448 ple) to reduce downstream dimensionality. 449 450 Identification of candidate features 451 To test the distribution of m/z separated features, the m/z intensity matrix was log + 1 trans-452 formed as well as Pareto scaled and a PCA was fit using vegan::rda() (v. 2.7.1; (74)). To de-453 termine whether clustering was explained by treatment, a PERMANOVA was fit using ve-454 gan::adonis2() (v. 2.7.1; (74)). 455 456 To identify m/z features of interest, PCAs were performed on the full dataset as well as sub-457 sets comprising a single treatment or a singl e timepoint and the top 100 loadings were ex-458 tracted. Furthermore, as a complementary approach, further m/z peaks of interest were 459 identified by using a random forest model fit with ranger (v. 0.17.0; (78)) using 70% of the 460 data frame at each permutation to allow for estimating permutation importance of each fea-461 ture (n = 1000), selecting features present in the in top 2x of the elbow of the mean impor-462 tance curve in at least 50% of all permutations. An LMM (intensity ~ treatment * transfer + 463 (1|bio_rep)) was fit to each m/z of interest lme4::lmer()(v. 1.1.37; (62)), and m/z peaks 464 where treatment or the interaction with transfer had a significant effect were shortlisted. The 465 shortlist was further filtered based on there bei ng a significant difference as well as a mini-466 mum 1.5 fold change between at least two of the treatments by the final transfer. 467 468 MS/MS methods 469 The top features from the statistical analysis were targeted for DESI-MS/MS analysis. Data 470 was acquired using MassLynx V4.2 (Waters Cor poration, Wilmslow, UK) by rastering over 471 the sample slides. Trap collision voltage was optimised for each quadruple isolated m/z 472 value, ranging from 4-45 V, with all other MS conditions maintained as those used during the 473 full scan DESI-MS experiment. Putative ID’s were assigned using product ion m/z and 474 MS/MS fragment, comparing to literature (79–83), the LIPIDMAPS and HMBD database. 475 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 19, 2026. ; https://doi.org/10.64898/2026.04.16.718893doi: bioRxiv preprint Characteristic headgroup fragments were used to confirm lipid class, m/z 225 [C 6H9O7S]− 476 for sulfonquinovosyl diagcylglycerol (SQDG) (79), m/z 241 [C 6H10O8P]- for phosphatidy-477 linositol (PI) and m/z 153 [C 3H6O5P]- for all glycerophospholipids (84). Additionally, fatty 478 acyl chain fragments and neutral loss were used for further confirmation. MS/MS data for 479 features that could not be assigned an identity can be found in Supplementary Materials 3. 480 481

482 1. Wernegreen JJ. Endosymbiosis. Curr Biol. 2012 Jul;22(14):R555–61. 483 doi:10.1016/j.cub.2012.06.010 484 2. Zook DP. Prioritizing Symbiosis to Sustain Biodiversity: Are Symbionts Keystone 485 Species? In: Seckbach J, editor. Symbiosis [Internet]. Dordrecht: Springer Nether-486 lands; 2001 [cited 2025 Nov 4]. p. 3–12. (Cellular Origin, Life in Extreme Habitats 487 and Astrobiology). Available from: https://link.springer.com/10.1007/0-306-48173-488 1_1 doi:10.1007/0-306-48173-1_1 489 3. Sudakaran S, Kost C, Kaltenpoth M. Symbiont Acquisition and Replacement as a 490 Source of Ecological Innovation. Trends Microbiol. 2017 May;25(5):375–90. 491 doi:10.1016/j.tim.2017.02.014 492 4. López-García P, Moreira D. The symbiotic origin of the eukaryotic cell. C R Biol. 493 2023 May 30;346(G1):55–74. doi:10.5802/crbiol.118 494 5. Kiers ET, West SA. Evolving new organisms via symbiosis. Science. 2015 Apr 495 24;348(6233):392–4. doi:10.1126/science.aaa9605 496 6. Powell JR, Rillig MC. Biodiversity of arbuscular mycorrhizal fungi and ecosystem 497 function. New Phytol. 2018 Dec;220(4):1059–75. doi:10.1111/nph.15119 498 7. Frank SA. Models of Symbiosis. Am Nat. 1997 Jul;150(S1):S80–99. 499 doi:10.1086/286051 500 8. Law R, Dieckmann U. Symbiosis through exploitation and the merger of lineages 501 in evolution. Proc R Soc Lond B Biol Sci. 1998 Jul 7;265(1402):1245–53. 502 doi:10.1098/rspb.1998.0426 503 9. Brockhurst MA, Cameron DD, Beckerman AP. Fitness trade-offs and the origins of 504 endosymbiosis. PLOS Biol. 2024 Apr 12;22(4):e3002580. 505 doi:10.1371/journal.pbio.3002580 506 10. Sørensen MES, Wood AJ, Minter EJA, Lowe CD, Cameron DD, Brockhurst 507 MA. Comparison of Independent Evolutionary Origins Reveals Both Convergence 508 and Divergence in the Metabolic Mechanisms of Symbiosis. Curr Biol. 2020 509 Jan;30(2):328-334.e4. doi:10.1016/j.cub.2019.11.053 510 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 19, 2026. ; https://doi.org/10.64898/2026.04.16.718893doi: bioRxiv preprint 11. He M, Wang J, Fan X, Liu X, Shi W, Huang N, et al. Genetic basis for the es-511 tablishment of endosymbiosis in Paramecium. ISME J. 2019 May 1;13(5):1360–9. 512 doi:10.1038/s41396-018-0341-4 513 12. Wetzel RG. Limnology [Internet]. Elsevier; 2001 [cited 2026 Apr 8]. Available 514 from: https://linkinghub.elsevier.com/retrieve/pii/C20090021126 515 doi:10.1016/C2009-0-02112-6 516 13. Fujishima M, Kodama Y . Mechanisms for establishing primary and secondary 517 endosymbiosis in Paramecium. J Eukaryot Microbiol. 2022 Sep 19;69(5). 518 doi:10.1111/jeu.12901 519 14. Jenkins BH. Mutualism on the edge: Understanding the Paramecium–520 Chlorella symbiosis. PLOS Biol. 2024 Apr 4;22(4):e3002563. 521 doi:10.1371/journal.pbio.3002563 522 15. Vitonyt ė I, Hansson EM, Malumphy-Montesdeoca D, Leonard G, Savory FR, 523 Dunning LT, et al. Intracellular niche specialisation drives evolutionary entrapment 524 of endosymbiotic algae [Internet]. Evolutionary Biology; 2025 [cited 2026 Apr 8]. 525 Available from: http://biorxiv.org/lookup/doi/10.64898/2025.12.22.695913 526 doi:10.64898/2025.12.22.695913 527 16. Lowe CD, Minter EJ, Cameron DD, Brockhurst MA. Shining a Light on Ex-528 ploitative Host Control in a Photosynthetic Endosymbiosis. Curr Biol. 2016 529 Jan;26(2):207–11. doi:10.1016/j.cub.2015.11.052 530 17. Dorling M, McAuley PJ, Hodge H. Effect of pH on growth and carbon metabo-531 lism of maltose-releasing Chlorella (Chlorophyta). Eur J Phycol. 1997 532 Feb;32(1):19–24. doi:10.1080/09541449710001719335 533 18. Kodama Y, Fujishima M. Secondary Symbiosis Between Paramecium and 534 Chlorella Cells. In: International Review of Cell and Molecular Biology [Internet]. 535 Elsevier; 2010 [cited 2026 Apr 8]. p. 33–77. Available from: 536 https://linkinghub.elsevier.com/retrieve/pii/S193764481079002X 537 doi:10.1016/S1937-6448(10)79002-X 538 19. Leonard G, Vitonyte I, Savory FR, Hansson EM, Cameron DD, Brockhurst M, 539 et al. De novo genome sequence assembly of the model algal endosymbiont Mi-540 cractinium conductrix derived from its host Paramecium bursaria 186b [Internet]. 541 Genomics; 2025 [cited 2026 Apr 8]. Available from: 542 http://biorxiv.org/lookup/doi/10.64898/2025.12.15.694297 543 doi:10.64898/2025.12.15.694297 544 20. Moreira E, Coimbra S, Melo P . Glutamine synthetase: an unlikely case of 545 functional redundancy in Arabidopsis thaliana. Wicke S, editor. Plant Biol. 2022 546 Aug;24(5):713–20. doi:10.1111/plb.13408 547 21. Miyawaki K, Tarkowski P , Matsumoto-Kitano M, Kato T, Sato S, Tarkowska D, 548 et al. Roles of Arabidopsis ATP/ADP isopentenyltransferases and tRNA isopen-549 tenyltransferases in cytokinin biosynthesis. Proc Natl Acad Sci. 2006 Oct 550 31;103(44):16598–603. doi:10.1073/pnas.0603522103 551 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 19, 2026. ; https://doi.org/10.64898/2026.04.16.718893doi: bioRxiv preprint 22. Garbarino JE, Gibbons IR. Expression and genomic analysis of midasin, a 552 novel and highly conserved AAA protein distantly related to dynein. BMC Genom-553 ics. 2002 Jul 8;3(1):18. doi:10.1186/1471-2164-3-18 554 23. Apel K, Hirt H. Reactive Oxygen Species: Metabolism, Oxidative Stress, and 555 Signal Transduction. Annu Rev Plant Biol. 2004 Jun 2;55(1):373–99. 556 doi:10.1146/annurev.arplant.55.031903.141701 557 24. Singla ‐ Pareek SL, Kaur C, Kumar B, Pareek A, Sopory SK. Reassessing plant 558 glyoxalases: large family and expanding functions. New Phytol. 2020 559 Aug;227(3):714–21. doi:10.1111/nph.16576 560 25. El-Sappah AH, Li J, Yan K, Zhu C, Huang Q, Zhu Y, et al. Fibrillin gene family 561 and its role in plant growth, development, and abiotic stress. Front Plant Sci. 2024 562 Oct 29;15:1453974. doi:10.3389/fpls.2024.1453974 563 26. Durante L, Hübner W, Lauersen KJ, Remacle C. Characterization of the 564 GPR1/FUN34/YaaH protein family in the green microalga Chlamydomonas sug-565 gests their role as intracellular membrane acetate channels. Plant Direct. 2019 566 Jun;3(6):e00148. doi:10.1002/pld3.148 567 27. Pao SS, Paulsen IT, Saier MH. Major Facilitator Superfamily. Microbiol Mol 568 Biol Rev. 1998 Mar;62(1):1–34. doi:10.1128/MMBR.62.1.1-34.1998 569 28. Yan N. Structural investigation of the proton-coupled secondary transporters. 570 Curr Opin Struct Biol. 2013 Aug;23(4):483–91. doi:10.1016/j.sbi.2013.04.011 571 29. Mollá ‐ Morales A, Sarmiento‐ Mañús R, Robles P , Quesada V, Pérez‐ Pérez JM, 572 González‐ Bayón R, et al. Analysis of ven3 and ven6 reticulate mutants reveals the 573 importance of arginine biosynthesis in Arabidopsis leaf development. Plant J. 2011 574 Feb;65(3):335–45. doi:10.1111/j.1365-313X.2010.04425.x 575 30. Gaufichon L, Reisdorf-Cren M, Rothstein SJ, Chardon F, Suzuki A. Biological 576 functions of asparagine synthetase in plants. Plant Sci. 2010 Sep;179(3):141–53. 577 doi:10.1016/j.plantsci.2010.04.010 578 31. Wu S, Yu Z, Wang F, Li W, Ye C, Li J, et al. Cloning, characterization, and 579 transformation of the phosphoethanolamine N-methyltransferase gene 580 (ZmPEAMT1) in maize (Zea mays L.). Mol Biotechnol. 2007 May 31;36(2):102–581 12. doi:10.1007/s12033-007-0009-1 582 32. Niu L, Lu F, Pei Y, Liu C, Cao X. Regulation of flowering time by the protein 583 arginine methyltransferase AtPRMT10. EMBO Rep. 2007 Dec;8(12):1190–5. 584 doi:10.1038/sj.embor.7401111 585 33. Kabelitz T, Brzezinka K, Friedrich T, Górka M, Graf A, Kappel C, et al. A JU-586 MONJI Protein with E3 Ligase and Histone H3 Binding Activities Affects Transpo-587 son Silencing in Arabidopsis. Plant Physiol. 2016 May;171(1):344–58. 588 doi:10.1104/pp.15.01688 589 34. Tanaka K. The proteasome: Overview of structure and functions. Proc Jpn 590 Acad Ser B. 2009;85(1):12–36. doi:10.2183/pjab.85.12 591 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 19, 2026. ; https://doi.org/10.64898/2026.04.16.718893doi: bioRxiv preprint 35. Hoecker U. The activities of the E3 ubiquitin ligase COP1/SPA, a key repres-592 sor in light signaling. Curr Opin Plant Biol. 2017 Jun;37:63–9. 593 doi:10.1016/j.pbi.2017.03.015 594 36. Cardozo T, Pagano M. The SCF ubiquitin ligase: insights into a molecular ma-595 chine. Nat Rev Mol Cell Biol. 2004 Sep 1;5(9):739–51. doi:10.1038/nrm1471 596 37. Jin H, Song Z, Nikolau BJ. Reverse genetic characterization of two paralo-597 gous acetoacetyl CoA thiolase genes in Arabidopsis reveals their importance in 598 plant growth and development. Plant J. 2012 Jun;70(6):1015–32. 599 doi:10.1111/j.1365-313X.2012.04942.x 600 38. Schelbert S, Aubry S, Burla B, Agne B, Kessler F, Krupinska K, et al. Pheo-601 phytin Pheophorbide Hydrolase (Pheophytinase) Is Involved in Chlorophyll Break-602 down during Leaf Senescence in Arabidopsis. Plant Cell. 2009 Apr 28;21(3):767–603 85. doi:10.1105/tpc.108.064089 604 39. Kang J, Park J, Choi H, Burla B, Kretzschmar T, Lee Y , et al. Plant ABC 605 Transporters. Arab Book. 2011 Jan;9:e0153. doi:10.1199/tab.0153 606 40. Dahuja A, Kumar RR, Sakhare A, Watts A, Singh B, Goswami S, et al. Role of 607 ATP‐ binding cassette transporters in maintaining plant homeostasis under abiotic 608 and biotic stresses. Physiol Plant. 2021 Apr;171(4):785–801. 609 doi:10.1111/ppl.13302 610 41. Zhou T, Yue C peng, Huang J yong, Cui J qian, Liu Y , Wang W ming, et al. 611 Genome-wide identification of the amino acid permease genes and molecular 612 characterization of their transcriptional responses to various nutrient stresses in al-613 lotetraploid rapeseed. BMC Plant Biol. 2020 Dec;20(1):151. doi:10.1186/s12870-614 020-02367-7 615 42. Wada H, Murata N. The essential role of phosphatidylglycerol in photosynthe-616 sis. Photosynth Res. 2007 Aug 2;92(2):205–15. doi:10.1007/s11120-007-9203-z 617 43. Upchurch RG. Fatty acid unsaturation, mobilization, and regulation in the re-618 sponse of plants to stress. Biotechnol Lett. 2008 Jun;30(6):967–77. 619 doi:10.1007/s10529-008-9639-z 620 44. Benning C. Biosynthesis and Function of the Sulfolipid Sulfoquinovosyl dia-621 cylglycerol. Annu Rev Plant Physiol Plant Mol Biol. 1998 Jun;49(1):53–75. 622 doi:10.1146/annurev.arplant.49.1.53 623 45. Heilmann I. Phosphoinositide signaling in plant development. Development. 624 2016 Jun 15;143(12):2044–55. doi:10.1242/dev.136432 625 46. Ramanan R, Tran QG, Cho DH, Jung JE, Kim BH, Shin SY , et al. The Ancient 626 Phosphatidylinositol 3-Kinase Signaling System Is a Master Regulator of Energy 627 and Carbon Metabolism in Algae. Plant Physiol. 2018 Jul;177(3):1050–65. 628 doi:10.1104/pp.17.01780 629 47. Horvath SE, Daum G. Lipids of mitochondria. Prog Lipid Res. 2013 630 Oct;52(4):590–614. doi:10.1016/j.plipres.2013.07.002 631 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 19, 2026. ; https://doi.org/10.64898/2026.04.16.718893doi: bioRxiv preprint 48. Yu L, Zhou C, Fan J, Shanklin J, Xu C. Mechanisms and functions of mem-632 brane lipid remodeling in plants. Plant J. 2021 Jul;107(1):37–53. 633 doi:10.1111/tpj.15273 634 49. Sanz-Luque E, Chamizo-Ampudia A, Llamas A, Galvan A, Fernandez E. Un-635 derstanding nitrate assimilation and its regulation in microalgae. Front Plant Sci. 636 2015 Oct 26;6. doi:10.3389/fpls.2015.00899 637 50. Fortunato S, Nigro D, Lasorella C, Marcotuli I, Gadaleta A, De Pinto MC. The 638 Role of Glutamine Synthetase (GS) and Glutamate Synthase (GOGAT) in the Im-639 provement of Nitrogen Use Efficiency in Cereals. Biomolecules. 2023 Dec 640 10;13(12):1771. doi:10.3390/biom13121771 641 51. Majumdar R, Barchi B, Turlapati SA, Gagne M, Minocha R, Long S, et al. Glu-642 tamate, Ornithine, Arginine, Proline, and Polyamine Metabolic Interactions: The 643 Pathway Is Regulated at the Post-Transcriptional Level. Front Plant Sci. 2016 Feb 644 16;7. doi:10.3389/fpls.2016.00078 645 52. Lavandosque LL, Vischi Winck F. Polyamine-Mediated Growth Regulation in 646 Microalgae: Integrating Redox Balance and Amino Acids Pathway into Metabolic 647 Engineering. SynBio. 2025 May 28;3(2):8. doi:10.3390/synbio3020008 648 53. Fernandez E, Galvan A. Nitrate Assimilation in Chlamydomonas. Eukaryot 649 Cell. 2008 Apr;7(4):555–9. doi:10.1128/EC.00431-07 650 54. Kato Y, Ueno S, Imamura N. Studies on the nitrogen utilization of endosymbi-651 otic algae isolated from Japanese Paramecium bursaria. Plant Sci. 2006 652 Mar;170(3):481–6. doi:10.1016/j.plantsci.2005.09.018 653 55. Pankey MS, Foxall RL, Ster IM, Perry LA, Schuster BM, Donner RA, et al. 654 Host-selected mutations converging on a global regulator drive an adaptive leap 655 towards symbiosis in bacteria. eLife. 2017 Apr 27;6:e24414. 656 doi:10.7554/eLife.24414 657 56. Koga R, Moriyama M, Onodera-Tanifuji N, Ishii Y , Takai H, Mizutani M, et al. 658 Single mutation makes Escherichia coli an insect mutualist. Nat Microbiol. 2022 659 Aug 4;7(8):1141–50. doi:10.1038/s41564-022-01179-9 660 57. Siozios S, Nadal-Jimenez P, Azagi T, Sprong H, Frost CL, Parratt SR, et al. 661 Genome dynamics across the evolutionary transition to endosymbiosis. Curr Biol. 662 2024 Dec;34(24):5659-5670.e7. doi:10.1016/j.cub.2024.10.044 663 58. Prescott DM, James TW. Culturing of Amoeba proteus on Tetrahymena. Exp 664 Cell Res. 1955 Jan;8(1):256–8. doi:10.1016/0014-4827(55)90067-7 665 59. Lebret K, Kritzberg ES, Figueroa R, Rengefors K. Genetic diversity within and 666 genetic differentiation between blooms of a microalgal species. Environ Microbiol. 667 2012 Sep;14(9):2395–404. doi:10.1111/j.1462-2920.2012.02769.x 668 60. Hoshina R, Iwataki M, Imamura N. Chlorella variabilis and Micractinium reis-669 seri sp. nov. (Chlorellaceae, Trebouxiophyceae): Redescription of the endosymbi-670 otic green algae of Paramecium bursaria (Peniculia, Oligohymenophorea) in the 671 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 19, 2026. ; https://doi.org/10.64898/2026.04.16.718893doi: bioRxiv preprint 120th year: Taxonomy of the photobionts of Paramecium. Phycol Res. 2010 Jun 672 24;58(3):188–201. doi:10.1111/j.1440-1835.2010.00579.x 673 61. R Core Team. R: A Language and Environment for Statistical Computing 674 [Internet]. Vienna, Austria: R Foundation for Statistical Computing; 2025. Available 675 from: https://www.R-project.org/ 676 62. Bates D, Mächler M, Bolker B, Walker S. Fitting Linear Mixed-Effects Models 677 Using lme4. J Stat Softw. 2015;67(1). doi:10.18637/jss.v067.i01 678 63. Ellis B, Haaland P , Hahne F, Meur NL, Gopalakrishnan N, Spidlen J, et al. 679 flowCore: flowCore: Basic structures for flow cytometry data [Internet]. 2024. 680 Available from: https://bioconductor.org/packages/flowCore 681 doi:10.18129/B9.bioc.flowCore 682 64. Fox J, Weisberg S. An R Companion to Applied Regression [Internet]. Third. 683 Thousand Oaks CA: Sage; 2019. Available from: https://www.john-684 fox.ca/Companion/ 685 65. Lenth RV, Piaskowski J. emmeans: Estimated Marginal Means, aka Least-686 Squares Means [Internet]. 2026. Available from: https://CRAN.R-687 project.org/package=emmeans 688 66. Bürkner PC. Bayesian Item Response Modeling in R with brms and Stan. J 689 Stat Softw. 2021;100(5):1–54. doi:10.18637/jss.v100.i05 690 67. Krueger F. TrimGalore! [Internet]. Available from: 691 https://github.com/FelixKrueger/TrimGalore 692 68. Li H. Minimap2: pairwise alignment for nucleotide sequences. Birol I, editor. 693 Bioinformatics. 2018 Sep 15;34(18):3094–100. doi:10.1093/bioinformatics/bty191 694 69. Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, et al. 695 Twelve years of SAMtools and BCFtools. GigaScience. 2021 Jan 696 29;10(2):giab008. doi:10.1093/gigascience/giab008 697 70. Auwera G van der, O’Connor BD. Genomics in the cloud: using Docker, 698 GATK, and WDL in Terra. First edition. Sebastopol, CA: O’Reilly Media; 2020. 467 699 p. 700 71. Wickham H, Averick M, Bryan J, Chang W, McGowan L, François R, et al. 701 Welcome to the Tidyverse. J Open Source Softw. 2019 Nov 21;4(43):1686. 702 doi:10.21105/joss.01686 703 72. Cingolani P, Platts A, Wang LL, Coon M, Nguyen T, Wang L, et al. A program 704 for annotating and predicting the effects of single nucleotide polymorphisms, 705 SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118/i2 ; iso-2; iso-706 3. Fly (Austin). 2012 Apr;6(2):80–92. doi:10.4161/fly.19695 707 73. Jones P, Binns D, Chang HY , Fraser M, Li W, McAnulla C, et al. InterProScan 708 5: genome-scale protein function classification. Bioinformatics. 2014 May 709 1;30(9):1236–40. doi:10.1093/bioinformatics/btu031 710 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 19, 2026. ; https://doi.org/10.64898/2026.04.16.718893doi: bioRxiv preprint 74. Oksanen J, Simpson GL, Blanchet FG, Kindt R, Legendre P , Minchin PR, et 711 al. vegan: Community Ecology Package [Internet]. 2025. Available from: 712 https://CRAN.R-project.org/package=vegan 713 75. Cáceres MD, Legendre P . Associations between species and groups of sites: 714 indices and statistical inference. Ecology. 2009;90:3566–74. doi:10.1890/08-715 1823.1 716 76. Alexa A, Rahnenfuhrer J. topGO: Enrichment Analysis for Gene Ontology 717 [Internet]. 2024. Available from: https://bioconductor.org/packages/topGO 718 doi:10.18129/B9.bioc.topGO 719 77. Sayols S. rrvgo: a Bioconductor package to reduce and visualize Gene Ontol-720 ogy terms. MicroPublication Biol. 2023. doi:10.17912/micropub.biology.000811 721 78. Wright MN, Ziegler A. ranger: A Fast Implementation of Random Forests for 722 High Dimensional Data in C++ and R. J Stat Softw. 2017;77(1):1–17. 723 doi:10.18637/jss.v077.i01 724 79. Da Costa E, Silva J, Mendonça S, Abreu M, Domingues M. Lipidomic Ap-725 proaches towards Deciphering Glycolipids from Microalgae as a Reservoir of Bio-726 active Lipids. Mar Drugs. 2016 May 19;14(5):101. doi:10.3390/md14050101 727 80. Rey F, Lopes D, Maciel E, Monteiro J, Skjermo J, Funderud J, et al. Polar lipid 728 profile of Saccharina latissima, a functional food from the sea. Algal Res. 2019 729 May;39:101473. doi:10.1016/j.algal.2019.101473 730 81. Conde TA, Couto D, Melo T, Costa M, Silva J, Domingues MR, et al. Polar 731 lipidomic profile shows Chlorococcum amblystomatis as a promising source of 732 value-added lipids. Sci Rep. 2021 Feb 23;11(1):4355. doi:10.1038/s41598-021-733 83455-y 734 82. Marques F, Lopes D, Conde T, Melo T, Silva J, Abreu MH, et al. Lipidomic 735 Characterization and Antioxidant Activity of Macro- and Microalgae Blend. Life. 736 2023 Jan 13;13(1):231. doi:10.3390/life13010231 737 83. Chadova K. Algal Adaptation to Environmental Stresses: Lipidomics Re-738 search. Int J Plant Biol. 2024 Jul 22;15(3):719–32. doi:10.3390/ijpb15030052 739 84. Pi J, Wu X, Feng Y. Fragmentation patterns of five types of phospholipids by 740 ultra-high-performance liquid chromatography electrospray ionization quadrupole 741 time-of-flight tandem mass spectrometry. Anal Methods. 2016;8(6):1319–32. 742 doi:10.1039/C5AY00776C 743 744 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 19, 2026. ; https://doi.org/10.64898/2026.04.16.718893doi: bioRxiv preprint