Ultrastructural, Cytochemical, and Comparative Genomic Evidence of Peroxisomes in Three Genera of Pathogenic Free-Living Amoebae, Including the First Morphological Data for the Presence of This Organelle in Heteroloboseans

Abstract Peroxisomes perform various metabolic processes that are primarily related to the elimination of reactive oxygen species and oxidative lipid metabolism. These organelles are present in all major eukaryotic lineages, nevertheless, information regarding the presence of peroxisomes in opportunistic parasitic protozoa is scarce and in many cases it is still unknown whether these organisms have peroxisomes at all. Here, we performed ultrastructural, cytochemical, and bioinformatic studies to investigate the presence of peroxisomes in three genera of free-living amoebae from two different taxonomic groups that are known to cause fatal infections in humans. By transmission electron microscopy, round structures with a granular content limited by a single membrane were observed in Acanthamoeba castellanii, Acanthamoeba griffini, Acanthamoeba polyphaga, Acanthamoeba royreba, Balamuthia mandrillaris (Amoebozoa), and Naegleria fowleri (Heterolobosea). Further confirmation for the presence of peroxisomes was obtained by treating trophozoites in situ with diaminobenzidine and hydrogen peroxide, which showed positive reaction products for the presence of catalase. We then performed comparative genomic analyses to identify predicted peroxin homologues in these organisms. Our results demonstrate that a complete set of peroxins—which are essential for peroxisome biogenesis, proliferation, and protein import—are present in all of these amoebae. Likewise, our in silico analyses allowed us to identify a complete set of peroxins in Naegleria lovaniensis and three novel peroxin homologues in Naegleria gruberi. Thus, our results indicate that peroxisomes are present in these three genera of free-living amoebae and that they have a similar peroxin complement despite belonging to different evolutionary lineages.


