*2.1. Identification and Phylogenetic Analysis of F. hepatica Superoxide Dismutases*

The *F. hepatica* SOD gene sequences were identified by BLAST analysis using the previously reported *F. hepatica* SOD sequence (Kim et al., 2000 [16]; AF071229) against the *F. hepatica genome* (WormBase ParaSite Version WBPS16 (WS279): PRJEB6687 and PR-JEB25283; Cwiklinski et al., 2015 [5]; FhSOD1: BN1106\_s3189B000243; maker-scaffold10x\_61\_ pilon-snap-gene-0.36, FhSOD2: snap\_masked-scaffold10x\_1664\_pilon-processed-gene-0.2, FhSOD3: BN1106\_s4478B000037; maker-scaffold10x\_713\_pilon-snap-gene-0.105). Annota-

tion of the resulting sequences was confirmed using in silico tools (Uniprot, Gene Ontology (GO), and InterProScan).

Homologous trematode SOD DNA sequences were identified and retrieved using BLAST analysis of publicly available genome databases at WormBase Parasite (http:// parasite.wormbase.org/index.html Version WBPS16 (WS280), accessed on 7 April 2021). Genomic DNA sequences were imported, manually inspected, translated, and aligned in CLC Main Workbench 21.0.3 (Table S1). Identification of N-terminal signal peptide sequences was carried out using the SignalP and TMHMM plugin (Version 21.0) in CLC Main Workbench (Version 21.0.3). Previously characterised mammalian SOD sequences (*Ovis aries, Bos taurus, Bos indicus*, and *Homo sapiens*) were downloaded from GenBank (National Center for Biotechnology Information, NCBI) (Table S1). Ambiguous regions (i.e., containing gaps and/or poorly aligned) were removed from the resultant amino acid alignment of all mammalian and trematode sequences with Gblocks (v0.91b) using the following parameters: minimum length of a block after gap cleaning: 10; positions with a gap in less than 50% of the sequences were selected in the final alignment if they were within an appropriate block; all segments with contiguous non-conserved positions longer than 8 were rejected; minimum number of sequences for a flank position: 85% [18,19]. The resultant sequence alignment spanned 107 amino acids (Gly 67-Gly173, relative to FhSOD1), containing 43 amino acid sequences, which was submitted for Smart Model Selection (SMS) using PhyML 3.0 (Figure S1) [20]. Evolutionary history was inferred using the maximum likelihood method and the Whelan and Goldman (WAG + G+I) model, with 1000 bootstrap support in MEGA11 [21,22]. Initial tree(s) for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using the JTT model, and then selecting the topology with superior log likelihood value. A discrete Gamma distribution was used to model evolutionary rate differences among sites (5 categories (+G, parameter = 1.3434)). The rate variation model allowed for some sites to be evolutionarily invariable ([+I], 11.21% sites). Amino acid sequence similarity and identity was determined using UniProt ClustalO [23].
