frequency of events at specific time points were carried out using *ad hoc* R scripts (R version 4.4.1) and the *tidyverse*¹⁶⁷ package (version 2.0.0). Graphs were generated using the *ggplot2*¹⁶⁶ package (version 3.5.1). For each gene, a COG functional category was retrieved from the *eggNOG* mapper output and the COG enrichment analysis was carried out in R using the *clusterProfiler*¹⁶⁵ package (version 4.6.2) by comparing each set of gene families with the full set of gene families.

## Detection of viral genome insertions and viral regions in algal genomes

To reduce the dataset size for analysis, 64 Phaeoexplorer and eight public genomes were initially filtered to retain only contigs that were more than 10 kbp in length. Gene prediction was then carried out on all contigs using *Prodigal*¹²⁶ (V2.6.3, settings: default, meta) and the resulting proteins were used as queries against the *NCVOG*¹⁷¹ and *VOGDB*¹⁸⁵ databases using *hmmscan* (HMMER 3.3.2 with default settings). The contigs detected by *hmmscan* were then filtered to retain only sequences with at least one match to either viral database at a defined e-value cutoff ($1e^{-20}$ for *NCVOG*, and $1e^{-80}$ for *VOGDB*). The resulting positive 4,951 contigs were then analysed using *ViralRecall*¹²⁷ version 2.0 with settings -w 50 -g 1 -b -f -m 2 using the built-in *Nucleocytoviricota* (NCV) database *GVOG* and a window size of 50 kbp. To ensure that viral genes were not missed because they had not been annotated by *Prodigal*, six-frame translations of the contigs were generated using *esl-translate* (version 0.48 with default settings), and the resulting proteins queried against the same databases used by *ViralRecall* using *hmmssearch* (HMMER 3.3.2, settings: -E 1e-10). The *ViralRecall* results were then parsed using an in-house workflow. Six-frame translations were removed from the results if they overlapped (even partially) with any *Prodigal* gene prediction, as identified using *bedtools*¹²⁸ (v2.29.2; intersect). Likewise, overlapping six-frame translations and gene predictions with the same *NCVOG* match were removed to reduce redundancy. Based on the distance between query sequences with the same *GVOG* hit, queries were flagged as frame-shifted (less than 100 bp gap), intron-containing (100-5,000 bp gap) or mono-exonic (greater than 5,000 bp gap). All queries were also checked for overlaps with multi-exonic genes that had been annotated by the Phaeoexplorer gene prediction procedure (using *Gmove*¹⁰⁹), and flagged if they did. All queries were then filtered to retain only those that matched a set of key NCV marker genes, identified by *NCVOG* code (A32, D5 helicase, D5 DNA primase, MCP, DNA polymerase B, SFII and VLTF3) or some Phaeovirus integrase genes (integrase recombinase, integrase resolvase and RNR). The marker gene proteins were clustered with the protein sequences of *NCVOGs* using *MMseqs* cluster¹²⁹ (version 13.45111 with settings -min-seq-id 0.3-c 0.8). Finally, the parsed results of the NCV marker gene set identified by the *ViralRecall* screen were manually curated, retaining only those queries with varying combinations of the following properties: placement within a viral region as identified by *ViralRecall*, similar *hmmssearch* results (score and e-value) and gene length to that of known NCV genes, not part of a multi-exonic gene, lack of Pfam HMM matches to cellular domains sharing homology to the marker gene (specific to certain marker genes), and clustered with an *NCVOG* in the *MMseqs* analysis. We noted that the median number of viral regions found in genomes assembled with long reads was very similar to that for genomes assembled with short reads (9 and 10, respectively). The marker gene content of the viral regions was manually assessed to estimate the number of complete or partial inserted viruses in each genome. VRs were considered to be complete proviruses if they contained all seven of the key NCV marker genes listed above. VRs were classed as partial proviruses if they only contained a subset of the seven key NCV marker genes, the presence of the MCP and DNA polymerase B genes being particularly strong indicators of a partial provirus.

To classify genes in VRs (Figure 6B), viral sequences were removed for the NBCI RefSeq non-redundant protein database (NR) by removing proteins assigned to the "Viruses" category and by comparing the database with RVDB using *BLASTp* and removing any proteins that matched with an e-value cut-off of < 1-e40 to create a "virus-free NR" database. Deduced proteins were then compared with the RVDB and the virus-free NR databases using *BLASTp* and relative bitscores (rbitscores) were calculated by dividing the *BLASTp* bitscore for the best match in each database by the query protein's self-hit bitscore.¹⁸⁶ Self-hit scores were acquired by comparing the complete deduced proteomes with themselves using *BLASTp*. Proteins with a RVDB rbitscore at least 20% greater than its virus-free NR rbitscore were designated as "viral". Proteins with a virus-free NR rbitscore at least 20% greater than its RVDB rbitscore were designated as "cellular" (i.e. corresponding to a gene from a cellular organism). Ambiguous cases without a 20% differential were designated as "viral or cellular" and proteins with no significant matches were designated as ORFans (i.e. unknown proteins).

The presence of host regions flanking the viral regions was evaluated based on the *ViralRecall* output (Table S5C). The percentages of viral regions with two, one or zero flanking regions (longer than 2 kbp) were 25.8%, 15.0% and 59.2%, respectively (i.e. 40.8% of viral regions had at least one flanking region). Of the viral regions that had two flanking regions, 89.5%, 7.0% and 3.5% had flanking regions with a total length of >200 kbp, between 20 and 200 kbp or between 2 and 20 kbp, respectively. For the viral regions that had one flanking region, the corresponding percentages were 25.3%, 36.7% and 38.0%.

## Phylogenetic analysis of viral genes

Amino acid sequences of manually-curated collections of major capsid protein (MCP) and DNA polymerase B proteins were aligned using *MAFFT* (v7.520, settings: -adjustdirectionaccurately -auto -maxiterate 1000) and phylogenetic trees were generated using *IQ-TREE* (v 2.2.2.3, settings: -m MFP -B 1000).

## Metabolic networks