Introduction
Free-living amoebae are widely distributed in nature. These cosmopolitan protozoa have been found in soil (Reyes-Batlle et al. 2014) and aquatic environments such as ponds and swimming pools (Aghajani et al. 2016). Of the many freeliving amoeba described to date, three genera from two different taxonomic groups have gained medical importance, as they have been found to be the etiological agents of various human infections. Naegleria fowleri (Excavates: Discoba; Heterolobosea; Vahlkampfiidae) (Adl et al. 2019) can produce a highly lethal infection of the central nervous system (CNS) known as primary amoebic meningoencephalitis (Carter 1972). Balamuthia mandrillaris (Amoebozoa: Discosea; Centramoebia; Acanthopodida) (Adl et al. 2019) is responsible for the so-called balamuthiasis, a lethal lesion of the CNS that occurs in various animals and can occasionally be observed in humans (Visvesvara et al. 1993). In addition, different species of Acanthamoeba (Amoebozoa: Discosea; Centramoebia; Acanthopodida) (Adl et al. 2019) such as A. castellanii, A. polyphaga, A. culbertsoni, and A. palestinensis can produce a subacute or chronic CNS infection known as granulomatous amoebic encephalitis (Mart ınez 1991;Marciano-Cabral et al. 2000). Acanthamoeba castellanii can also produce acanthamoebic keratitis, a sightthreatening infection of the cornea which can occur both in immunocompromised and healthy individuals (Niederkorn et al. 1999;Sharma et al. 2004). Other Acanthamoeba species such as A. royreba (Gonz alez-Robles et al. 2013;Omaña-Molina et al. 2016), A. polyphaga (Dini et al. 2000;da Rocha-Azevedo and Costa e Silva-Filho 2007), and A. griffini (Ledee et al. 1996;Gonz alez-Robles et al. 2014) have also been isolated from keratitis cases. The leading risk factor for acanthamoebic keratitis is the use of contact lens, as more than 85% of the patients diagnosed with this infection are contact lens users (Cheung et al. 2016;Padhi et al. 2017).
Peroxisomes, also referred to as microbodies, are small organelles that are present in almost all eukaryotic cells. The basic function of peroxisomes is the b-oxidation of fatty acids and the elimination of reactive oxygen species; however, there is a considerable variation in the enzymatic content of these organelles in different organisms. For example, the enzymes of fungal peroxisomes are involved in the synthesis of antibiotics, whereas the peroxisomes of trypanosomatids contain enzymes that perform various glycolytic reactions (Pieuchot and Jedd 2012).
Although the functions and enzymatic content of peroxisomes can vary greatly, all peroxisomes share the same biogenesis route. The different processes involved in the biogenesis of peroxisomes, such as the insertion of proteins into the membrane of the organelle and the import of matrix proteins, are carried out by a group of proteins known as peroxins (Distel et al. 1996). To date, 36 different peroxins have been identified in protists, fungi, plants, and mammals (Galland and Michels 2010;Islinger et al. 2018).
Morphologically, peroxisomes are usually ovoid to spherical, with a size ranging from 0.1 to 1 mm or even 1.5 mm in diameter. These organelles are delimited by a single membrane that encloses a fine fibrogranular matrix which contains numerous enzymes (Smith and Aitchison 2013). The high protein concentration in the peroxisomal matrix can sometimes result in the formation of crystalline inclusions (Sibirny 2016). Biochemically, peroxisomes were first identified by the association of various H 2 O 2 -producing oxidases with catalase, an enzyme that degrades hydrogen peroxide into water and molecular oxygen (De Duve and Baudhuin 1966). The alkaline 3,3 0 -diaminobenzidine (DAB) reaction for catalase is considered as a specific cytochemical staining for the presence of peroxisomes, and this technique has routinely been used to detect peroxisomes by light and electron microscopy (Schrader and Fahimi 2008).
Data regarding the presence of peroxisomes in free-living amoebae are variable. Mü ller and MØller (1969) and Childs (1973a) performed biochemical assays that demonstrated the presence of peroxisomes in Acanthamoeba and Hartmanella culbertsoni, but there is no information regarding the existence of this organelle in B. mandrillaris. As for N. fowleri, there is no biochemical or morphological evidence for the presence of peroxisomes in this Heterolobosean. However, a couple of studies identified several proteins with peroxisomal targeting signals as well as various putative peroxin homologues (Pex1,3,4,5,6,7,10,12,and 19) in Naegleria gruberi, a nonpathogenic relative of N. fowleri (Fritz-Laylin et al. 2010;Opperdoes et al. 2011). Schlü ter et al. (2006 proposed that four of these peroxins (Pex3, 10, 12, and 19) could be considered as unequivocal markers for the in silico identification of peroxisomes in eukaryotic cells after observing that they are present in all organisms with peroxisomes and absent in those that lack these organelles. Nonetheless, recent studies in various protozoan lineages have found that at least one of these markers (Pex3) is absent in organisms that are predicted to have peroxisomes (Opperdoes et al. 2016;Ludewig-Klingner et al. 2018;Paight et al. 2019), which suggests that this peroxin may not be an unequivocal indicator for the presence of peroxisomes.
The present study investigated the existence of peroxisomes in A. castellanii, A. griffini, A. polyphaga, A. royreba, and B. mandrillaris (Amoebozoa) as well as in N. fowleri and Naegleria lovaniensis (Heterolobosea). Our morphological and cytochemical observations demonstrate that peroxisomes are present in these organisms. These experimental results are further supported by comparative genomic analyses, which allowed us to identify a complete set of predicted peroxins that are necessary for peroxisome biogenesis and maintenance in these Amoebozoans and Heteroloboseans. Additionally, we identified several putative catalase homologues with predicted peroxisomal targeting signals in all of these amoebae, further supporting our cytochemical results. Thus, although they are phylogenetically distant, these freeliving amoebae contain peroxisomes and a complete set of peroxins.

Materials and Methods
Amoebae Naegleria fowleri and four Acanthamoeba species (A. castellanii, A. griffini, A. polyphaga and A. royreba) were grown in borosilicate tubes (Pyrex, Mexico) at 37 and 30 C, respectively, and maintained in axenic culture in 2% Bacto Casitone culture medium (pancreatic digest of casein, Becton Dickinson, Sparks, MD) supplemented with 10% fetal bovine serum (JR Scientific Inc). Trophozoites were harvested during the logarithmic phase of growth (72 h) by chilling the culture tubes in an ice-water bath for 5 min and they were pelleted by centrifugation at 280 Â g for 5 min.
Axenic Balamuthia mandrillaris trophozoites from the ITSON-1 strain (Lares-Jim enez et al. 2014) were grown in 25 cm 2 cell culture flasks (TPP, Switzerland) in 4.5 ml of BM-3 culture medium (Schuster and Visvesvara 1996) supplemented with 10% fetal serum. Cultures were incubated at 37 C in a 5% CO 2 atmosphere until they reached confluence. Afterwards, the culture medium was carefully removed and the amoebae were washed with phosphate buffered saline pH 7.2 which was previously heated to incubation temperature.

Cytochemical Staining
To detect catalase activity, we used a modified technique proposed by Beard and Novikoff (1969). Briefly, samples were fixed with 2.5% glutaraldehyde in 0.1 M sodium cacodylate buffer pH 7.2 for 1 h at room temperature. After fixation, the cells were washed twice with 0.1 M cacodylate buffer and incubated for 60 min at 37 C in an incubation medium containing 0.05 M propanediol buffer pH 10 (Merck-Schuchardt) supplemented with 20 mg/10 ml of 3,3 0 -diaminobenzidine tetrahydrochloride (Sigma-Aldrich) and 0.2% H 2 O 2 . This solution was filtered and the pH was adjusted to 9 prior to incubation. Mouse liver samples were used as a positive control to detect peroxisome activity. To verify the specificity of the reaction, mouse liver samples were subjected to the same procedure without the addition of 0.2% H 2 O 2 .

Electron Microscopy
Following incubation, cells were washed twice with 0.1 M cacodylate buffer and postfixed for 1 h with 1% osmium tetroxide in cacodylate buffer. Trophozoites were pelleted by centrifugation, dehydrated, and embedded in polybed resins. Afterwards, samples were polymerized at 60 C for 24 h. Thin sections (70 nm) were stained with uranyl acetate and lead citrate and examined in a JEOL JEM-1011 transmission electron microscope.
The protein sequences of the putative A. castellanii and N. fowleri peroxin homologues were retrieved directly from AmoebaDB. For A. polyphaga, A. royreba, B. mandrillaris and N. lovaniensis, the GeneWise algorithm (Birney et al. 2004) from the European Bioinformatics Institute (https://www.ebi. ac.uk/Tools/psa/genewise/; last accessed July 1, 2020) was used to obtain the protein sequence for each predicted peroxin. Using this algorithm, we compared the sequence of the query protein with the scaffold where the TBlastN analysis had identified the presence of the putative homologue. For sequences predicted to be located on the reverse strand, the reverse complement sequence of the scaffold was obtained before using it as input for the comparison.
EMBOSS Needle (https://www.ebi.ac.uk/Tools/psa/emboss_ needle/; last accessed July 1, 2020) (Needleman and Wunsch 1970) was used to determine the amino acid identity between the different members of the Pex11 family found in the Amoebozoans and Heteroloboseans included in this study (supplementary text 1, Supplementary Material online).

Assessment of Genomic Completeness
For the phylogenetic analyses, a data set of protein sequences from various protozoans was assembled by performing BlastP searches against the NCBI and the EuPathDB (https://eupathdb. org/eupathdb/; release 46; last accessed July 1, 2020) (Aurrecoechea et al. 2017) databases, using the A. castellanii and N. gruberi Pex1 and Pex5 sequences as queries for the analyses. To avoid losing distant homologues, sequences from protozoans with E-values 0.05 were included in the data set. The sequences were then analyzed with NCBI's Conserved Domain Search to verify that they contained the expected domains. Sequences were aligned with ProbCons 1.12 (Do et al. 2005) (http://www.phylogeny.fr/one_task.cgi?task_ty-pe¼probcons; last accessed July 1, 2020) using the default parameters. Alignments were trimmed with BMGE 1.12.1 (Criscuolo and Gribaldo 2010), which is available on the NGPhylogeny platform (https://ngphylogeny.fr/; last accessed July 1, 2020) (Lemoine et al. 2019) using the default parameters (BLOSUM62 matrix; minimum block size of 5). RAxML-HPC2 8.2.12 (Stamatakis 2014) was used to perform maximum likelihood analyses using the protein GAMMA model and the LG amino acid substitution matrix (Le and Gascuel 2008). Bootstrap analyses (1,000 replicates) were performed to estimate branch support. These analyses were run using the CIPRES Science Gateway v3.3 (https://www.phylo.org/portal2/login!input.action; last accessed July 1, 2020) (Miller et al. 2010). The resulting trees were visualized using iTOL 5.3 (https://itol.embl.de/; last accessed July 1, 2020) (Letunic and Bork 2019) and rooted midpoint.

Results
Of the numerous species of free-living amoebae described to date, only Naegleria fowleri (Heterolobosea),

Balamuthia mandrillaris, Sappinia pedata and several
Acanthamoeba species (Amoebozoa) are known to cause fatal infections of the CNS in humans (Visvesvara et al. 2007;Visvesvara 2013). In most of these amoebae, peroxisomes have not been described, possibly because of their variable size (0.1-1.0 lm) and scarce number. Thus, to determine if peroxisomes are present in three of these four genera of pathogenic free-living amoebae, we first performed morphological and cytochemical observations of trophozoites by transmission electron microscopy. Round structures limited by a single membrane and varying in size from 0.3 to 0.6 mm were observed in thin sections of A. castellanii, B. mandrillaris, and N. fowleri ( fig. 1D-F). These structures had a dark granular content similar to the one found in mouse liver control samples, which were used as positive controls for the identification of peroxisomes Fig. 1-Identification of peroxisomes in three genera of pathogenic free-living amoebae by transmission electron microscopy. Spherical structures with a dark granular content and limited by a single membrane were observed in Acanthamoeba griffini (A), Acanthamoeba polyphaga (B), Acanthamoeba royreba (C), Acanthamoeba castellanii (D), Balamuthia mandrillaris (E), and Naegleria fowleri (F). These structures ranged in size from 0.2 to 0.9 mm and they were morphologically similar to peroxisomes from mouse liver samples, which are shown for comparison (G, H). M, mitochondrion; Nu, nucleus. Bar ¼ 0.5 mm.
Next, trophozoites were treated with alkaline DAB and hydrogen peroxide to detect catalase activity. Examination of thin sections by transmission electron microscopy revealed the presence of distinct round particles measuring between 0.12 and 0.4 mm which showed positive reaction products for this peroxisome-specific cytochemical staining (Baumgart et al. 2003). These uniformly stained particles were clearly observed in A. castellanii, B. mandrillaris, and N. fowleri, but the number of peroxisomes per cell was low ( fig. 2D-F). In other Acanthamoeba species incubated with DAB and H 2 O 2 (A. griffini, A. polyphaga, and A. royreba), these structures were also clearly observed, varying only in size (0.1-0.4 mm) and in electron density ( fig. 2A-C). Importantly, labeling was dependent on the presence of hydrogen peroxide, as no electrondense deposits were observed when control samples were incubated without this reagent ( fig. 2H). The positive reaction products from this cytochemical staining strongly suggest that peroxisomes are present in N. fowleri, B. mandrillaris, and in various Acanthamoeba species.
Several authors have suggested that the identification of a number of peroxins in the genome of a given organism can be indicative of the presence of peroxisomes, as these proteins are present in all peroxisomes despite their metabolic content (Gabald on 2010(Gabald on , 2018. Therefore, to support our previous morphological and cytochemical observations, we performed systematic comparative genomic analyses to determine if orthologues of different peroxin proteins were present in Fig. 2-Peroxisome identification in three genera of pathogenic free-living amoebae by cytochemical staining. Thin sections of trophozoites were treated with diaminobenzidine and hydrogen peroxide to detect catalase activity. Positive reaction products were clearly seen in Acanthamoeba griffini (A), Acanthamoeba polyphaga (B), Acanthamoeba royreba (C), Acanthamoeba castellanii (D), Balamuthia mandrillaris (E) and Naegleria fowleri (F). In all of these amoebae, the labeling was deposited in round structures with a uniform electrondense content. Mouse liver samples were used as a positive control for the cytochemical staining (G). As a negative control for the cytochemical reaction, a mouse liver sample was incubated without H 2 O 2 (H). ER, endoplasmic reticulum; M, mitochondrion; Nu, nucleus. Bar ¼ 0.5 mm.
As shown in table 1 and in the supplementary tables 2-8, Supplementary Material online, a complete set of peroxins is predicted to be present in the genomes of these free-living amoebae.
The sequences of the 15 putative peroxin homologues found in A. castellanii were identified by using the corresponding H. sapiens and N. gruberi protein sequences as queries for BlastP analyses against the A. castellanii protein database stored in AmoebaDB (table 1 and supplementary  table 2, Supplementary Material online). All of these putative peroxin sequences were identified with the human queries (except for Pex4 and Pex22, which are not present in humans) and with the N. gruberi queries (except for Pex2, Pex13, Pex14, and Pex16, which have not been described in this amoeba). We also identified three distinct members of the Pex11 family (labeled Pex11a, Pex11b, and Pex11c) using the human queries and the presence of two of them was confirmed with the N. gruberi query. At the amino acid level, Pex11a was 20% identical to Pex11b but only 3.5% identical to Pex11c, whereas Pex11b and Pex11c shared a 4% identity at the protein level (supplementary text 1, Supplementary Material online). Interestingly, we also found a potential Pex26 homologue in A. castellanii. Although the BlastP and the reciprocal BlastP E-values were not significant, we did find the presence of the Pex26 domain in this sequence (supplementary table 2, Supplementary Material online).
The draft genomes of A. polyphaga and A. royreba are available on the NCBI WGS database as part of the Phylogenomics of Acanthamoeba species project (accession PRJEB7687). There are two draft genomes for A. polyphaga (CDFK00000000.1 and LQHA00000000.1) and one draft genome for A. royreba (CDEZ00000000.1). By performing TBlastN searches using the protein sequences of the peroxins identified in A. castellanii as queries, we identified 15 putative peroxin homologues in these amoebae (table 1 and supplementary tables 3

and 4, Supplementary Material online).
Similarly to A. castellanii, we found three distinct members of the Pex11 family (also labeled as Pex11a, Pex11b, and Pex11c) in these organisms. For A. polyphaga, Pex11a is 20% identical to Pex11b and 5% identical to Pex11c, whereas Pex11b and Pex11c share a 6% identity at the amino acid level (supplementary text 1, Supplementary Material online). For A. royreba, the amino acid identity between Pex11a, and Pex11b, and Pex11c is 19% and 23%, respectively, whereas the predicted Pex11b and Pex11c protein sequences are 20% identical (supplementary text 1, Supplementary Material online). In contrast to A. castellanii, we found more than one putative homologue for most peroxins in both of these Acanthamoeba species. For A. polyphaga, we found two or more putative peroxin homologues in either the CDFK or the LQHA assemblies, with the exception of Pex2, Pex10, Pex13, and Pex14, which only have one predicted homologue in each assembly (supplementary table 3, Supplementary Material online). All of these putative homologues have the expected functional domains except for the Pex12 sequences, which have either the mRING_PEX12 domain or the PEX2_PEX12 domain. To determine if this could be due to a high level of fragmentation of the A. polyphaga assemblies, we examined the putative Pex12 sequences and we found that one of them was located in a very short contig (2,250 nt) and the rest were located in the terminal part of their respective contigs, which explains why we were unable to find a putative Pex12 sequence with both domains. Importantly, when the sequences bearing the different PEX12 domains are aligned, they appear to be part of the same sequence (supplementary text 2, Supplementary Material online).
Similarly, we also found more than one putative homologue for each peroxin in A. royreba, except for Pex2, Pex7, Pex19, and Pex22 (supplementary table 4, Supplementary Material online). Additionally, we found various potential Pex26 homologues in these organisms (supplementary tables 3 and 4, Supplementary Material online), but since we could only identify the corresponding domain in one of the sequences from A. polyphaga (supplementary table 3, Supplementary Material online), further studies are necessary to confirm that these sequences are indeed Pex26 homologues. Overall, our bioinformatic results, together with the previous ultrastructural and cytochemical observations made in this study, indicate that peroxisomes are also present in A. polyphaga and A. royreba.
To determine whether the presence of more than one sequence for most of the putative peroxins identified in A. polyphaga and A. royreba could be due to problems during the assembly procedure, we evaluated the quality of these genome assemblies using the Benchmarking Universal Single-Copy Orthologs (BUSCO) assessment tool, which measures assembly completeness based on expected gene content. BUSCOs are expected to be found as single-copy orthologs and duplication of these genes is considered as a rare event . When the assemblies were analyzed using the 303 single-copy orthologs from the eukaryotic data set, only six duplicated BUSCOs were found in the A. polyphaga LQHA assembly, whereas 66 duplicated BUSCOs were found in the A. polyphaga CDFK assembly (supplementary fig. 1,  Supplementary Material online). Similarly, 132 duplicated BUSCOs were found in the A. royreba CDEZ assembly. The large number of duplicated BUSCOs found in A. royreba and in the A. polyphaga CDFK assembly suggests that the genome assembly procedure could have failed to collapse the sequenced haplotypes ).
Next, we performed in silico analyses to identify putative peroxin homologues in B. mandrillaris. There are two genome assemblies from different B. mandrillaris strains deposited in the NCBI WGS database, one corresponding to the CDC-V039 isolate under accession number LFUI00000000.1 (Detering et al. 2015) and another corresponding to the 2046 strain under accession number LEOU00000000.1 (Greninger et al. 2015). Neither genome assembly has been fully annotated, thus, we performed TBlastN searches against both B. mandrillaris draft genomes using the protein sequences of various peroxins from H. sapiens and N. gruberi as queries.
As shown in table 1 and in the supplementary table 5, Supplementary Material online, 15 putative peroxins were identified in both B. mandrillaris genome assemblies. The conserved domains specific to each peroxin were also predicted to be present in the proteins identified in B. mandrillaris, and the results of the reciprocal BlastP analyses against the query proteomes confirmed the identity of the putative homologues.
In contrast to A. castellanii, only Pex5, Pex6, Pex7, Pex10, and Pex12 were identified in B. mandrillaris with both the human and the N. gruberi protein sequences. Pex1 and Pex3 were only identified with the human queries, whereas Pex11 and Pex19 were only identified with the N. gruberi queries. Pex2 and Pex16 were also identified with the human queries, but we were unable to identify putative homologues for Pex13 or Pex14 when the respective human sequences were used as queries for the TBlastN analyses. Similarly, Pex4 was identified with the corresponding N. gruberi protein sequence, but not Pex22.
Because the human and the N. gruberi queries appeared to be too divergent from B. mandrillaris to detect peroxin homologues with confidence, we decided to confirm our results by using the sequences of the peroxins identified in A. castellanii, which is more closely related to B. mandrillaris.
When the A. castellanii sequences were used as queries, we were able to confirm the presence of all of the putative peroxins identified with either the human queries, the N. gruberi queries, or both. Additionally, the A. castellanii sequences allowed us to identify putative homologues for Pex13 and Pex14-which had not been found with the human queries-, a putative Pex22 homologue-which had not been identified with the N. gruberi query-and a second member of the Pex11 family, which had not been identified with either query. We identified two sequences for Pex1, Pex7, Pex12, and Pex13 in the LFUI release, whereas two Pex14 sequences were identified in the LEOU release. For Pex11, we found two different sequences in the LFUI release, which represent distinct members of this family (table 1 and  supplementary table 5, Supplementary Material online), as they are only 15% identical at the protein level (supplementary text 1, Supplementary Material online). We also found three potential Pex26 homologues, but given that the sequences were very short, we could not confirm the presence of the Pex26 domain. The fact that only the A. castellanii queries allowed us to identify all of the putative peroxin sequences in B. mandrillaris underscores the importance of using closely related sequences to confirm the presence or absence of putative homologues.
In contrast to B. mandrillaris and the different Acanthamoeba species, only 13 putative peroxins were identified in N. fowleri, as Pex13 and Pex14 were not found in this amoeba (table 1 and supplementary table 6, Supplementary Material online). The putative peroxin sequences were identified by using the corresponding H. sapiens and N. gruberi protein sequences as queries for BlastP analyses against the N. fowleri protein database stored in AmoebaDB. All these sequences were identified with the human queries (except for Pex4 and Pex22, which are absent in humans) and with the N. gruberi queries (except for Pex2 and Pex16, which have not been reported in this amoeba). Intriguingly, a longer version of the Pex2 sequence than the one reported in AmoebaDB was identified when the D. discoideum Pex2 sequence was used as query for TBlastN analyses against the N. fowleri NCBI WGS database. This longer sequence includes both the PEX2_PEX12 and the RING-HC_PEX2 domains (supplementary table 6, Supplementary Material online). Finally, in contrast to N. gruberi-in which only one member of the Pex11 family has been described-two putative members of this family (labeled Pex11a and Pex11b) were identified in N. fowleri. Interestingly, both Pex11 sequences were identified with the N. gruberi query but only one of them was found with the human queries. These putative Pex11a and Pex11b sequences are 19% identical at the amino acid level (supplementary text 1, Supplementary Material online).
Recently, Liechti et al. (2018) characterized the genome of N. lovaniensis which, of the 47 species of Naegleria described to date, is the one that is most closely related to the pathogenic N. fowleri (De Jonckheere 2014). Using the peroxin sequences from D. discoideum, N. gruberi, and N. fowleri as queries for TBlastN searches against the N. lovaniensis WGS database, we identified 13 putative peroxin homologues in this amoeba (table 1 and supplementary table 7, Supplementary Material online). All of these putative peroxin sequences were identified with the D. discoideum queries (except for Pex3, Pex16, and Pex19), the N. gruberi queries (except for Pex2 and Pex16, which have not been described in this amoeba) and with the N. fowleri queries. Importantly, the putative Pex16 homologue, which was not found with the D. discoideum query, was identified when the N. fowleri Pex16 sequence was used as a query for the TBlastN analysis. All of the putative peroxin homologues detected in N. lovaniensis have the expected functional domains and, similarly to N. fowleri, two putative members of the Pex11 family were also identified in this amoeba, which are 20% identical at the amino acid level (supplementary text 1, Supplementary Material online). Also, as observed in N. fowleri, no Pex13 or Pex14 homologues were found in N. lovaniensis.
The fact that N. fowleri and N. lovaniensis had the same number of putative peroxin homologues prompted us to check the N. gruberi genome for sequences coding for Pex2, Pex16, and for another member of the Pex11 family, which had not been reported in this amoeba thus far. Using the corresponding N. fowleri sequences as queries for TBlastN analyses against the N. gruberi NCBI WGS database, we were able to identify putative homologues for these peroxins, and we also found longer versions of Pex1 and Pex6 (table 1 and supplementary table 8, Supplementary Material online). All of these putative peroxins have the expected functional domains and the reciprocal BlastP analyses against the query database confirmed the homology between the sequences. These results constitute the first report on the presence of Pex2 and Pex16 in N. gruberi. As for the two members of the Pex11 family, one of them (ACER01000200.1) has not been reported previously, whereas the second one (ACER01000043.1) has been described as part of the N. gruberi proteome (UniProt ID: D2V0G7). When compared at the amino acid level, these two putative members of the Pex11 family are 22% identical (supplementary text 1, Supplementary Material online). As expected, we were unable to identify homologues for Pex13 or Pex14 in this amoeba, which suggests that they are missing in this genus.
Based on the results of our cytochemical assays, we then performed BLAST analyses to confirm the presence of catalase in the different Amoebozoans and Heteroloboseans included in our study. As expected, we were able to identify several putative homologues of this enzyme in all of these organisms (table 1 and supplementary tables 2-7, Supplementary Material online). Our in silico analyses also predicted the presence of a peroxisomal targeting sequence (PTS1) in at least one of the putative catalase enzymes found in each of the different Amoebozoans included in our study (table 1 and supplementary table 9, Supplementary Material online). For N. fowleri and N. lovaniesis, one of the two putative catalase sequences identified in these amoebae was predicted to have a PTS1 signal, but none of the putative catalase sequences from N. gruberi was predicted to contain this peroxisomal targeting sequence (table 1  and supplementary table 9, Supplementary Material online).
Finally, we performed two phylogenetic analyses using the ATPase Pex1 and Pex5, the receptor for the PTS1 protein import pathway. The phylogeny of these peroxins suggests that A. castellanii is more closely related to A. polyphaga than to A. royreba. In both phylogenies, this grouping had a high statistical support. These organisms were then grouped with B. mandrillaris and this cluster was then grouped with the other Amoebozoans included in the phylogenetic analyses ( figs. 3 and 4). Similarly, the three Naegleria species included in our study were grouped in a distinct clade; however, this cluster was only grouped with the other Excavates in the Pex1 phylogeny ( fig. 3). In both phylogenetic analyses, the pathogenic N. fowleri was found to be more closely related to N. lovaniensis than to N. gruberi (figs. 3 and 4).

Discussion
Peroxisomes are metabolically diverse organelles that are present in most eukaryotic cells. The interest in peroxisomes and in their physiological functions has increased significantly since numerous studies have demonstrated their active participation in crucial metabolic processes such as the b-oxidation of fatty acids, the biosynthesis of ether phospholipids, and the metabolism of reactive oxygen species. These organelles have also been shown to rapidly assemble, multiply, and degrade in response to metabolic needs . The importance of peroxisomes in maintaining cellular homeostasis is also highlighted by the fact that defects in proteins involved in the biogenesis of peroxisomes or in single peroxisomal enzymes can result in inherited peroxisomal disorders (Platta and Erdmann 2007).
Although peroxisomes are involved in numerous cellular processes, some protists can survive without them, including Apicomplexan parasites such as Plasmodium and Cryptosporidium ( Z arsk y and Tachezy 2015) and parasites which lack canonical mitochondria-but which contain mitochondrion-related organelles-such as Entamoeba histolytica, Giardia lamblia, and Trichomonas vaginalis (Makiuchi and Nozaki 2014;Z arsk y and Tachezy 2015). Nonetheless, peroxisomes have been observed in other species of protozoa including Tetrahymena pyriformis, Paramecium caudatum, Euglena gracilis, D. discoideum as well as in Trypanosoma and Leishmania among others (Mü ller 1975;Gabald on 2010;Gabald on et al. 2016). Despite this remarkable diversity regarding the presence or absence of peroxisomes in freeliving and parasitic protists, many protozoan lineages remain understudied. In this respect, information regarding the presence of peroxisomes in opportunistic free-living amoebae is limited and in many cases it is still unknown if these organisms have peroxisomes at all.
By transmission electron microscopy, various round structures with a dark granular content and surrounded by a single membrane were observed in several Acanthamoeba species, as well as in B. mandrillaris and N. fowleri ( fig. 1A-F). Our cytochemical observations indicated that these structures contained catalase, as positive reaction products were observed when samples were incubated with DAB and hydrogen peroxide ( fig. 2A-F). DAB labeling has been used extensively to detect peroxisomes by light and electron microscopy and its widespread use has demonstrated the ubiquity of this organelle in different organisms (Fahimi and Baumgart 1999;Baumgart et al. 2003). Thus, the presence of positive reaction products for this cytochemical staining, together with our ultrastructural observations, provides experimental evidence for the existence of peroxisomes in these free-living amoebae. To the best of our knowledge, our ultrastructural observations in N. fowleri constitute the first morphological and cytochemical evidence for the presence of peroxisomes in a member of the Heterolobosea (Gabald on et al. 2016). Moreover, these results represent the first morphological and cytochemical evidence for the presence of peroxisomes in B. mandrillaris as well as in A. polyphaga, A. royreba, and A. griffini.
Our ultrastructural observations in the different Acanthamoeba species are complemented by earlier biochemical studies which found that peroxisomes are present in this genus. Childs (1973aChilds ( , 1973b reported the presence of peroxisomes in Hartmannella culbertsoni (now known as Acanthamoeba culbertsoni) by using DAB and hydrogen peroxide for catalase localization, whereas Mü ller and MØller (1969) were the first to describe the presence of peroxisome-like particles containing urate oxidase and catalase in a free-living amoeba by isopycnic centrifugation of Acanthamoeba trophozoites. These enzymes were also found in the slime mold D. discoideum (Parish 1975;Hayashi and Suga 1978) which, together with Acanthamoeba, was one of the first members of the Amoebozoa in which the presence of peroxisomes was described by biochemical assays.
In recent years, the genome sequences of numerous organisms from different taxonomic groups have become available, allowing the use of comparative genomics to identify putative peroxisomal proteins based on a high degree of sequence similarity (Wanders and Waterham 2006;Gabald on 2018). Because peroxins are common to all peroxisomes despite their enzymatic content, their genomic identification can be considered as an indicator for the existence of peroxisomes in a particular organism (Gabald on 2018). This bioinformatic approach has provided compelling evidence regarding the presence or absence of this organelle in different organisms, including several protist lineages (Moog et al. 2017;Ludewig-Klingner et al. 2018).
Our comparative genomic analyses allowed us to identify a group of 15 putative peroxins in the Amoebozoans A. castellanii, A. polyphaga, A. royreba, and B. mandrillaris (table 1). Among these putative peroxins, we identified Pex5, which is a cytosolic receptor for proteins destined to the peroxisome that bear the peroxisomal targeting signal type 1 (PTS1) and Pex7, which is a soluble receptor for proteins that carry the PTS2 sequence. We also found putative homologues for Pex13 and Pex14, which are membrane proteins that form a docking complex that interacts with the receptors to translocate their cargo proteins into the peroxisomal matrix. The ubiquitin conjugating enzyme Pex4 and its membrane anchor, Pex22, as well as the RING ubiquitin ligases Pex2, Pex10, and Pex12-which monoubiquitinate Pex5 for recycling or polyubiquitinate it for degradation in the proteasome-were also predicted to be present in these amoebae. Likewise, the AAA-type ATPases Pex1 and Pex6, which recycle the monoubiquitinated Pex5 back to the cytoplasm, were identified in these organisms. We also found putative homologues for Pex3, Pex16, and Pex19, which are involved in the insertion of membrane proteins into the peroxisomal membrane. Finally, several members of the Pex11 family, which is essential for peroxisome division, were also identified in these Amoebozoans (peroxins functions are described in Girzalsky et al. [2010], Kim and Hettema [2015], and Platta et al. [2016]). The fact that at least two distinct members of the Pex11 family were found in these Amoebozoans is not uncommon, as several Pex11 variants have been described in other organisms. Arabidopsis thaliana has five Pex11-related proteins (Lingard and Trelease 2006), whereas three distinct members of this family have been identified in mammals and filamentous fungi (Thoms and Erdmann 2005;Kiel et al. 2006). Besides being fundamental for peroxisome proliferation, Pex11 also regulates peroxisome morphology, number, and size . Although their functions are partially redundant, the different members of this family appear to have some degree of functional variability: in plants, two Pex11 variants have been shown to increase the number of peroxisomes, whereas others induce elongation or peroxisome aggregation (Lingard and Trelease 2006;Schrader et al. 2012). It will be interesting to determine whether the different members of the Pex11 family identified in these Amoebozoans have overlapping or independent functions, as well as their specific role during peroxisome proliferation.
Besides these 15 putative peroxins, a possible homologue for Pex26 was also identified in these Amoebozoans. Pex26 is present in mammals but it is absent in plants as well as in Saccharomyces cerevisiae and related yeasts (Islinger et al. 2010). In mammals, Pex26 anchors the Pex1-Pex6 complex to the peroxisomal membrane (Matsumoto et al. 2003), a function that is performed by Pex15 in yeasts (Kalel and Erdmann 2018). The Pex26 domain was only identified in the putative A. castellanii homologue and in just one of the three putative homologues found in A. polyphaga (supplementary tables 2 and 3, Supplementary Material online). None of the potential Pex26 homologues described in A. royreba bears the corresponding domain, despite the fact that these sequences had high TBlastN and reciprocal BlastP E-values (supplementary table 4, Supplementary Material online). It is possible that the sequence for this domain has diverged significantly from that of other organisms, which would explain why we were unable to identify it. Alternatively, these putative homologues might have become redundant and may be evolving new functions or they may not be bona fide Pex26 homologues. Importantly, we could not identify Pex15 homologues in any of these amoebae. Overall, this group of predicted peroxins is probably sufficient for the generation and maintenance of functional peroxisomes in A. castellanii, A. polyphaga, A. royreba, and B. mandrillaris.
Our bioinformatic data regarding the presence of peroxins in Acanthamoeba and B. mandrillaris indicate that the peroxin content of these organisms is similar to that of other Amoebozoans, because all of the putative peroxins homologues found in Acanthamoeba and B. mandrillaris-with the exception of Pex22 and Pex26-have also been described in D. discoideum (supplementary table 1, Supplementary Material online). As more genomic data from other Amoebozoan lineages become available, future studies will provide additional insights regarding the peroxin content of this taxonomic group.
With the exception of Pex13 and Pex14, a complete set of putative peroxins was also identified in the pathogenic N. fowleri and in the closely related but nonpathogenic N. lovaniensis (table 1). Our bioinformatic analyses also found putative homologues for the previously unidentified Pex2, Pex16, and for a second member of the Pex11 family in the nonpathogenic N. gruberi (supplementary table 8, Supplementary Material online), confirming that all three Naegleria species share the same set of peroxins regardless of their pathogenic capacity. We also failed to identify a Pex26 homologue in these Naegleria species, which suggests that this peroxin might not be necessary for peroxisome function in these organisms. Thus, this group of Heteroloboseans appear to have a slightly reduced set of peroxins with respect to the Amoebozoans analyzed in this study.
The absence of Pex13 and Pex14 in N. gruberi, N. fowleri, and N. lovaniensis raises the question of whether these peroxins were present in Naegleria and were subsequently lost during evolution or whether these amoebae developed an import mechanism that does not require this particular docking complex. Typically, Pex13 and Pex14 form a membranebound protein complex that interacts with Pex5 and Pex7 to import proteins into the peroxisomal matrix. Given that Pex5 and Pex7 are present in all three Naegleria species, future studies should aim to understand how this import mechanism functions in the absence of the Pex13-Pex14 docking complex. Alternatively, the Pex13 and Pex14 sequences could have diverged significantly in these Heteroloboseans, to the point that they could not be detected by our in silico analyses. Importantly, the absence of Pex13 and Pex14 has been reported in other organisms that contain peroxisomes such as the diatoms Thalassiosira pseudonana (Schlü ter et al. 2006) and Phaeodactylum tricornutum (Gonz alez et al. 2011). However, these organisms also lack Pex7, which indicates that they have lost the PTS2 import pathway (Gonz alez et al. 2011). Because all three Naegleria species are predicted to have Pex7 homologues and N. gruberi has been reported to have at least one peroxisomal protein with a predicted PTS2 signal (Gonz alez et al. 2011), it seems unlikely that this import pathway is absent in Naegleria.
As shown in table 1, with the exception of Pex13 and Pex14, our comparative genomic results support previous observations which indicate that, of the more than 30 peroxins described to date, at least 13 of them (5)(6)(7)(10)(11)(12)(13)(14)16,and 19) are present in all major eukaryotic lineages ( Z arsk y and Tachezy 2015). Besides this core group of peroxins, we found putative homologues for Pex4 and Pex22 in all of the Amoebozoans and Heteroloboseans included in our study. Both of these peroxins are absent in mammals, but they are expressed in plants, filamentous fungi, and yeasts (Islinger et al. 2010;Platta et al. 2016).
In contrast to the other free-living amoebae analyzed in this study, two or more putative copies for almost every peroxin were found in A. polyphaga (supplementary table 3, Supplementary Material online) and A. royreba (supplementary table 4, Supplementary Material online). When the completeness of the genome assemblies used in this study was evaluated using the BUSCO software, we found that the A. polyphaga CDFK assembly, the A. royreba CDEZ assembly, and the B. mandrillaris LFUI assembly had a large number of duplicated BUSCOs compared with the other assemblies (supplementary fig. 1, Supplementary Material online). Although the precise ploidy levels of these Amoebozoans have not been determined, there are strong indications that they are polyploid (Matsunaga et al. 1998;Detering et al. 2015). Thus, the multiple gene copies that were found in these assemblies might be heterozygous alleles that the assembly procedure failed to collapse . Nonetheless, additional studies such as genomic and expression analyses are needed to confirm that these putative copies are indeed sequence assembly artifacts.
To the best of our knowledge, none of the putative peroxins identified in this study has been previously reported in A. polyphaga, A. royreba, B. mandrillaris, or N. lovaniensis. In contrast, most of the putative peroxins described in N. fowleri are properly annotated in AmoebaDB or in the UniProt database, but most of the A. castellanii and N. gruberi peroxins are not annotated as such and instead they are classified according to their domains or as hypothetical proteins.
Our bioinformatic analyses also allowed us to confirm the presence of several putative catalase sequences in all of these organisms (table 1 and supplementary tables 2-7, Supplementary Material online). For each of these organisms, at least one putative catalase sequence was predicted to be located in the peroxisome, as they contain a PTS1 signal (supplementary table 9, Supplementary Material online). These results support our previous cytochemical observations which indicated that catalase was present in the peroxisomes of A. castellanii, A. polyphaga, A. royreba, B. mandrillaris, and N. fowleri ( fig. 2). All of the putative catalase sequences from B. mandrillaris and from the different Acanthamoeba species had the same PTS1 signal (-AKL), which suggests that this sequence is conserved in these Amoebozoans. As for the Heteroloboseans analyzed in this study, only one putative catalase sequence from N. fowleri and one putative catalase sequence from N. lovaniesis were predicted to have a PTS1 signal (-NNL). The -NNL tripeptide has been identified as a putative PTS1 sequence in the sterol carrier protein 2 (SCP2) from the yeast Yarrowia lipolytica (Ferreyra et al. 2006). Although this tripeptide has not been experimentally validated as a functional PTS1, fungal SCP2 is strictly peroxisomal and it is believed to be involved in the peroxisomal oxidation of long-chain fatty acids (Ferreyra et al. 2006). In contrast to N. fowleri and N. lovaniensis, none of the putative catalase sequences reported in N. gruberi appears to have a PTS1, although one of them has a C-terminal tripeptide that is similar to the one identified in the other Naegleria species (-QNL). Importantly, Opperdoes et al. (2011) also found that the predicted N. gruberi catalase sequence lacks a PTS signal. Interestingly, the PTS1 sequence identified in N. fowleri and N. lovaniensis (-NNL) differs from the predicted PTS1 signals identified by Opperdoes et al. in other N. gruberi peroxisomal enzymes, which were found to be predominantly -SKL or a variant thereof (-NKL/-SKM).
It should be noted that the bioinformatic approach used for identifying peroxisomal homologues in these Amoebozoans and Heteroloboseans has some limitations. The putative protein sequences of the peroxins and catalases identified in A. polyphaga, A. royreba, B. mandrillaris, N. gruberi, and N. lovaniensis were determined using the GeneWise algorithm, which has been used both as a standalone tool and as part of a gene prediction pipeline (Birney et al. 2004). GeneWise identifies potentially homologous genes by aligning a genomic DNA sequence to a homologous protein sequence (Birney et al. 2004). This tool is able to identify complete genes when the homologues are very similar, but when the homologues are more distant, GeneWise is only able to detect part of the exon structure (Guig o et al. 2000;Yeh et al. 2001). In fact, a known disadvantage of this algorithm is that it fails to detect complete terminal exons, as GeneWise predictions are limited to sites where there is evidence of similarity (Birney et al. 2004). Most of the putative protein sequences identified in A. polyphaga, A. royreba, N. lovaniensis, and N. gruberi appear to be fairly complete, as these sequences have a high level of similarity to their respective queries. However, in some cases, GeneWise could only find fragments of these putative proteins, as observed for the A. royreba Pex22 sequence and for the B. mandrillaris Pex14, Pex16, Pex19, and Pex22 sequences (supplementary tables 4 and 5, Supplementary Material online). Despite being incomplete, all of these putative sequences bear the expected PEX domains, which suggests that the aligned region was limited to the conserved domains and that the queries and the putative proteins have a low level of similarity (Guig o et al. 2000). Thus, further analyses are needed to obtain the complete sequences of the putative peroxisomal proteins identified in our bioinformatic surveys. Likewise, as our studies were performed using genomic data, transcriptomic and proteomic studies are needed to confirm the expression of these putative peroxisomal proteins.
With respect to the phylogenetic analyses of Pex1 and Pex5, A. castellanii was found to be more closely related to A. polyphaga than to A. royreba ( figs. 3 and 4). Importantly, a previous phylogenetic reconstruction based on the 18S ribosomal genes of various Acanthamoeba species also grouped A. castellanii with A. polyphaga, whereas A. royreba was clustered with other Acanthamoeba species (Chelkha et al. 2018). The Acanthamoeba cluster was then grouped with B. mandrillaris, confirming the close phylogenetic relationship between Acanthamoeba and Balamuthia ( figs. 3 and 4). These results are also in agreement with previous phylogenetic analyses of the mitochondrial cox1 protein and the 28S rRNA gene, which also clustered these two genera together (Greninger et al. 2015). As for the Heteroloboseans included in our study, the pathogenic N. fowleri was found to be more closely related to N. lovaniensis than to N. gruberi, as reported in previous studies (De Jonckheere 2014;Liechti et al. 2018). Importantly, only the Pex1 phylogeny grouped the Naegleria cluster with the other Excavates included in the analysis (figs. 3 and 4).
In summary, our morphological, cytochemical, and genomic data demonstrate that peroxisomes are present in B. mandrillaris, A. castellanii, A. polyphaga, and A. royreba. Our cytochemical and morphological data also indicate that A. griffini has peroxisomes, but because the genome of this amoeba has not been sequenced, we could not determine if some or all of the peroxins and catalase homologues that were identified in the other Acanthamoeba species are also present in A. griffini.
Our ultrastructural, cytochemical, and comparative genomic data also indicate that N. fowleri has bona fide peroxisomes. Similarly, our comparative genomic results indicate that peroxisomes are also present in the nonpathogenic N. lovaniensis, although morphological observations are still needed to confirm the presence of this organelle. Additionally, it will be important to explore the enzymatic content of the peroxisomes of these Amoebozoans and Heteroloboseans to establish their precise role in their metabolism and to increase our understanding of the diversity of peroxisomal functions in different eukaryotic lineages.

Supplementary Material
Supplementary data are available at Genome Biology and Evolution online.