{'fraction_non_alphanumeric': 0.024926319996003796, 'fraction_numerical': 0.013087566811529047, 'mean_word_length': 5.378684084753863, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A large scale, multiple genome comparison of acidophilic Archaea (pH ≤ 5.0) extends our understanding of oxidative stress responses in polyextreme environments

Acidophilic Archaea thrive in anaerobic and aerobic low pH environments (<pH 5) rich in dissolved heavy metals that exacerbate stress caused by the production of reactive oxygen species (ROS) such as hydrogen peroxide (H2O2), hydroxyl radical (•OH) and superoxide (O2 •−). ROS react with lipids, proteins and nucleic acids causing oxidative stress and damage that can lead to cell death. Herein, genes and mechanisms potentially involved in ROS mitigation are predicted in over 200 genomes of acidophilic Archaea with sequenced genomes. These organisms can be subjected to simultaneous multiple stresses such as high temperature, high salinity, low pH and high heavy metal loads. Some of the topics addressed include: (1) the phylogenomic distribution of these genes and what can this tell us about the evolution of these mechanisms in acidophilic Archaea; (2) key differences in genes and mechanisms used by acidophilic versus non-acidophilic Archaea and between acidophilic Archaea and acidophilic Bacteria and (3) how comparative genomic analysis predicts novel genes or pathways involved in oxidative stress responses in Archaea and possible Horizontal Gene Transfer (HGT) events.


Introduction
Reactive oxygen species (ROS) are formed as by-products of aerobic metabolism including, for example, Fenton chemistry [1,2] and can involve molecules such as the superoxide anion (O 2 -), hydrogen peroxide (H 2 O 2 ), and the hydroxyl radical (OH•). The accumulation of these molecules can lead to oxidative stress and produce damage to many cellular macromolecules such as DNA, RNA and proteins that, in turn can impair the function of important cellular structures such as membranes [3][4][5][6]. To mitigate these possible damages, different mechanisms have been developed by prokaryotes. Examples are superoxide scavengers such as superoxide dismutase (SOD) and superoxide reductase (SOR) that can diminish the quantity of superoxide anion [7]and a number of peroxide scavengers like peroxiredoxins (Prx) and catalases as examples that reduce the levels of hydrogen peroxide [8].
Whereas extensive research has been carried out on oxidative stress mechanisms in the bacterial domain [9,10], including the report of an inventory of genes across the whole domain [11], less is known about stress responses in Archaea. We aim to identify the potential mechanisms of oxidative stress response in extremely acidophilic Archaea and to analyze the phylogenetic distribution of predicted response mechanisms across the domain.
Since Archaea have been proposed to be linked with early Earth environments [12,13], their study could provide useful insight into the evolution of oxidative stress responses. Early evolution most likely proceeded in anaerobic conditions [14] but Archaea could be exposed to "whiffs" of oxygen that could lead to the need of the first oxidative stress response mechanisms [15,16], then additional mechanisms could have evolved as they transitioned from the anaerobic world to an aerobic atmosphere product of the great oxidation event (GOE) 2.4 billion year ago [17].
The archaeal domain comprises a wide range of organisms, many of which thrive in environments that are considered polyextreme, such as extremely high temperatures, low and high pHs and high salinity conditions [18][19][20]. Within these groups, stress response mechanisms have been described for hyperthermophiles [21] that constitute the majority of all identified Archaeal organisms. Acidophiles, are commonly defined as organisms that grow optimally at pH lower than 5 [22,23]. Another useful classification includes the categories of moderate acidophiles and extreme acidophiles, with the first group including organisms that grow optimally between pH 3 and pH 5 [24], and the second emphasizes organisms that grow optimally at pH values lower than 3 [25] as the environmental constraints that they face are particular to this pH range [26].
However, despite this progress, a large scale, more comprehensive understanding of stressassociated mechanisms in acidophilic Archaea has not been undertaken. Our report begins to fill in this lacuna in our knowledge, providing a more extensive inventory of oxidative stress responses and helping to identify mechanisms provoked by multiple simultaneous stresses such as high temperature, extreme acidity, and high metal loads. Since Archaea have been proposed to root at the base of the Tree of Life, their study could provide useful insight into the evolution of oxidative stress responses on exposure to regular oxygen intrusion as the world transitioned from an anaerobic to aerobic atmosphere.

Genome collection
A dataset comprising 234 acidophilic (optimal growth <pH 5) archaeal genomes was obtained from the database AciDB [36]. All genomes were downloaded from the National Center for Biotechnology Information (NCBI) and the Joint Genome Institute (JGI). Only high-quality genomes (draft or complete) were considered for the analyses, defined as a >80% completeness and <5% contamination calculated by CheckM v1.0.12 with the standard lineage workflow [37]. Existing genome annotations in Refseq [38] were used when available and, for genomes with no annotation, Prokka v1.13.3 was used to identify genomic features [39]. Growth condition metadata for each genome was obtained from AciDB using the optimum values for pH and temperature. Oxygen requirement information was collected from bibliographic research associated with each species. If the data for a specific species was not available, it was extrapolated from the genus description. This leads to a classification of all organisms into aerobic, anaerobic, or facultative. A summary of genomic information is provided in Table S1.

Phylogeny
Phylogenetic reconstruction was performed using Phylophlan 3.0 using the provided dataset of 400 conserved proteins [40]. The parameters were set as following: diversity as high (as recommended for phylum wide phylogenies); marker search performed using diamond (v2.0.2) [41]; sequence alignment by MAFFT (v7.310) with the L-INSI iterative refinement option [42] and alignment trimming by trimal (v1.2) [43]. Tree construction was performed using IQTREE v1.6.1. The best evolutionary model was identified by ModelFinder according to the Bayesian (BIC) and Akaike information criterion (AIC) [44]. To assess the support of each branch, the non-parametric ultrafast bootstrap method (with 1000 replicates) was used [45]. Phylogenetic trees were visualized and annotated using iTOL [46]. A species clustering analysis was performed using Fast-Ani with a value of 95% identity to classify genomes as the same or different species [47].

Functional annotation and orthologous groups identification
A set of 52 proteins associated with oxidative stress responses was identified from the literature. Functional annotations based on Pfam [48] and Interprofamilies (domain and superfamily) [49] for each protein were obtained (Table S2). All genomes in the study were functionally annotated using the InterProScan pipeline with default parameters [50]. Orthogroups defined as proteins that have evolved from a common ancestor and include both orthologs (homologous protein from two species) and paralogs (a set of proteins that have a common origin in the same genome) were identified using OrthoFinder v2.3.3 [51]. Protein sequences were compared using Diamond in an all-versus-all search and orthogroups were inferred [41]. DendroBLAST [52] was used to generate unrooted gene trees for each of the identified orthogroups. MAFFT was used for multiple sequence alignments and FastTree [53] was used for the tree inference. All other steps were performed with default parameters. An inhouse Python script was used to identify oxidative stress response proteins in the acidophilic genomes based on the functional annotation. Pfam annotations provided a broad scope view of candidates which were subsequently hand curated by InterPro family and domain annotation to obtain additional specific information. Subcellular locations were assigned to predicted proteins and classified as cytoplasmatic, inner membrane, exported, outer membrane or periplasmic using PSORTb v3.0 [54] and accompanied by signal peptide identification performed by SignalP v5.0 [55].

Evolutionary trajectory analysis
Phylogenetic trees were constructed using the 50 best hits from a blastp search against the nonredundant database from all identified proteins in their respective orthogroups. Proteins from the blastp search were downloaded using the Batch Entrez web tool in NCBI [56]. Sequence alignments were performed using MAFFT and were visualized with Aliview v1.2.6 [57] to identify the conservation of key amino acids when needed. The phylogenetic tree was then constructed using IQTREE with an ultrafast bootstrap of 1000 replicates. Genome neighborhoods were analyzed using Gene Graphics with a region size of 5000bp [58].

Phylogeny and Species clustering
High quality draft and complete genomes from acidophilic Archaea were obtained from AciDB [36] with a total number of 180 different genomes that fulfill the over 80% completeness and less than 5% of contamination defined by CheckM [37]. Of these genomes, 47 are complete and 133 are permanent drafts (Table S1). A phylogenetic tree was inferred using a set of 400 conserved markers as explained in material and methods. Three organisms from the ARMAN clade (Ca. Mancarchaeum acidiphilum Mia14, Thermoplasmatales archaeon A_DKE and Ca. Micrarchaeum sp. AZ1) were excluded from the tree inference as these organisms lack the majority of the conserved markers. Four different phyla are represented in the phylogenetic tree Crenarchaeota, Euryarchaeota, Marsarchaeota and Thaumarchaeota (Figure 1). At the genus level, Crenarchaeota and Euryarchaeota are the most diverse phyla in the dataset with ten and six different genera respectively. The species clustering analysis based on ANI values reveals that Cuniculiplasma divulgatum is the Euryarchaeota with the most genomes associated (19 in total), as a series of genomes annotated as "Thermoplasmatales archaeon" were identified as part of this species. Sulfolobus acidocaldarius from Crenarchaeota has the most genomes of any species with 54 genomes associated with it. These results are summarized in the Table S3. The lack of markers in the ARMAN clade organisms is not surprising as this clade is characterized by their small genome sizes derived from their symbiotic lifestyle [59]. The species clustering highlights the importance for the development of new standards in the naming of prokaryotic genomes [60,61] that are not from isolated sources, as the adoption of higher rank taxonomic names could lead to confusion in subsequent analysis.

Figure 1.
Phylogenetic tree reconstruction of acidophilic Archaea with genomes available in AciDB using Phylophlan 400 markers. Branches are collapsed at genus level, except for Sulfulobus which is collapsed at species level. Colors represent different phyla associations with yellow for Crenarchaeota, red for Thaumarchaeota, green for Marsarchaeota and blue for Euryarchaeota. Bootstrap values over 60% are represented by black dots. Abbreviations: Uncultured (Un.); Candidatus (Ca.). Scale bar represent 0.1 amino acid substitution per site.

Superoxide defense proteins
Superoxide anion (O 2 •) is a highly reactive species that can cause severe damage inside the cell. The first system used in the process of removing O 2 • is superoxide dismutase (SOD) that produces hydrogen peroxide and oxygen. Superoxide dismutase is classified in four different categories depending on the cofactor used (Mn, Fe, Cu/Zn and Ni). The second system is associated with the superoxide reductase (SOR) which produces hydrogen peroxide and protons in the reduction of superoxide. Details on how both mechanisms work are reviewed by Sheng and colleagues [62].

Superoxide dismutase (SOD)
We found a wide distribution of Fe-SOD throughout several acidophilic archaeal genera from all phyla in study, despite being previously identified in only a few organisms such as Saccharolobus solfataricus [63], Ferroplasma acidiphilum [64] and Metallosphaera sedula [65] (Figure S1). The abundance of this category of SOD (Fe-SOD type) could be of evolutionary importance given that the iron containing isoform has been suggested to be an ancient type of SOD [7]. The higher availability of iron in low pH environments, together with the notion that Archaea are ancestral organisms [12], highlights these acidophiles as prime candidates to be the original source of one of the most relevant oxidative stress mechanisms used across all domains of life for evolutionary studies and biotechnological applications. Phylogenetic analysis of the predicted sequences was performed to study the evolution of SOD in these organisms.The sequences were divided in four different clades, two clades that grouped most of the sequences associated with the phylum Crenarchaeota and Euryarchaeota and two smaller clades formed by sequences of organisms from the phylum Marsarchaeota and Thaumarchaeota, the latter including outliers of the Crenarchaeota phylum from the Acidilobus and Caldisphaera genera. The distribution in each of the separated clades appears to be associated with phylogenetic relatedness, as each predicted protein branch is formed mainly by organisms from the same genus with a clear distinction between the sequences from Euryarchaeota and those from the phyla Crenarchaeota, Thaumarchaeota and Marsarchaeota which are phylogenetically associated as part of the TACK superphylum [66] (Figure 2).
An exception to the phylogenetic association can be observed, specifically, Ca. Micrarchaeum sp. AZ1 and Ca. Mancarchaeum acidiphilum Mia14 sequences, which are both inside the Crenarchaeota clade with the former close to Thermoproteus and the latter to Acidianus (Figure 2, red asterisks). This difference could be explained by the fact that Micrarchaeota organisms have reduced genome size and limited metabolic capabilities, relying on other members of the community to supply these deficiencies by a symbiotic approach [67]. The observation that all SODs from these organisms are closely related with organisms that thrive in the same environments suggest that this oxidative stress mechanism probably was gained by HGT, a process which has been observed to be essential in the evolution of other members of the DPANN clade [68].

Figure 2.
Superoxide dismutase phylogenetic tree from predicted protein sequences for gene sodA. Two Candidatus Micrarchaeota leaves are marked with red asterisks which are in the bigger clade of Crenarchaeota. Colors represent different phyla associations with yellow for Crenarchaeota, red for Thaumarchaeota, green for Marsarchaeota and blue for Euryarchaeota. Scale bar represents 1 amino acid substitution per site.

Superoxide reductase (SOR)
The superoxide reductase (SOR) was predicted only in the genus Aciduliprofundum, which is the only genus that lacks any SOD sequence. Phylogenetic analysis of the sequences identified in Aciduliprofundum and archaeal proteins from SORGOdb [69] shows a clade formed with other strict anaerobes from Euryarchaeota (non-acidophiles), suggesting that vertical transmission could play a key role in the evolution of this protein ( Figure S2). Previous studies have shown that SOR is mainly found in anaerobic organisms [70] which is the case of Aciduliprofundum, increasing the interest for other anaerobic organisms in the study (for example Vulcanisaeta and Acidilobus) which use Fe-SOD instead as main medium to get rid of O 2 •. Rubredoxin (Rb) has been identified as an electron donor for the reduction of superoxide via SOR in Pyrococcus furiosus (also from the Euryarchaeota phylum) [71] and was also identified only in Aciduliprofundum, supporting the idea that this organism is the only one in dataset able to reduce superoxide. One possible explanation for this contradiction with previous studies that claim that the main mechanism of coping with superoxide in anaerobic organisms is the use of SOR (avoiding the production of O 2 resulting from SOD) are the special conditions in which acidophiles thrive. SOR in the reduction of O2• produces extra protons (which instead are not produced by SOD), leading to additional stress for acidophilic organisms that need to maintain a neutral intracellular pH and balance the extremely high concentration of protons in the environment. This produces a disadvantage to the specific use of SOR in low pH environments that could favor the use of SOD, which has been previously observed as one of the main oxidative stress response mechanisms in acidophiles [72], resulting in the loss of the SOR gene through adaptation to the harsh environmental conditions.

Catalase
Catalases convert hydrogen peroxide into water and oxygen and work as a common oxidative stress mechanism in organisms from all domains of life [11]. In acidophilic Archaea, only genomes from Thermoplasmatales and Aciduliprofundum sp. MAR08-339 contain predicted catalase-peroxidase (with the same functional annotation as KatG). Phylogenetic analyses with proteins from NCBI show that the sequences from Thermoplasmatales form a clade close to extremely acidophilic bacteria from the genus Methylacidimicrobium which are also moderate thermophilic organisms. In the case of the KatG identified for Aciduliprofundum sp. MAR08-33, this forms a separated clade with other Euryarchaeota organisms (not acidophiles) and some Deltaproteobacteria ( Figure S3). The low presence in Archaea has been reported previously, as organisms that have SOD usually do not have catalases [21] and as previously discussed we identified SOD in most acidophilic genomes with the exception of Aciduliprofundum. The result of mixed origin clades in the phylogenetic reconstruction where Archaea and Bacteria share the same clade suggests that HGT between domains could play a major role in the evolution of catalases. This goes in line with studies that show that the phylogeny of catalases is complex and often does not follow a pattern of vertical inheritance [73,74].

Osmotically inducible protein C (OsmC)
Another protein with peroxide scavenging activity is the osmotically inducible protein C (OsmC) which has been shown to be involved in the elimination of organic peroxides in Bacteria [75]. OsmC was found to be ubiquitous in acidophilic Archaea with most organisms in Crenarchaeota, Euryarchaeota and Marsarchaeota predicted to have one or more copies, even reaching 3 different copies in Metallosphaera organisms. Phylogenetic tree analysis of the identified sequences reveals a clear division between the sequences that evolved from Euryarchaeota and Crenarchaeota. Nine different orthologous groups were formed by the predicted OsmC sequences. The orthogroups mostly match the clade as observed in Figure S4, where each orthogroup contains proteins from the same phylum except for OG2461 and OG3045 which also contain sequences from Marsarchaeota. Protein sequences from Marsarchaeota were observed to be distributed across clades with sequences from both Crenarchaeota and Euryarchaeota, forming separated clades suggesting HGT between these phyla. KAD1 in Thermococcus kodakaraensis, is an OsmC protein that can reduce inorganic peroxides as well as organic ones [76] and it is interesting to note that the different orthogroups identified across acidophilic Archaea could be related with functional distinction to the classic OsmC. Further experimental investigation is needed to unravel the functional activities of these different identified predicted proteins.

DNA-Binding Protein from Starved Cells
Another way to prevent the creation of ROS is sequestering iron, thus decreasing the Fenton reactions. A protein that fulfills this function is the DNA-Binding Protein from Starved Cells (Dps), a member from the ferritin-like superfamily which can also reduce hydrogen peroxide in the process of oxidation of free ferrous iron [77]. Dps in Archaea can be divided into two different functional annotations, the first one with the InterPro family IPR002177 group from the model protein for Dps in Bacteria and a second InterPro family IPR014490 domain that represents Dps-like proteins from Archaea first identified in Saccharolobus solfataricus [78]. Using InterPro family IPR002177 domain as a probe (Bacteria associated), we predicted Dps in Marsarchaeota and Sulfolobales MK5 as part of a single orthogroup. Unlike the previous case, the Dps-like family (Archaea associated) was identified in three different phyla, Crenarchaeota (Metallosphaera, Saccharolobus and Sulfobolus), Euryarchaeota (Aciduliprofundum) and Thaumarchaeota (Ca. Nitrosotalea). To study the possible evolutionary origin for the Dps protein, a phylogenetic tree was constructed using the best hits from a blastp search. The Dps phylogeny shows their closest sequences include other acidophiles coming from Acidobacteria and several Proteobacteria suggesting a likely cross domain HGT as the main evolutionary force for this protein in Archaea that could be associated with a specific eco-niche adaptation [79] (Figure 3). In contrast, Dps-like protein is widely found in Archaea from Crenarchaeota and Thaumarchaeota. The functionality of the Dps-like protein as an actual Dps protein (instead of other ferritin-like superfamily protein) could be inferred from the high similarity of the sequences with the functionally active protein previously identified in Saccharolobus solfataricus [78]. In addition, we could not identify any homologous proteins to either Dps or Dps-like to be widely present in Euryarchaeota, suggesting that this mechanism could be replaced in the oxidative stress response by another protein of similar function.

Rubrerythrin (Rbr)
Rubrerythrin (Rbr, rbr), another member of the ferritin-like superfamily that consists of varied non-heme di-iron proteins, has an important role in oxidative stress response by reducing hydrogen peroxide [80]. Rbr was identified in almost all acidophilic Archaea both in aerobic and anaerobic organisms. The identified proteins with rubrerythrin domains are annotated in several different orthogroups, forming distinct clades in the phylogenetic analysis. Five different domain architectures were identified across the different proteins. The first group (associated with OG0767) is characterized by a single and isolated Rbr domain with conserved di iron metal motifs, this group was identified across all major phyla in study. The second group (associated with OG0732) is characterized by a domain architecture of a full Rbr domain with di-iron metal motifs together with a VIT1 domain [81]. This group is also conserved across all major phyla in study. A third domain (associated with OG6126, OG7923 and OG3502) similar to the first group identified contains a small truncated Rbr domain at the C-terminal end and is exclusive to Euryarchaeota. The fourth group (associated with OG1026 and OG3704), presents a truncated Rbr domain with missing di-iron metal motifs and with the addition of Linocin M18 domain. The final domain architecture group (associated with OG1052) presents a small Rbr domain without any di-iron metal motifs. The presence of the di-iron metal motifs is important in the context of oxidative stress response, as these are essential in the action of Rbr as peroxidase [82]. From these analyses, we propose that group 1 and 2 could be directly related with hydrogen peroxide reduction. In particular, the addition of the domain VIT1 in group 2 has been previously reported to aid in iron detoxification via transportation of the metals, a process that prevents oxidative stress via Fenton reaction products in plants, yeast, and fungi [83,84]. The presence of a Linocin M18 domain in group 3 could help to the formation of nano compartments that encapsulate peroxidases [85]. Previously, this variant was only identified in the extreme thermophile Pyrococcus furiosus [80], where the specific action has been hypothesized to be related with the hyperthermophilic trait of these organisms, in our analysis this variant appears not to be exclusive to hyperthermophiles; Vulcanisaeta and Caldivirga also have this domain architecture and grow at a more moderate temperature (45-50 ºC). (Figure  4).

Peroxiredoxins
Peroxiredoxins (Prx) have been identified as the major mechanism of hydrogen peroxide elimination in Archaea as reviewed in Pedone et al. 2020. These enzymes have been classified in 6 different subfamilies (AhpC-Prx1, Bcp-PrxQ, Tpx, Prx5, Prx6 and AhpE) using structural information near the active site [86]. A previous study in Saccharolobus solfataricus characterized four different peroxiredoxins from these organisms, three of them (described as Bcp1-3-4) belonging to the Bcp-PrxQ subfamily and one (Bcp2) to the Prx6 [87]. All these sequences have their crystal structure solved, providing experimentally proved and classified examples of peroxiredoxins in Archaea [88]. Peroxiredoxins are widely distributed across acidophilic Archaea with a varying number of copies depending on the organisms (number of copies mostly conserved at genus level). The phylogenetic distribution of the different enzymes identified across the acidophilic organisms is presented in Figure 5a showing four larger clades that include the proteins characterized in Saccharolobus (Prx6/Bcp2, Bcp1, Bcp3 and Bcp4) and a fifth new clade not previously reported in acidophilic Archaea with sequences corresponding to the Prx1 subfamily. The clades formed by Bcp1 (PrxQ family) and Bcp 2 (Prx6 family) sequences are widely distributed across organisms from all the phyla in study, which shows a high degree of conservation of both families of Prx across all acidophilic Archaea. In contrast, Bcp3 is exclusive to organisms of the Crenarchaeota phylum. The Bcp4 clade is divided into three different sub-clades with the first divergent clade showing a similar distribution to the Bcp3 clade with Vulcanisaeta, Thermoproteus and Caldivirga, a second clade containing sequences from Euryarchaeota (that includes the ARMAN organism Thermoplasmatales archaeon B_DKE) and finally a clade that include the rest of Crenarcheota and Marsarchaeota. One last clade, composed of sequences identified as Prx1 and formed by organisms of the phylum Euryarchaeota, Crenarchaeota and Micrarchaeota, represents a divergent branch that does not have association with any of the previously mentioned Saccharolobus sequences. To study this in further detail multiple sequence alignment of all the sequences in this clade was performed. We identified the presence of the conserved motif PxDFTFVC [86], characteristic of the subfamily Prx1 (Figure 5b). Prx1 from the 2-cys subfamily has been reported to act with a dual enzymatic activity both as peroxiredoxin and catalase [89]. All organisms in which we identify this protein are either anaerobes or facultative anaerobes, suggesting that they could be utilizing this additional capability of Prx1 in specific conditions and ranges of H 2 O 2 as has been previously shown for organisms that have different peroxiredoxins [90]. In addition, all organisms with Prx1 have been isolated from similar environments (hotsprings or hydrothermal vents) which suggests an eco-niche adaptation that explains the narrow distribution of Prx1 in comparison with the other Prx conserved across a larger number of acidophilic Archaea. The unique divergence of this clade in the phylogenetic tree alongside the identification of the specific motifs confirms the idea that this may be a Prx1 type previously not reported in acidophilic Archaea, and further experimental validation is needed to confirm this finding.

Thioredoxin
The mechanism of action of thioredoxins (Trx, trxA) and thioredoxin reductases (TR, trxB) is similar to that of peroxiredoxins as they are needed to complete the redox cycle and reduce the disulfide bonds in Prx to continue with the peroxide elimination process. Our results show that most organisms with Trx have two copies, with both copies being part of the same orthogroup (paralogs) suggesting a shared evolutionary story that provides clues of a gene duplication event that could have occurred at an ancestral level. An interesting finding is a hypothetical protein (defined by the orthogroup OG0923) that appears contiguous to trxA, and was identified in members of Crenarchaeota, Euryarchaeota and Marsarchaeota with high conservation (Figure S5). Sequences from this orthogroup were analyzed by subcellular localization prediction tools which predict these proteins are in either the cellular membrane or extracellular space (Table S4). This highly conserved hypothetical protein could be important for these organisms as proteins in the cellular membrane of acidophiles are exposed directly to the physicochemical constraints of a low pH environment. Experimental analyses of this kind of protein of unknown function could help us reveal new proteins or mechanisms of oxidative response that are unique to Archaea in extreme environments. Predicted TR were also identified widely across acidophilic Archaea as part of the same orthogroup, with most members in Acidianus, Saccharolobus, Sulfolobus and Marsarchaeota having an extra copy ( Table S5). The presence of both Trx and TR in addition to different types of Prx suggests that the complete redox thiol cycle could play an important role in oxidative stress response in Archaea as has been previously demonstrated in Bacteria [91].

Methionine sulfoxide reductases
Another target sensitive to ROS oxidation is methionine, a sulfur-containing amino acid. The first stage of oxidation of this amino acid is the formation of methionine sulfoxide with two different stereospecific forms (Met-S-O and Met-R-O), process that could lead to the protein being aggregated and eliminated or even further oxidized to methionine sulfone (an irreversible state). MsrA can reduce Met-S-O and MsrB reduces the Met-R-O stereospecific form. The first protein, MsrA, was identified in several genera of the acidophilic Archaea. Previous studies have reported this to be mainly absent in hyperthermophiles [92]. In contrast, our results include hyperthermophiles such as Metallosphaera, Acidianus and Sulfolobus, all of which have optimal growth over 70ºC, pointing out a previous underestimation of the abundance of MsrA in Archaea. The second protein, MsrB, was only identified in Thermoplasmatales archaeon and Cuniculiplasma in Euryarchaeota, Marsarchaeota and Thaumarchaeota, completely absent from Crenarchaeota. The protein fRMSR also can reduce the same isoform as MsrB, but in free amino acids (not as part of a protein) and is characterized by a GAF-like domain, identified in Picrophilus and Thermoplasma, both organisms that also have MsrB. All these proteins share similarity with the fRMSR from Thermoplasma acidophilum, where structural and biochemical analysis has been performed to confirm its functional activity [93]. Given the high degree of similarity of the sequences, we suggest that all other copies in the same orthogroup could also be functionally active. Finally, both stereospecific forms can also be reduced by the molybdopterin-dependent sulfite oxidase family MrsP. We identified this protein in acidophiles from all phyla. Some organisms from Crenarchaeota (Saccharolobus, Metallosphaera, Vulcanisaeta and Thermocladium) have more than one copy, which could act as counterbalance to the fact they lack the other common mechanism for reducing Met-R-O. These results contrast with previous observations that found MsrP to be mostly absent in Archaea [92]. MsrQ is commonly used as an electron donor associated with MsrP in Bacteria [94], however, it was not identified in any of the organisms under study, suggesting that an alternative pathway to obtain electrons could be in use by Archaea. This may have special relevance when considering the differences in cell wall components compared with Bacteria. An overview of the described system is shown in Figure 5.

General Discussion
Our study provides valuable clues for the prediction of the mechanisms used by acidophilic Archaea to cope with oxidative stress and represents one of the most extensive comparative genomics analyses on this topic to date. An overview of all the analyzed mechanisms identified is presented in Figure 6. A comparison of the predicted proteins with those from nonacidophilic Archaea, indicates that most of the mechanisms are conserved widely across the different acidophilic and no-acidophilic phyla in this study. However, there are some interesting deviations from this observation. Superoxide dismutase (SOD) constitutes an interesting case, where all the structurally described Fe-SOD are from acidophilic Archaea, while other non-acidophilic Archaea have either Mn-SOD or the cambialistic Mn/Fe-SOD. This difference could be associated with acidophiles maintaining the ancestral version of Fe-SOD [7], given the high availability of iron in acidic environments even in current times, while non-acidophilic Archaea prefer SOD that could use other elements as catalytic centers [95]. Peroxiredoxins are also predicted both in acidophilic and non-acidophilic Archaea, but one specific family identified in this study (Prx1) has not been previously described in the literature in acidophilic Archaea, which highlights the importance of the use of comparative genomics strategies to help unravel the distribution of different proteins across a big dataset of organisms. Lastly, a hypothetical protein described by OG0923 was identified to be highly conserved in the genomic context of Trx and to be restricted to acidophiles. OG0923 is predicted to be transported to the membrane. We hypothesize that it carries out a role in acid stress as membrane proteins are directly confronted by the extremely low pH of the environment. This serves as an example of the use of bioinformatics in predicting potential new proteins associated with the oxidative stress response that could be particular to Archaea.
A notable difference that we detect between Archaea and Bacteria is that both aerobic and anaerobic representatives of acidophilic Archaea use SOD to deal with superoxide whereas anaerobic Bacteria use SOR [96]. The Archaeal preference for SOD could be related to the issue of the increased availability of Fe in acidophilic environments.
It has been observed that HGT events play an important role in adaptation and speciation in Archaea, and that many of these events imply the acquisition of genes from Bacteria (mostly acidophilic Bacteria) to Archaea [97,98]. The high similarity between the sequences from Archaea and acidophilic Bacteria also points out that the HGT event could be associated with specific niche adaptations that confer ecological advantage to these organisms. The importance of HGT events in acidophiles has been previously highlighted in several works [99][100][101]. Cross domain events of HGT were also predicted and could play a role in the evolution of catalases and Dps. In both cases of HGT, we identify a restricted presence in specific lineages of Archaea in contrast with what is found in Bacteria, where both proteins are widely distributed across the domain [102,103], which makes us suggest that the most probable case for direction of this HGT event is from Bacteria to Archaea.
The association of Archaea and early anoxic Earth environments and their deep rooting in the tree of life is commonly discussed in literature in the context of Archaea as candidates to harbor ancestral mechanisms and the possibility that they may illuminate the origin of different traits [104,105]. The wide conservation across organisms of superoxide dismutase, peroxiredoxin, and rubrerythrin in both anaerobic and aerobic organisms, could suggest an ancestral origin of these proteins [82,106,107], and that they were evolved in similar environments to the early LUCA in response to local concentrations of oxygen before the GOE [108]. It has been previously reported that multiple mechanisms for oxidative stress response are found in hyperthermophiles [21] and we identified a similar pattern in acidophiles with no major differences across lifestyles (such as aerobic versus anaerobic or thermophiles versus mesophiles). Instead, the distribution of oxidative stress response mechanisms present in each organism seems to be largely determined by phylogeny. Interestingly, the same type of distribution was observed in an extensive study of the stress enzymes in Bacteria [11], showing that this could be a trend globally associated with oxidative stress response. Figure 6. Summary of predicted oxidative stress response proteins in acidophilic Archaea. The heatmap depicts the presence of each protein as a percentage of the genomes it was identified in (see insert scale-bar). Genomes were grouped at the genus level. Colors correspond to phylum level with yellow for Crenarchaeota, red for Thaumarchaeota, green for Marsarchaeota and blue for Euryarchaeota.

Conclusions
Through comparative genomics of more than 200 genomes from acidophilic Archaea, we identified a wide range of mechanisms to cope with oxidative stress mostly conserved across all 5 different phyla showing no clear difference between aerobic and anaerobic lifestyles. Nonetheless, with this study we could unravel previously unidentified protein families, missed in individual studies, as the Prx1 family present in several acidophiles. On the other hand, the conservation of two proteins with iron coordinating centers, SOD and Rbr, across almost all organisms in study, suggests an early evolutionary origin that could have arisen in an anaerobic early Earth environments and be maintained in acidic environments. Furthermore, HGT events from Bacteria to Archaea, exemplified in the case of catalases and Dps highlight the importance of these events in prokaryotic adaptation to oxidative stress.

Supplementary Materials
The following are available online Figure S1: Multiple sequence alignment of Crenarchaeota Fe-SOD, Figure S2: Unrooted phylogenetic tree constructed from the predicted amino acid sequences of SorA in Aciduliprofundum and Archaeal sequences from sorGOdb, Figure S3: Unrooted phylogenetic tree constructed from the predicted amino acid sequences of KatG in Aciduliprofundum, Cuniculiplasma and Thermoplasmatales and their best hits in the NCBI, Figure S4: Unrooted phylogenetic tree constructed from the predicted amino acid sequences of OsmC in the different acidophilic Archaea, Figure S5: Genomic context of trxA and the conserved hypothetical protein OG0923 in several representative genomes, Table S1: Genomic features for acidophilic genomes in the study, Table S2: Functional annotations used for the identification of each protein associated with oxidative stress response, Table S3: Species clustering analysis for each of the genomes in the study done by Fast-Ani, Table S4: Subcellular location prediction results for the hypothetical protein OG0923 and MsrP, Table S5: Orthogroups for all the identified oxidative stress response associated protein with the corresponding protein id for each of the genomes in the study.

Author Contributions
GN designed the research. GN and EV performed the research. GN, EV, DC and DH analyzed the data. GN and DH wrote the paper. All authors participated in the construction of the final manuscript. All authors read and approved the final manuscript.

Institutional Review Board Statement
Not applicable.

Informed Consent Statement
Not applicable.

Data Availability Statement
The data presented in this study are available in article.

Conflicts of Interest
The authors declare no conflict of interest.