Genome-scale metabolic networks were reconstructed using *AuCoMe*⁴⁶ version 0.5.1 using the *MetaCyc*¹⁸⁷ version 26 database. A first dataset, consisting of the 60 species listed in Table S1F (column "Metabolic networks") plus two public diatom genomes
already used in the initial AuCoMe study (*Fragilariopsis cylindrus* and *Fistulifera solaris*) was processed to build the largest possible database (phaeogem) for exploratory comparisons (https://gem-aureme.genouest.org/phaeogem/). Then, a second comparison was performed on all long-read species plus outgroups. Based on Multidimensional-scaling (MDS) analyses, the most divergent long-read species (*Choristocarpus tenellus*, *Laminaria digitata*, *Phaeothamnion wetherbeei* and the public genome of *Sargassum fusiforme*) were excluded to construct a 16 species dataset, balancing assembly quality and phylogenetic coverage (https://gem-aureme.genouest.org/16bestgem/). MDS plots were build using the vegan package, version 2.6-4 (https://github.com/vegandevs/vegan) with R 4.1.2,¹⁶² using Jaccard distances. A third stricter dataset (fwgem), enriched in high-quality long-read Ectocarpales, was built to address questions related to freshwater adaptation (https://gem-aureme.genouest.org/fwgem/). A set of reactions that were overrepresented in brown algae compared to the outgroup was created by taking reactions present in 100% of brown algae and less than 70% of outgroups. Reactions corresponding to genes lost in freshwater species were also extracted. These reaction sets were extracted from all the networks using the Aucomana library (https://github.com/PaulineGHG/aucomana). Online wikis (phaeogem, 16bestgem and fwgem) were generated using AuReMe.¹⁸⁸

## CAZymes

CAZyme genes were identified based on shared homology with biochemically characterised proteins, either individually or as hidden Markov model (HMM) profiles. For phylogenetic analyses, proteins were aligned using MAFFT¹³⁰ with the iterative refinement method and the scoring matrix Blosum62. The alignments were manually refined and trees were constructed using the maximum likelihood approach. Alignment reliability was tested by a bootstrap analysis using 100 resamplings of the dataset. Only bootstrap values above 60% are shown. The phylogenetic trees were displayed with MEGA.¹³¹ The annotated genes are listed in Table S4B with accession numbers.

## Sulfatases

The sulfatases encoded by each brown algal genome were identified and assigned to their respective family and subfamily using the SulfAtlas database¹⁷³,¹⁷⁴ (https://sulfatlas.sb-roscoff.fr/). Each predicted proteome was first submitted to the SulfAtlas HHM server (https://sulfatlas.sb-roscoff.fr/sulfatashmm/), which allows rapid identification of sulfatase candidates and (sub)family assignment using hidden Markov model profiles for each SulfAtlas (sub)family. Each sulfatase candidate sequence was then used as a query in a BLASTp⁸⁴ search against the SulfAtlas database (https://blast.sb-roscoff.fr/sulfatlas/). Sequences with at least 50% identity with sulfatases from marine bacteria or other marine microorganisms were considered to be contaminants. Below this threshold, additional examination of the predicted gene structure and genomic context of the candidate sequence was undertaken to identify possible horizontal gene transfers.

## Haloperoxidases

vHPO genes were identified based on sequence homology and active site conservation. Maximum likelihood phylogenetic analyses were carried out using the NGphylogeny platform at https://ngphylogeny.fr/. MAFFT was used to align vHPO sequences and alignments were automatically curated with TrimAl,¹³³ leading to the selection of 444 informative positions from the initial 1450 positions for the algal-type vHPOs and 402 informative positions from the initial 1078 positions for the bacterial-type vHPOs. Maximum likelihood trees were constructed using FastTree with the WAG+G gene model and 1000 bootstrap replicates. Maximum likelihood Newick files were formatted as circular representations using iTOL. Only bootstrap values between 0.7 and 1 were conserved. The lists of annotated vHPO genes are in Tables S4C and S4D.

## Ion channels

A search was carried out for 12 classes of ion channel in the predicted proteomes of the 21 Phaeoexplorer reference genomes plus those of two diatoms, *Phaeodactylum tricornutum* and *Thalassiosira pseudonana*. Predicted proteomes were screened using BLASTp⁸⁴ and query sequences from *Ectocarpus* species 7 and seven other species from diverse eukaryotic taxa.

## Membrane-localised proteins

Membrane protein family genes were identified either by carrying out BLASTp⁸⁴ searches of the predicted Phaeoexplorer proteomes using *Ectocarpus* species 7 sequences as queries or by recovering orthogroups containing the relevant *Ectocarpus* species 7 sequences as members. The BLASTp approach was used for DEK1-like calpains, fasciclims, tetraspanins, CHASE, ethylene-binding-domain-like and MASE1 domain histidine kinases whereas the orthogroup approach was used to recover other members of the histidine kinase family. Both approaches were used to search for integrins and transmembrane receptor kinases. For integrins the two methods detected exactly the same set of proteins. For receptor kinases the BLASTp and orthogroup analyses detected 99.3% and 98.3% of the 269 genes, respectively. For these analyses, either the whole genome dataset was analysed or only the set of 21 reference genomes (Table S1F), depending on the size of the gene family.

Manually-curated histidine kinase protein families were aligned with Muscle¹¹⁵ before phylogenetic tree construction using IQ-TREE 2¹⁶³ (version 2.3.4) with automatic model selection and 1000 bootstraps.
## Transcription-associated proteins

TAPscan v4¹³⁵ was used to analyse the transcription-associated protein (TAP) complements of 21 species. TAPscan¹³⁴ is a comprehensive tool for annotating TAPs based on the detection of highly conserved protein domains using HMM profiles with specific thresholds and coverage cut-offs. Following detection, specialised rules are applied to assign protein sequences to TAP families based on the detected domains. TAPscan v4 can assign proteins to 138 different TAP (sub)families with high accuracy.

## EsV-1-7 domain proteins

EsV-1-7 domain proteins were identified in the 31 brown algal and sister taxa genomes (Table S1F) by recovering the members of all orthogroups (with the exception of OG0000001, which is a very large OG that consisting principally of transposon sequences) that either contained one or more of a curated set of 101 EsV-1-7 domain proteins⁶² for *Ectocarpus* species 7 or contained an EsV-1-7 domain protein based on a match to the Pfam EsV-1-7 motif PF19114. The recovered proteins were screened manually for the presence of at least one EsV-1-7 domain and a total of 2018 were finally identified as members of the EsV-1-7 family.

To identify orthologues of the EsV-1-7 protein IMMEDIATE UPRIGHT⁶² (IMM), BLASTp searches of 25 brown algal and four sister taxa predicted proteomes were carried out with the amino-terminal domain of the IMM protein minus the five EsV-1-7 repeats as this domain is unique to IMM. Proteins were retained as IMM orthologues if they were more similar to IMM than to the most closely-related protein in *Ectocarpus* species 7, Ec-17_002150.

## Histones

Histone protein sequences were analysed in *Ascophyllum nodosum*, *Chordaria linearis*, *Chrysoparadoxa australica*, *Desmarestia herbacea*, *Dictyota dichotoma*, *Discosporangium mesarthrocarpum*, *Ectocarpus crouaniorum*, *Ectocarpus fasciculatus*, *Ectocarpus siliculosus*, *Fucus serratus*, *Heterosigma akashiwo*, *Pleurocladia lacustris*, *Porterinema fluviatile*, *Pylaiella littoralis*, *Saccharina latissima*, *Sargassum fusiform*, *Schizocladia ischiensis*, *Scytosiphon promiscuus*, *Sphacelaria rigidula*, *Tribonema minus* and *Undaria pinnatifida* using BLASTp against the complete predicted proteomes (https://blast.sb-roscoff.fr/phaeoexplorer/) with the histone protein sequences from the diatom *Phaeodactylum tricornutum* as queries. The genes and transcripts coding for the identified histones were then retrieved from the genomes and predicted transcripts using BLAST (https://blast.sb-roscoff.fr/phaeoexplorer/). The proteins encoded by the identified genes and transcripts were predicted with the Expasy web translator (https://web.expasy.org/translate/). In order to identify truncated proteins or incorrect start codons, the following constraints were applied: H2A proteins must start with the SGKGKGGR sequence, H2B with AKTP, canonical H3.1 and variants H3.3 with ARTKQT and H4 with SGRGKGGKGLGKGG. For the linker histone H1, protein sequences had to be lysine-rich and sequences with incorrect start codons were determined by alignments of all identified H1 proteins. For proteins with incorrect start codons, the region upstream of the correct start codon was removed. For truncated proteins, i.e. proteins whose transcripts lacked either the start (no methionine) or stop codons, the protein sequence was completed based on alignment with the corresponding genomic region using the Geneious 11.0.5 software. When the sequence could not be completed, a BLAST was performed against the Phaeoexplorer *de novo* transcriptomes (https://blast.sb-roscoff.fr/phaeoexplorer_denovo/) when this data was available (this was not possible for the public genomes *T. minus*, *U. pinnatifida* and *S. fusiforme*). Based on the nomenclature established by,¹⁸⁹ H3 histones were classified as follows: canonical H3.1 proteins harbour AT residues at positions 31-32 while histone variants H3.3 harbour TA residues, H3 proteins with other residues at positions 31-32 were named H3.4 and so on. CenH3 variants of H3 were identified by analysis with Panther 17.0 (www.pantherdb.org/tools/sequenceSearchForm.jsp?) and/or Interproscan¹³⁷ 94.0 (www.ebi.ac.uk/interpro/search/sequence/).

Species abbreviations used in histone phylogenetic trees are: Atr, *Amborella trichopoda*; At, *Arabidopsis thaliana*; Ce, *Caenorhabditis elegans*; Di, *Dictyostellum discoideum*; Dr, *Danio rerio*; Dm, *Drosophila melanogaster*; Hs, *Homo sapiens*; Mm, *Mus musculus*; Pp, *Physarum polycephalum*; Ppa *Physcomitrium patens*; Sc, *Saccharomyces cerevisiae*; Tm, *Tetrahymena thermophila*; Zm, *Zea mays*; Mp, *Marchantia polymorpha* subsp. *Ruderalis*; Bd, *Brachypodium distachyon*; Ccr, *Chondrus crispus*; Gs, *Galdieria sulphuraria*; Cm, *Cyanidioschyzon merolae*; Cr, *Chlamydomonas reinhardtii*; Ol, *Ostreococcus lucimarinus*; Ot, *Ostreococcus tauri*; To, *Thalassiosira oceanica*; Pt, *Phaeodactylum tricornutum*; An, *Ascophyllum nodosum*; Cl, *Chordaria linearis*; Ca, *Chrysoparadoxa australica*; Dh, *Desmarestia herbacea*; Ddi, *Dictyota dichotoma*; Dme, *Discosporangium mesarthrocarpum*; Ec, *Ectocarpus crouaniorum*; Ef, *Ectocarpus fasciculatus*; Es, *Ectocarpus siliculosus*; Fse, *Fucus serratus*; Ha, *Heterosigma akashiwo*; Pla, *Pleurocladia lacustris*; Pf, *Porterinema fluviatile*; Pli, *Pylaiella littoralis*; Sl, *Saccharina latissima*; Sf, *Sargassum fusiform*; Si, *Schizocladia ischiensis*; Sp, *Scytosiphon promiscuus*; Sri, *Sphacelaria rigidula*; Tm, *Tribonema minus*; Up, *Undaria pinnatifida*.

## DNA methyltransferases

Searches were carried out for methyltransferases and demethylases in the predicted proteomes of 20 of the high quality brown algal reference genome assemblies (based on Nanopore long-read sequence) plus the sister taxa *Chrysoparadoxa australica* and *Schizocladia ischiensis* using BLASTp (Table S1F). A methyltransferase reference database was constructed by recovering sequences from NCBI, ENSEMBL and UniProtKB. Methyltransferase sequences were recovered for stramenopiles such as *Nannochloropsis gaditana*, the diatom *Phaeodactylum tricornutum*, the oomycete *Phytophthora infestans* and for species from more distant lineages including *Arabidopsis thaliana*, *Homo sapiens* and the fungus *Neurospora crassa*. The proteomes of the selected brown algal strains were then queried against this database using BLASTp and matches with an e-value of < 0.001, a bit score > 70, a maximum gap of 5 and percentage identity of >30% were retained. The retained matches were screened against the NCBI, UniProt and
SwissProt databases to identify and remove contaminating bacterial or viral proteins. Methyltransferase domains were detected in the retained matches using the Simple Modular Architecture Research Tool (SMART)¹⁷⁷ domain architecture analysis and InterPro searches (https://www.ebi.ac.uk/interpro/). Sequences with methyltransferase domains were retained for further analysis. Validated brown algal methyltransferases were aligned with reference methyltransferases using ClustalW¹³⁸ 2.1.

## Spliceosome

Components of the Major Spliceosome were identified using a reference set of 147 human components (https://www.genenames.org/data/genegroup/#!/group/1518), excluding the five small nuclear RNAs (snRNAs). Including isoforms, this query set consisted of 626 proteins. These proteins were used to screen the predicted proteomes of 54 genomes (Table S1F) using BLASTp and matches were retained if they had an e-value of at most $1e^{-40}$ and coverage >30%. Searches were also carried out for components of LSm and Sm complexes which have roles as scaffolds in the formation of ribonucleoprotein particles (RNPs), in the maturation of mRNAs (including splicing, such as the cytoplasmic complex LSm1-7, LSm2-8 which is part of the core U6 snRNP and other complexes important for the formation of the 3' ends of histone transcripts), in the assembly of P-Bodies and in the maintenance of telomeres.

## Flagella proteins

A previous proteomic analysis of anterior and posterior flagella of the brown alga *Colpomenia bullosa* identified a total of 592 proteins across the two proteomes.⁴¹ Here the *Ectocarpus* species 7 orthologues of 70 of these proteins that had been detected with a very high level of confidence were used to identify the corresponding orthogroups and the presence or absence of these orthogroups was scored for seven representative species (Table S1F).

## Detection of *Porterinema fluviatile* genes differentially expressed in freshwater and seawater

Six independent cultures of *Porterinema fluviatile* were cultivated for four weeks in 140 mm Petri dishes with Provasoli-enriched culture medium,¹⁹⁰ which was renewed every two weeks. For three Petri dishes, the culture medium was based on autoclaved natural seawater (high salinity treatment), for the other three Petri Dishes natural seawater was diluted 1:19 vol/vol with distilled water (low salinity treatment). Cultures were harvested with 40 µm nylon sieves, dried with a paper towel, and immediately frozen in liquid nitrogen. RNA extraction library construction and sequencing were carried out as described in section "RNA extraction, Illumina RNA-seq library preparation and sequencing". RNA-seq reads were cleaned with Trimmomatic⁹⁹ V0.38 and then mapped to the *P. fluviatile* genome using Kallisto¹⁴⁰ version 0.44.0. Differentially expressed genes were identified using the DESeq2 package¹⁴¹ included in Bioconductor version 3.11, considering genes with an adjusted *p* < 0.05 and a log₂ fold-change > 1 as differentially expressed. To compare the differentially expressed genes in *P. fluviatile* with an equivalent set previously identified for *Ectocarpus subulatus* in a microarray experiment using nearly identical growth conditions,¹⁹¹ orthologues in the two species were detected using Orthofinder version 2.3.3. Of the 10,066 shared orthogroups, 6,606 had microarray expression data for *E. subulatus*. This information was used to classify differentially expressed genes for the two species as either shared orthologues or as lineage-specific.

## Identification of genes with generation-biased expression patterns

RNA-seq data (two to five replicates per condition) was recovered for gametophyte and sporophyte generations of ten species (Table S1F). Data quality was assessed with FastQC¹⁴² version 0.11.9 and sequences were then trimmed with Trim Galore version 0.6.5 with the parameters -length 50, - quality 24, -stringency 6, -max_n 3. The cleaned reads were mapped onto the corresponding genome for each species using HISAT2 version 2.1.0 with default options. Counting was carried out with featureCounts¹⁴⁵ from the subread package (version 2.0.1) on CDS features grouped by Parent. Transcript Per Kilobase Million (TPM) tables were generated for all conditions and differentially expressed genes were detected using DESeq2¹⁴¹ version 1.30.1. Genes were classified into six categories based on the differential expression analysis and the TPM values: gametophyte-biased, mean TPM ≥ 1 in gametophyte and sporophyte, log₂(fold change) ≥ 1, adjusted *p*-value <0.05; sporophyte-biased: mean TPM ≥ 1 in gametophyte and sporophyte, log₂(fold change) ≤ -1, adjusted *p*-value <0.05; gametophyte-specific, mean TPM <1 in sporophyte and ≥ 1 in gametophyte, log₂(fold change) ≥ 1, adjusted *p*-value <0.05; sporophyte-specific, mean TPM <1 in sporophyte and ≥ 1 in gametophyte, log₂(-fold change) ≤ -1, adjusted *p*-value <0.05; unbiased genes: mean gametophyte and sporophyte TPMs ≥ 1, log₂(fold change) <1 or >-1 and/or adjusted *p*-value ≥ 0.05; unexpressed genes, mean gametophyte and sporophyte TPM <1.

## Life cycle and thallus architecture

### Genome dataset and traits

To study the impact of body architecture, the brown algae were divided into three categories: 22 filamentous species, eight simple parenchymatous species and 13 species with elaborate thalli (Table S1F). For the life-cycle-based assessment, the groups were: 30 haploid-diploid species and six diploid species (Table S1F). Body architecture information was available for 43 species, and life cycle information was available for 36 species; species without body plan or life cycle information were not used in subsequent analyses. Two approaches were used to estimate selection intensity across the phylogeny, (i) a model-based method, and (ii) by evaluating codon usage bias and nucleotide composition. Two evolutionary models were used, one based on architecture and the other based on life cycle. For model-based methods the phylogeny was categorised based on the above traits, and selection intensity parameters were estimated using PAML¹⁴⁶ version 4.9i. Rate estimates were obtained for non-synonymous substitutions (dN), synonymous
substitutions (dS) and omega (dN/dS) for the multiple sequence alignments of all genes within each orthogroup using the variable-ratio model of CODEML from PAML, which allows different omegas for different branch categories. The traits were assigned to the branches of the phylogeny using ancestral state estimation by stochastic mapping with the phytools R package.<sup>147,162</sup>

### Evolutionary models to study impacts of body architecture

To study variation in selection intensity as a function of body architecture, we devised a model with the following trait categories: filamentous/pseudoparenchymatous (simple cell division and organisation on a single plane), parenchymatous (cell division and organisation on multiple planes) and elaborate thallus (tissue differentiation). To ensure that at least 50% of the species in each category were used in the analysis, we selected orthogroups (OGs) that contained at least 11 members for filamentous, at least four members for parenchymatous and at least six members for elaborate thallus algae. Using this filter, 1068 OGs were obtained, on which the model based on body architecture was fitted. Selection intensity parameters [rate of non-synonymous substitution (dN), rate of synonymous substitution (dS) and omega (dN/dS)] were estimated for the three trait categories for each gene alignment. We used the Wilcoxon signed-rank test to evaluate the statistical significance of differences between the selection intensity parameters (dN, dS and dN/dS) for each category.

### Evolutionary models to study the impacts of life cycle

The impact of life cycle on molecular evolution was assessed using a model with two categories consisting of diplontic and haplodiplontic species. For this model we used 1,058 OGs that contained at least three members for diploid species and at least 15 members for haploid-diploid species. Using alignments of the gene within the OGs, we estimated the selection intensity parameters for the different categories and applied the Wilcoxon signed-rank test to assess the statistical significance of differences in selection intensity between the diploid and haploid-diploid life cycles.

### Selection of intensity parameters

Omega (dN/dS) provides an estimate of the ratio of substitutions at sites under selection compared to neutral sites, and is generally used to infer the strength of purifying selection. Omega needs to be interpreted with caution because not all synonymous sites are neutral<sup>192</sup> and also synonymous substitutions are often underestimated due to saturation of synonymous sites, which might in turn impact the omega ratios.<sup>193</sup> Omega values lower than one indicate substitutions are less frequent at sites under selection compared to neutral sites and are characteristic of highly conserved genes or genes evolving under strong purifying selection. As we used primarily low copy number genes in this study, the analysed genes were expected to evolve under strong purifying selection, with omega values much lower than one. Using omega for near neutral studies is challenging because near neutral sites are determined by effective population size, that is to say, sites under mild selection constraint in larger populations can behave as neutral sites in smaller populations. It is therefore difficult to infer the amount of mutation from relative values of omega. In order to obtain better insight into selection intensity, mutation accumulation was not only investigated using rates of synonymous (dS) and non-synonymous (dN) substitutions but also by estimating codon bias and nucleotide composition. Codon usage bias was used, in addition to omega, to infer selection intensity across species as the former reflects selection efficacy at synonymous sites.<sup>194–196</sup> We inferred codon usage bias by estimating the effective number of codons (ENC) for each species using the enc method from the VHICA package.<sup>148,182</sup> The effective number of codons (ENC) quantifies the extent of deviation of codon usage of a gene from equal usage of synonymous codons. For the standard genetic code, ENC values range from 20 (where a single codon is used per amino acid implying strong codon usage bias) to 61 (implies that all synonymous codons are equally used for each amino acid<sup>197</sup>). Low ENC indicates constrained use of codons, which potentially highlights stronger codon bias due to stronger selection at synonymous sites. As nucleotide composition can also influence codon bias, we calculated the overall GC composition, GC at the third codon position (GC3) and the theoretical expected ENC (EENC) based on GC3 using local R scripts. The lower the observed ENC (OENC, estimated from the gene sequence) relative to EENC, the stronger the influence of selection due to translation on codon usage. This was studied by estimating the difference (DENC = EENC - OENC) between the expected ENC and the observed ENC.<sup>198</sup> Positive DENC indicates a role for selection constraints on codon usage in addition to the influence of nucleotide composition. DENC values of zero or less indicate that codon bias is entirely driven by nucleotide composition. DENC values were used to study the influence of translation selection and nucleotide composition on codon usage bias.

### Assembly and analysis of organellar genomes

Plastid and mitochondrial genomes were assembled *de novo* using NOVOPlasty<sup>149</sup> v3.7 and *rbcL* and *cox1* nucleotide sequences as seeds. Assembled genomes were checked by aligning reads using Bowtie2<sup>86</sup> v2.3.5.1 and processed with SAMtools<sup>150</sup> v1.5. Annotation of protein-coding genes was performed with GeSeq<sup>151</sup> v2.03. Annotation of tRNAs, tmRNAs and rRNAs was performed with ARAGORN<sup>152</sup> v1.2.38.

Maximum-likelihood (ML) phylogenetic trees were constructed using 92 plastid genomes (11 non-brown outgroup sequences) and 89 mitochondrial genomes (seven non-brown outgroup sequences). The conserved coding-region amino acid sequences of 139 plastid genes (31,159 amino acids) and 35 mitochondrial genes (7,461 amino acids) were used to construct these phylogenetic trees. The sequence for each gene was aligned individually using MAFFT<sup>130</sup> v7 (-maxiterate 1000) and then concatenated. Alignment partitions were assigned based on genes. Each of the aligned gene sequences was trimmed with *trimAl*<sup>133</sup> v1.2 (-automated1). ML phylogenetic trees were constructed with IQ-TREE 2.<sup>163</sup> The protein substitution models in each gene partition were selected using *ModelFinder*.<sup>153</sup> Statistical support for tree branches was assessed with 1,000 replicates of ultrafast bootstrap (UFBoot2).<sup>154</sup>
<figure><img src="image_32.png" /></figure>

## Analysis of *Ectocarpus* genome synteny

Global genome synteny analysis was performed using SynMap<sup>155</sup> on the CoGe platform (https://genomevolution.org/coge/) with the following genomes: *Ectocarpus crouaniorum* male, *Ectocarpus fasciculatus* male, *Ectocarpus siliculosus* male, *Ectocarpus* species 7 male and *Ectocarpus subulatus*. SynMap identifies syntenic regions between two or more genomes using a combination of sequence similarity and collinearity algorithms. Last<sup>199</sup> was used as the BLAST algorithm and syntenic gene pairs were identified using DAGChainer<sup>156</sup> with settings "Relative Gene Order", -D = 20, -A = 5. Neighbouring syntenic blocks were merged into larger blocks. Substitution rates between the synthetic CDS pairs were calculated using CodeML,<sup>146</sup> which was also implemented in SynMap, CoGe. In detail, protein sequences were aligned using the Needleman-Wunsch algorithm implemented in nalign (https://pypi.org/project/nwalign/) and then translated back to aligned codons. CodeML was run five times for each alignment using the default parameters and the lowest dS was retained, with the upper cutoff for dS values set at 2. *Ectocarpus* genes were grouped according to their age based on the phylostratigraphic analysis and by chromosomal location based on their chromosome position in *Ectocarpus* species 7. All plots and statistical analysis were carried out in R version v.4.3.1. Local synteny analysis was based on orthologous genes as identified by Orthofinder.

## Analysis of *Ectocarpus* gene evolution

Protein sequence alignments were used to remove gaps with trimAl<sup>133</sup> and then translated back to DNA with backtranseq.<sup>200</sup> Only DNA fasta files with a minimum of 70 bp were retained (831 single-copy orthologs). PhyML trees were built with Geneious v11.1.5 (https://www.geneious.com). Maximum likelihood analysis was carried out to detect site specific, branch-site specific and branch specific positive selection as well as sites under negative selection, using PAML.<sup>201</sup>

## Phylogenetic analysis of *Ectocarpus* species

Phylogenetic analysis was carried out for 11 *Ectocarpus* species plus *Scytosiphon promiscuus* as an outgroup (Table S1F). Of the 933 single-copy orthogroups identified for these 12 species, 261 high-confidence alignments were retained for gene tree and species tree inferences following the removal of low-quality alignments using BMGE.<sup>202</sup> Bayesian inference of the phylogeny of the *Ectocarpus* species complex was performed using BEAST<sup>157</sup> v2.7. The analysis was conducted under the multi-species coalescent (MSC) model, implemented in StarBEAST3<sup>158</sup> v1.1.7. The MSC model coestimates gene trees and the species tree within a multispecies coalescent framework, enabling the assessment of incongruences among genes with respect to the species tree. To account for substitution model uncertainty, bModelTest<sup>159</sup> was employed to average over a set of substitution models for each alignment. StarBEAST3 was run under both the Yule model and the strict clock model. A total of 300,000,000 Markov Chain Monte Carlo (MCMC) generations were conducted, with tree states stored every 50,000 iterations. Posterior tree samples were combined, discarding the initial 10% burn-in, using LogCombiner v2.4.7. A maximum clade credibility tree was generated using TreeAnnotator<sup>157</sup> v2.4.7.

## *Ectocarpus* introgression analysis

To distinguish introgression from shared ancestry, D estimates (i.e. ABBA-BABA tests) were generated from 36 four-taxon combinations<sup>203</sup>: four to test the level of introgression within clade 1 (i.e. *E. subulatus*, *E. crouaniorum*, *Ectocarpus* species 1, *Ectocarpus* species 2), 20 to test the level of introgression within clade 2 (i.e. *Ectocarpus* species 6, *Ectocarpus* species 7, *Ectocarpus* species 5, *Ectocarpus* species 9, *E. siliculosus*, *Ectocarpus* species 3) and 12 to test the level of introgression between these two clades. Tests were designed using a four-taxon fixed phylogeny (((P1,P2)P3)O), where P1 and P2 are closely related species from the same clade, P3 is a more divergent species that may have experienced admixture with one or both of the (P1,P2) taxa, and an out-group (O). *E. fasciculatus* was used as the out-group taxon for all ABBA-BABA tests. Details about how P1, P2 and P3 taxa were selected for each test are given in Table S6. Previous results of species tree inference were used to inform subsequent ABBA-BABA tests and to define the (((P1,P2)P3)O) phylogenies. ABBAs are sites at which the derived allele (called B) is shared between the taxa P2 and P3, whereas P1 carries the ancestral allele (called A), as defined by the outgroup while BABAs are sites at which the derived allele is shared between P1 and P3, whereas P2 carries the ancestral allele. Under incomplete lineage sorting, conflicting ABBA and BABA patterns should occur in equal frequencies, resulting in a D statistic equal to zero. Historical gene flow between P2 and P3 causes an excess of ABBA, generating positive values of D. Historical gene flow between P1 and P3 causes an excess of BABA, generating negative values of D. Patterson's D-statistic was calculated for the concatenated alignments of the 261 orthologroups. Significance was detected using a block-jackknifing approach,<sup>203-205</sup> with a block size of 5 kbp. For the jackknife procedure, one block of adjacent sites was removed n times. A Z-score was finally obtained by dividing the value of the D statistic by the standard error over n sequences of 5 kbp. The ParimonySplits network was reconstructed for the genus *Ectocarpus* using SplitsTree 4<sup>160</sup> (version 4.14.6) with 1000 bootstrap replicates.

## QUANTIFICATION AND STATISTICAL ANALYSIS

All statistical analyses are described in detail in the relevant sections of the "method details" section and the results of statistical tests are shown in the tables and figures.
ADDITIONAL RESOURCES

The Phaeoexplorer website (https://phaeoexplorer.sb-roscoff.fr) provides access to all the annotated genome assemblies described in this study as downloadable files. The output files from the Orthofinder,¹¹³ Interproscan,¹³⁷ Hectar¹⁶¹ and eggNOG-mapper¹²⁵ analyses, together with the results of the various DIAMOND¹⁰⁷ sequence similarity analyses (see section "Analyses aimed at deducing functional characteristics of predicted proteins"), can also be downloaded. In addition, the site provides genome browser interfaces for the genomes and multiple additional tools and resources including BLAST interfaces for genomes, proteomes and de novo transcriptomes, various experimental protocols, an AskOmics genomic data query interface (PhaeoAskOmics), an RShiny-based transcriptomic aggregator for the model brown alga *Ectocarpus* species 7 strain Ec32, a link to genome-wide metabolic networks for the Phaeoexplorer species and a list of project-related publications.

Additional data and results have been deposited in the CNRS Research Data depository under the title "Data for Phaeoexplorer publication: Evolutionary genomics of the emergence of brown algae as key components of coastal ecosystems" (DOI: https://doi.org/10.57745/9U1J85). Dataset description: "The Phaeoexplorer project sequenced 60 genomes corresponding to 44 brown algal and sister species. This dataset corresponds to supplementary information relating to the initial annotation of the Phaeoexplorer genomes and multiple analyses of the genome data. The dataset includes additional results of the project, together with accompanying additional figures and tables, (Additional_results.tar.gz), presubmission (v0) versions of the Phaeoexplorer genome annotation (GFF) files (GFF_v0.tar.gz) and genome-wide predicted proteomes as fasta files (Proteomes_v0.tar.gz), de novo transcriptome assemblies for the Phaeoexplorer species (RNA-seq data assembled with Trinity or rnaSPAdes; de-novo-transcriptomes.tar.gz), RepeatMasker analyses of repeat sequences (RepeatMasker.tar.gz), alignment files used to generate a phylogenetic tree for the Phaeoexplorer species (PhylogeneticTree.tar.gz), alignments used to build a densitree specifically for *Ectocarpus* species (Microevolution_Ectocarpus.tar.gz), an Orthofinder-based analysis of shared orthologues (Orthogroups.tar.gz) together with a Dollo-logic-based analysis of orthogroup gain and loss during evolution (Dollo_analysis.tar.gz), a Phylostratigraphy analysis of brown algal genes (Phylostratigraphy.tar.gz), an analysis of protein functional domain fissions and fusions (CompositeGenes.tar.gz), Interproscan analyses of protein domains (InterProScan.tar.gz), Hectar predictions of protein subcellular localisations (Hectar.tar.gz), eggNOG output providing information about predicted protein functions (eggNOG.tar.gz), RNA-seq-based data on gene expression levels (mRNAexpression.tar.gz), results of a search for genes acquired via horizontal gene transfer (HGT.tar.gz), analyses of intron conservation across genomes (Introns_conservation.tar.gz), an analysis of tandem gene duplications (Tandemely_duplicated_genes.-tar.gz), CAFE5 reconstruction of gene family amplifications (CAFE5.tar.gz), comparisons of CDS size with the *Ectocarpus* reference genome that were used to evaluate gene model completeness (CDS_size.tar.gz), a DESeq2 analysis of differential gene expression between the sporophyte and gametophyte generations of several brown algal species (DEG_LifeCycle.tar.gz), information about orthogroups selected to analyse the effects of morphological complexity and life cycle structure on gene evolution (Genes_selection.-tar.gz). Each individual dataset contains a README file explaining its content. Detailed information about the methodology used for each analysis can be found in the STAR Methods section of the manuscript preprint (https://doi.org/10.1101/2024.02.19.579948). The majority of these analyses and datasets can also be accessed via the Phaeoexplorer website (https://phaeoexplorer.sb-roscoff.fr/)."
Supplemental figures

<figure><img src="image_33.png" /><figcaption>B</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>C</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>D</figcaption></figure>

<figure><img src="image_36.png" /></figure>
<figure><img src="image_37.png" /></figure>

Figure S1. Taxonomic diversity and assembly quality of the Phaeoexplorer genomes, and geographic localization of species and strains, related to Figure 1

(A) Taxonomic distribution and assembly quality (contig N50) of the Phaeoexplorer genome dataset (blue) and previously published brown algal genomes (brown). "Reference" quality Phaeoexplorer genomes are circled in black.

(B) Number of genes annotated in each genome (upper). BUSCO scores for the predicted proteome of each genome (middle). Correlation of coding sequence (CDS) lengths for each species with the corresponding sequences from the *Ectocarpus* species 7 reference genome (lower). Previously published genomes are indicated in gray. longreads, genomes assembled using long reads. F, female; M, male.

(C) World map indicating the positions of the sampling sites for the strains sequenced in this study.

(D) North-south distributions of the species analyzed in this study. The data shown were obtained from the Global Biodiversity Information Facility (gbif.org). Asterisks indicate the sampling latitude for the sequenced strains of each species.
<figure><img src="image_38.png" /><figcaption>Figure S2. Fossil-calibrated time trees and Dollo logic analysis of orthogroups, related to Figures 1 and 5</figcaption></figure>

(A) Fossil-calibrated maximum-likelihood phylogenetic tree of the brown algal and closely related taxa analyzed in this study. Bayesian divergence time estimation using 32 nuclear protein sequences, together with two fossil calibrations and a calibration for the root based on Choi et al.¹ (numbered blue circles). Gray bars on the nodes show 95% highest posterior density region (HPD) intervals of the node ages. The Great Ordovician Biodiversification Event (GOBE) is shown as an orange box, the Permian-Triassic mass extinction as a purple band and Pangea rifting as a blue box.

(B) Dollo-logic-based analysis of gene family gain and loss during the emergence of the brown algae. Cladogram indicating orthogroup (OG) gain and loss during the emergence of the Phaeophyceae based on Dollo parsimony analysis. Taxonomic classes, orders and families are indicated in brown (brown algae) or gray (outgroup species) on the right. The nodes of the tree (n0-n22) are numbered in red, and listed on the left with the corresponding name, if one exists. The number of OGs predicted to be present at each node is indicated by the circles, and the branches are colored according to overall gene family gain.

(C) Fossil-calibrated phylogenetic tree for 11 *Ectocarpus* species. Extracted from the fossil-calibrated tree shown in (A). Numbered gray bars on the nodes show 95% HPD intervals of the node ages.
<figure><img src="image_39.png" /><figcaption>C</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>B</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>D</figcaption></figure>

<figure><img src="image_42.png" /></figure>

1. Phaeophyceae

<figure><img src="image_43.png" /></figure>

3. Fucales

<figure><img src="image_44.png" /></figure>

2. FDI clade

<figure><img src="image_45.png" /></figure>

4. Fucales

<figure><img src="image_46.png" /></figure>

5. Laminariales

<figure><img src="image_47.png" /></figure>

6. Fucales

<figure><img src="image_48.png" /></figure>

7. Laminariales

<figure><img src="image_49.png" /></figure>

5. Ectocarpales

<figure><img src="image_50.png" /><figcaption>E</figcaption></figure>

Gene family amplifications and reductions observed at index of the phylogenetic tree of COI-ITS analysis

<figure><img src="image_51.png" /></figure>
**Figure S3. Genome-wide analyses of gene family evolution, related to Figure 2**

(A) Functional and structural features of founder genes from three taxonomic levels. FDI, Fucophycideae/Dictyotales/Ishigeales; PS, Phaeophyceae plus Schizocladiophyceae; PX, Phaeophyceae plus Xanthophyceae.

(B) Inference of HGT origins from 74 species based on monophyletic most similar homolog (MMSH) analysis. The tree on the left, which is derived from the NCBI taxonomy tree, indicates the source taxa for HGT-derived genes. The tree at the bottom of the figure, which was constructed using single-copy genes, indicates the species that have received genes by HGT. The middle part of the figure indicates the number of HGT genes transferred from each source to the receiver genomes. The panel on the right illustrates the number of HGTs from each phylum. The legend in the lower left corner provides a reference for the circle size, which corresponds to 10, 100, or 300 HGT gene counts.

(C) Composite gene analysis. Phylogenetic distribution of fused (blue), split (green), and non-remodeled (gray) gene family originations across the evolution of brown algae (middle). Pie charts on each branch of the phylogeny indicate the relative contribution of gene fusion and fission to the overall emergence of novel gene families, quantified by the area of the circle. Brown algal species are indicated in brown and other stramenopiles in black. Note that only the topology of the species tree is displayed here, without specific branch lengths. Right: bar plot indicating the percentage of gene families retained in extant genomes among all gene families that emerged during the evolution of the species set. Left: bar plot representing the distribution of gene families in COG functional categories for functionally annotated fused, split, and non-remodeled orthologous groups. The functional annotation assigned to an orthogroup corresponded to the most frequent functional category annotated for the members of each orthogroup. Asterisks next to the bars indicate statistically significant differences between remodeled and non-remodeled gene families (p < 0.05, two-sided chi² test with Yates correction).

(D) Left: gene families (orthogroups) significantly amplified (binomial test) in the brown algae compared with outgroup taxa. Orthogroups amplified in specific groups of species are indicated by green rectangles. Right: pie charts representing the proportions of manually determined functional categories for each group of amplified gene families highlighted in the left.

(E) Plot showing the timing (phylogenetic clade) of the amplification of gene families (y axis) and the timing of the appearance (gain) of the amplified gene family (x axis, based on the Dollo analysis of orthogroups). Nodes (e.g., n2) are as indicated in Figure S2B.
A

<figure><img src="image_52.png" /></figure>

B

Intron gains

<figure><img src="image_53.png" /></figure>

<figure><img src="image_54.png" /></figure>

<figure><img src="image_55.png" /></figure>

<figure><img src="image_56.png" /><figcaption>**Figure S4.** Structural features of the predicted genes in the Phaeoexplorer genomes, related to Figure 1</figcaption></figure>

(A) Genome and gene statistics for the 21 reference genomes (Figure S1A; Table S1B). Violin plots display size distributions. Cross or triangle, mean; diamond or square, median. For intergenic regions, half violins on the left correspond to intergenic distances between adjacent genes on opposite strands, and half violins on the right to intergenic distances between adjacent genes on the same strand.

(B) Intron acquisition during the emergence of the brown algal lineage (left). The analysis examined 949 introns in 235 conserved orthologs. Colors correspond to the species distributions of introns. Hatched blue indicates introns that are shared with at least two outgroups (ancestral introns). The right panel shows an example of an intron conservation profile for the orthogroup OG0004854. Color code and species numbering as for the left. F, female; M, male.
A

S1s
BROWN ALGAE

<figure><img src="image_57.png" /></figure>

GTs

<figure><img src="image_58.png" /></figure>

GHs

<figure><img src="image_59.png" /></figure>

PLs total CAZymes

<figure><img src="image_60.png" /></figure>

number of genes

B

Nb. of genes

<figure><img src="image_61.png" /></figure>

C

Tree scale: 1

<figure><img src="image_62.png" /></figure>

Limited set of algal-type vHPOs
Frequent loss of one algal-type

Multiple independent
algal-type 1 and II
expansions

Algal-type
divergence or duplication

Independent acquisition
of algal-type and bacterial-type
vHPOs

D

Tree scale: 1

<figure><img src="image_63.png" /></figure>

Tree scale: 1

<figure><img src="image_64.png" /></figure>

E

Life cycle

dN/dS ratio

<figure><img src="image_65.png" /></figure>

<figure><img src="image_66.png" /></figure>

dN

<figure><img src="image_67.png" /></figure>

<figure><img src="image_68.png" /></figure>

dS

<figure><img src="image_69.png" /></figure>

Developmental complexity

dN

<figure><img src="image_70.png" /></figure>
**Figure S5. Analyses of metabolism gene families and life cycle-related gene evolution, related to Figures 1 and 3**

(A) Counts of numbers of genes predicted to encode glycosyltransferases (GTs), glycoside hydrolases (GHs), polysaccharide lyases (PLs), and all CAZymes (GTs, GHs, PLs, carbohydrate esterases [CEs], auxiliary activities [AAs], carbohydrate-binding modules [CBMs]), showing numbers for both full-length proteins (dark colors) and fragments (light colors). The data are averaged by order with the number of species analyzed per order in brackets.

(B) Number of genes for selected CAZyme families in brown algae. Only CAZyme families with at least three members per genome on average, are shown. The species analyzed are the same as in (A). Counts include full-length proteins and fragments.

(C) vHPO genes identified in the 21 reference genomes based on sequence homology and active site conservation. vHPO genes are indicated by a green cross or a number and absence by a red cross.

(D) Maximum-likelihood phylogenetic trees for 259 algal-type vHPOs (left) and for bacterial-type vHPOs (right). Algal-type I vHPOs are colored in blue and algal-type II vHPOs in violet. The clades that have structurally or biochemically characterized enzymes are highlighted in red for vBPOs and in yellow for vIPOs. Strongly supported representative branches have been collapsed. The clade names correspond either to taxa or to individual gene names. For the bacterial vHPOs, the FastTree reconstruction tool with 1,000 bootstraps was drawn as a circular representation taking the gammaproteobacterial group as the starting point to arbitrary root the tree. Brown algal branches are colored in brown. Green dots indicate bootstrap values of between 0.7 and 1.0 (1,000 replicates).

(E) Rates of gene evolution in relation to life cycle structure and developmental complexity. Left, violin plots showing the distribution of omega (dN/dS), rates of non-synonymous substitution (dN), and rates of synonymous substitution (dS) for brown algae with haplodiplontic or diplontic life cycles. Right, violin plots showing the distribution of rates of non-synonymous substitution (dN) for brown algal species with filamentous, simple parenchymatous, or elaborate parenchymatous thalli. The *p* values are for pairwise (gene-by-gene) Wilcoxon tests and the percentage of genes that exhibited the same patterns of differences in dN/dS, dN, or dS values as the median values are indicated. Significant differences are indicated by an asterisk.
**Figure S6. Phylogenetic trees of viral and putatively virus-derived genes, related to Figure 6**

(A and B) *Nucleocytoviricota* major capsid protein and DNA polymerase phylogenies, respectively. The trees were both generated from aligned amino acid sequences with the Q.yeast+F+R6 model and 1,000 bootstraps in IQ-TREE. Subgroup labels refer to *Phaeovirus* genotypes. All sequences identified by this study are in bold and highlighted in yellow if they clustered outside the *Phaeovirus*. *Phycodnaviridae*.

(C) Phylogenetic trees for three classes of histidine kinase. Membrane-associated sensor1 domain (MASE1) class: brown algal genes are more closely related to viral than to closely related outgroup species genes, EsuBft1789.4 is located within a viral clade. Ethylene-binding-domain-like (EBD) class: viral-related clade limited to Ectocarpales. Cyclases/histidine kinases associated sensory extracellular domain (CHASE) class: complex pattern suggesting possible multiple HGTs. Genes are classed as algal (green or light green background label, strong or weak prediction) or inserted viral sequences (blue or light blue background label, strong or weak prediction) based on exon number, expression level and genomic context (neighboring monoexonic or multiexonic genes). See Table S4J for the gene name abbreviations. Branch colors signify brown algae (green), closely related outgroup taxa (violet) or EsV-1 genes (blue). TPM, maximum TPM; asterisks, gene located in an identified VR; nd, not determined.