Two datasets were built, one including bacterial and archaeal FtsZs (BA reconstruction) and the other including only the bacterial sequences (B reconstruction). After removing redundancy by taxonomic classes and spurious sequences, the final multiple sequence alignments (MSAs) were treated with three different trimming methods, referred to as Manual-, BMGE-, and trimAl-gt70-trimming (refer to section “Materials and methods”). Subsequently, we performed phylogenetic reconstructions using empirical models automatically selected by Modelfinder (Kalyaanamoorthy et al., 2017) and complex models (LG4X and C20 + R + F). We used an ultrafast bootstrap for manual-trimmed reconstructions and nearest neighbor interchange (NNI) search analyses (Minh et al., 2020) for all the MSAs.

In the resulting reconstructions, most bacterial phyla show a monophyletic pattern, although the topology was sensitive to different parameters, such as taxonomic sampling, sequence redundancy threshold, and phylogenetic methods (Supplementary Figure 1). Basal nodes were weakly supported (<90%) in most reconstructions, except for UFBoot2 reconstructions. This low support in basal nodes can be explained by the antiquity and short length of these proteins (MSAs of 350 amino acids applying inclusive trimming). Therefore, despite the fact that UFBoot2 values lower than 95% should not be considered reliable nodes (Hoang et al., 2018), we interpreted these reconstructions with a relaxed threshold of confidence as the biological relationship between species is limited when collapsing nodes with 95% UFBoot2 values (Supplementary Figure 2). On the contrary, collapsing at 85% resulted in biologically meaningful clades (Figure 2 and refer to Supplementary Figures 2, 3 for an extended version).

Taking the Manual + UFBoot2 + MF reconstruction as a reference, the phylogeny was better resolved in the B reconstruction than in the BA reconstruction since the basal nodes better resolved the relationships between bacterial phyla (specifically in Gracilicutes, Supplementary Figure 2). This suggests that the addition of archaeal sequences promotes long branch attraction (LBA) artifacts, providing weaker topologies. This was observed in our reconstructions where archaeal FtsZ sequences attracted different bacterial FtsZ using various trimming methods and evolutionary models (refer to BA reconstructions in Supplementary Figure 1). Thus, B reconstruction was selected as a reference (Figure 2).

The bacterial FtsZ tree recovered the two main bacterial supergroups, Terrabacteria and Gracilicutes (Figure 2; Coleman et al., 2021). Between these two supergroups, intermediary branches are found, including Elusimicrobiota, Spirochaetota, Fusobacteriota, Synergistota, Armatimonadota, and Deinococcota phyla, or other phyla such as Firmicutes and Bacteroidota, forming paraphyletic groups and showing unstable locations in our reconstructions (Supplementary Figure 1). The FtsZ of Candidate Phyla Radiation (CPR, denoted as Patescibacteria in this GTDB version) is phylogenetically associated with Chloroflexota, consistent with their proposed origin (Coleman et al., 2021); however, this relationship was only observed in Manual-trimming reconstructions. Some CPR genomes also contain a divergent copy of FtsZ, which usually branches with Verrucomicrobiota, which are the most divergent groups (Figure 2 and Supplementary Figure 1). The verrucomicrobial FtsZ and the second copy of CPR had different topologies in the BA and B reconstructions (Supplementary Figures 1, 2), demonstrating their phylogenetic instability and suggesting strong divergence and artifacts usually associated with it, such as LBA. The divergent CPR FtsZ (together with that of Verrucomicrobiota) is nested in the canonical CPR FtsZ group in the BMGE-B-bnni reconstruction (Supplementary Figure 1), suggesting that it originated from the original FtsZ of CPR. In contrast, in the B reconstruction (refer to Figure 2 for UFBoot2 and Supplementary Figure 1 for NNI analysis), Verrucomicrobiota branches with other Gracilicutes near to Bacteroidetes (and other close relatives such as Cloacimonadota and WOR-3), consistent with their species relationship and supporting the idea that reconstructions without archaeal sequences solve more adequately bacterial evolution. Therefore, despite the phylogenetic instabilities or strong sequence divergence, the separation of Gracilicutes and Terrabacteria in most reconstructions (Supplementary Figures 1, 2 and in Manual-B-UFboot2 reconstruction supported by 95%; Figure 2) strongly suggests that FtsZ was present in the last bacterial common ancestor (LBCA).

We then inspected the genomic association of all these prokaryotic FtsZs (Figure 3). The most conserved Pfam domains present in the genes surrounding bacterial ftsZs are related to cell division, for example, FtsA, FtsW, or FtsQ, or peptidoglycan synthesis, for example, murein ligase, glycosyl transport, MreB, or MraY (Figure 3A), which form the dcw cluster. This genome context is well conserved in Gracilicutes, Terrabacteria, and CPR, showing that the dcw gene cluster, was already established in the LBCA. In the few bacteria having two copies of the ftsZ gene, such as Bipolaricuta, the additional copy is located outside the dcw cluster (some of them in conserved loci), indicating possible subfunctionalization due to gene duplication.

The genomes of most archaeal species present two copies of ftsZ, coding for FtsZ1 and FtsZ2, with the exception of Crenarchaeota archaea that mostly retained FtsZ2, or present divergent FtsZs branching outside of the canonical archaeal FtsZ groups (Figure 4 and Supplementary Figure 4, for an extended view). Alternatively, other Crenarchaeota classes do not present any FtsZ, and in the case of a few Euryarchaeota and a few Asgardarchaeota, only one copy of the gene was found, as previously reported (Pende et al., 2021).

The Crenarchaeota and Asgardarchaeota sequences are closely related (Figure 4), in agreement with their species relationships. The Euryarchaeota sequences (according to GTDB) tend to form monophyletic groups, but they are not monophyletic with other close relatives such as Halobacterota and Thermoplasmatota. However, Hydrothermarchaeota and Hadarchaeota sequences are closely related to Euryarchaeota (Figure 4), which is in agreement with the species’ relationships. The DPANN supergroup branches paraphyletically: Micrarchaeota and Altiarchaeota branch close to the Euryarchaeota supergroup but distally to other DPANN such as Nanoarchaeota and Huberarchaeaota that show different topologies of the FtsZ1 and FtsZ2 branches.

Acknowledging the low resolution in some basal nodes, and as previously suggested (Pende et al., 2021), the monophyletic branching of most archaeal phyla, combined with topological similarities with the archaeal species tree, suggests that FtsZ1 and FtsZ2 are derived from an ancestral duplication in the last archaeal common ancestor (LACA). In addition, FtsZ evolution in archaea is irregular as some archaeal groups have lost one or both FtsZ copies, for example, Crenarchaeota. However, others might have regained FtsZ by lateral gene transfer (LGT) from other archaea, as shown by suspicious phylogenetic groups in FtsZ2, such as Micrarchaeota and Altiarchaeota branching within Euryarcheaota, or the monophyletic branching of Crenarchaeota and Halobacterota sequences (refer to question marks in Figure 4), both of which could represent phylogenetic artifacts (assuming that FtsZ has been inherited vertically between the species) or actual LGT events (as suggested by the latter case with full UFBoot2 support). Therefore, although vertical inheritance seems to be the main mode of evolution of FtsZ in archaea, LGT events between them should not be discarded. However, it is important to remark that we did not detect notorious LGT of ftsZ genes between archaea and bacteria.

The genomic contexts of the archaeal ftsZ1 and ftsZ2 genes are different from those of their bacterial counterparts, but they are also different from each other (Figure 3B). The gene ftsZ1 has a conserved genome context, whereas the one of the ftsZ2 is not as well conserved (Figure 3B). In contrast to what is observed in bacteria, the domains associated with genes found in the genome context of archaeal ftsZ are related to informational housekeeping genes, mainly related to protein biosynthesis, including ribosomal and t-RNA genes (especially for the genome context of ftsZ1 but also in Crenarchaeota ftsZ2; Figure 3B). The genomic contexts of ftsZ2 usually contain a gene bearing the RHH_1 Pfam domain, which probably functions as a transcriptional regulator of the CopG family (Acebo et al., 1998), suggesting that the transcriptional expression of the ftsZ2 gene cluster could trigger a specific transcriptional regulation response.

Archaeal genomes often reflect functional expansions of the FtsZ/tubulin proteins (refer to arcs in Figure 5) as exemplified by the well-known CetZ and others, such as artubulins (Yutin and Koonin, 2012) and OdinTubulin from Odinarchaeota (Zaremba-Niedzwiedzka et al., 2017). In addition, we also identified other archaeal paralogs found in Asgardarchaeota (Asg.Tub) and another one in Halobacterota (Halo.Tub); although the phylogenetic positions of both were unstable, branching with tubulins (Figure 4) or with CetZ sequences (Figure 6; further discussed later).

CetZ was likely present in the Euryarchaeota ancestor as it is found in most classes, such as Archaeoglobi, Methanomicrobia, Halobacteria, Syntrophoarchaeia, and Hydrothermarchaeota (Figure 5). We found other CetZ-like proteins, including sequences from Halobacterota, Altiarchaeaota (DPANN), and sparse paraphyletic groups of Chloroflexota–Dehalococcoidia and Firmicutes–Bacilli/Clostridia (the latter was identified as TubZ; Figure 4).

As no recognizable domains were detected in the C-termini of the CetZ sequences (Figure 1), we attempted to clarify their evolutionary relationships by constructing protein models (hidden Markov model, HMM) of the C-termini of these proteins and performing hmmsearches against all our tubulin/FtsZ protein family collection (Supplementary Data). This CetZ_C HMM showed significant sequence similarities (e-value approximately 1e^–10) with the other CetZ-like proteins (including bacterial Chloroflexota but not Firmicutes TubZ sequences) but also weaker similarities (e-value approximately 1e^–3) with FtsZ, in particular with the archaeal sequences (see external green heatmap in Figure 5). CetZ_C HMM does not match with tubulins or archaeal tubulin-like proteins nor bacterial FtsZ. This possibility is more evident when visualizing the same tree unrooted (inset Figure 5). In this view, the closer proximity between the archaeal FtsZ sequences and tubulin-CetZ sequences is observed. Therefore, this phylogenetic signal, together with the homology at C-terminus, suggests that CetZ originated from the archaeal FtsZ.

The C-terminus of TubZ is substantially divergent, even between TubZ sequences. The HMM of CetZ_C did not match with that of Firmicutes TubZs, but the C-terminus HMM of TubZ matched with CetZ and archaeal FtsZ sequences (refer to external purple heatmap in Figure 5, e-value in the range of 1e^–2), providing consistent support for the evolutionary relationship observed in the phylogenetic reconstructions (Figure 4). Thus, considering the potential antiquity of archaeal FtsZ compared with TubZ, it is likely that TubZ originated from CetZ or archaeal FtsZ sequences.

Next, we inspected the C-termini of the archaeal groups branching basally to the tubulins, one of which is formed by Asgardarchaeota (Asg.Tub), another one by Halobacterota (Halo.Tub), and the third one by artubulins (Figure 4 and Supplementary Figure 4). The intermediary positions of these groups between CetZ and eukaryotic tubulins were stable when different phylogenetic methods were applied (Supplementary Figure 5). The monophyletic branching of Heimdallarchaeia and Lokiarchaeia is remarkable since it could illustrate the deep evolution of these paralogs in Asgardarchaeota. Here, it is important to mention that the scarcity of these groups translates into strict models that do not contain divergent sequences (note that the sequence conservation at the C-terminus is usually weak). The searches of their respective C-termini HMMs did not match with any other tubulin/FtsZ protein (refer to external red, black, and blue heatmaps in Figure 5 and Supplementary Data), even with each other. However, by inspecting the MSA of the selected sequences, we observed that their C-termini show limited similarity to CetZ sequences (Supplementary Figure 6). Conversely, we found that artubulins are more closely related to eukaryotic tubulins than to Halo.Tub and Asg.Tub since these contain important indels that characterize archaeal CetZ (and CetZ-like) sequences from eukaryotic tubulins (black highlight in Supplementary Figure 6), which is also supported by the mid-point rooting in the respective phylogeny (Figure 6). However, it is important to note that some of these CetZ-like sequences presented few conserved blocks at the N-terminus that resemble eukaryotic sequences (particularly in Altiarchaeota and different Halobacterota classes such as Methanosarcina, and Halobacteria, Halo.Tub, 5). Therefore, although our analyses suggest that Halo.Tub and Asg.Tub seems to be phylogenetically associated with the tubulins (Figure 4); inspections of the MSA and mid-point rooting of the selected set of sequences (Figure 6) argue against an apparent tubulin-like affiliation.

The genome context of these archaeal “tubulins like” is characterized by unknown genes that we manually annotated with HHPRED (Söding et al., 2005). We observed that most of these contain coiled-coil regions with different functions (refer to gray triangles in Figure 6). In cases such as Halo.Tub, their genome context contains several copies in a quasi-tandem disposition. In contrast, artubulin from Nitrososphaera korensis is associated with a potential Snf7 protein upstream, as previously observed (Yutin and Koonin, 2012). This abundance of unknown and unrelated genes could illustrate the diversity of molecular systems and functions in which these archaeal CetZ/Tubulin homologs are involved.

Therefore, the divergence of their C-termini, together with the sparse distribution of each phylogenetic group, shows the phylogenetic positions of these Halo.Tub and Asg.Tub in relation to each other and with eukaryotic tubulins must be interpreted cautiously. These proteins appear to be more closely related to CetZ sequences than eukaryotic tubulins. The genomic context between the CetZs and these CetZ-like sequences is not related (Figure 6), hindering inferences regarding the origin of these proteins. Indeed, the idea of LGTs between archaea, followed by neofunctionalization, should not be discarded for these specific proteins. Nevertheless, our results suggest that artubulins are the only ones of these proteins that could be related to ancestral eukaryotic tubulins, as previously pointed out (Yutin and Koonin, 2012).

The monophyletic branching of eukaryotic tubulins (Figure 6) confirms that eukaryotic tubulins had a single origin and suggests that its paralogs (α, β, γ, δ, ε, and ζ) were already established in the LECA as previously proposed (Findeisen et al., 2014). Other divergent tubulin-related proteins such as Misato, were not included in this analysis because of their strong sequence divergence. Tubulins have a very restricted distribution in prokaryotes and are only found in one bacterial species and one archaeal metagenome. OdinTubulin is only found in Odinarchaeota from the Asgardarchaeota phylum and a pair of tubulins, BtubA and BtubB, is present in various Prostecobacter species from the Verrucomicrobiota phylum. These prokaryotic sequences branch basally and paraphyletically to the α- and β-tubulins, although BtubA and BtubB can also branch basally to α- and β-tubulins, respectively (Figures 4, 6). These phylogenetic instabilities of BtubA/B have been observed previously (Findeisen et al., 2014). In our reconstructions, their monophyletic branching was exclusively observed in the inclusive trimming, removing only gap positions (Supplementary Figure 5 right panel). In contrast, using a stricter trimming with BMGE, BtubA and BtubB were paraphyletic (Supplementary Figure 5 right panel; note that the topology of eukaryotic tubulins also differed and branched monophyletically to the β/γ group in BMGE trimming). In contrast, using an inclusive trimming (trimAl -gt 0.10) and considering only tubulin and CetZ(–like) sequences, the paraphyly of BtubA/B was again recovered (Figure 6). Therefore, the mono/paraphyly of BtubA/B seems to be sensitive to the trimming methods and the diversity of sequences included in the MSA. Similarly, the locations of the OdinTubulin in our reconstructions (Figures 4, 6) differed from those obtained in previous studies, in which the Odinarchaeota sequence branched basally to the whole eukaryotic tubulin family (Zaremba-Niedzwiedzka et al., 2017), instead of basally to α- and β-tubulins. These differences could be due to different parameters of the reconstructions (taxonomic sampling, the use of different outgroups, or different methods for the MSA trimming), showing the phylogenetic instability of these prokaryotic tubulins and, therefore, the difficulty in inferring their evolutionary placement. Nevertheless, from our reconstructions, it is clear that OdinTubulins and BtubA/B are closer to α- and β-tubulins, in agreement with recent analyses (Akıl et al., 2022).

In eukaryotes, the α- and β-tubulins generally form 13 subunits filaments as imposed by the nucleation of the γ-tubulin ring composed of 13 γ-tubulins (Chaaban and Brouhard, 2017; Cavalier-Smith and Chao, 2020). The Prosthecobacter’s BtubA/B microtubules are formed by four or five protofilaments, presenting heterodimer polarity which can polymerize without the need of the microtubule processing centers, tubulin γ, never found in prokaryotes. In addition, btubA/B genes form an operon together with the btubC gene, a “bacterial kinesin light chain” that has been shown to stabilize the filament (Deng et al., 2017). Therefore, the functional similarities between BtubA/B and α- and β-tubulin supports their close relationships.

In contrast to the dual BtubA/B, there is a unique OdinTubulin in the Odinarchaeota metagenome sequenced from hydrothermal environments (Zaremba-Niedzwiedzka et al., 2017). Despite the possible incompleteness of the Asgardarchaeota metagenomes, the sequencing of several Asgardarchaeota clades (Liu et al., 2021) confirmed that this OdinTubulin is still uniquely found in this archaeon. In vitro and at high temperatures, this OdinTubulin forms tubules with short curved protofilaments coiling around the tubule circumference, more similar to FtsZ, rather than running parallel to its length as in eukaryotic microtubules (Akıl et al., 2022). Nevertheless, the genome context of this OdinTubulin displays interesting features related to eukaryotes. A hypothetical protein showing similarity to Prc1 at the N-terminus and Snf8 at the C-terminus (containing coiled-coil segments) is located upstream of the OdinTubulin gene. Prc1 is a key regulator of cytokinesis that cross-links antiparallel microtubules, while Snf8 is a member of the ESCRT-II system. Downstream of the OdinTubulin, we found a gene containing coiled-coil segments weakly related to the ATG16 Pfam domain, which belongs to the eukaryotic autophagy system. However, according to the NCBI annotation, this gene is annotated as chromosome partitioning protein Smc. The two genes surrounding the OdinTubulin are specific to this genome as protein searches do not retrieve any homolog sequences in the three domains. Together, these observations suggest that this OdinTubulin is specific to species related to Odinarchaeota and that it could be functionally related to a eukaryotic-like system involved in endomembrane system and cell division (exemplified by the hypothetical prc1-snf8 gene), in agreement with its phylogenetic position.

Therefore, even though BtubA/B and OdinTubulin are closely related at the sequence level, as they are associated with α- and β-tubulins, they seem to be involved in different molecular systems or complexes, where the BtubA/B pair functions as eukaryotic-like microtubule filaments and OdinTubulin could function as monomeric filaments in combination with the Pcr1-Snf8 protein.

To collect homologous sequences from the tubulin/FtsZ protein family, we performed HMMSEARCH (HMMER 3.1b2) (e-value threshold 1e^–3; Prakash et al., 2017) against a local database comprising all the NCBI genomes taxonomically annotated with GTDB (release bac120 and arc122; Parks et al., 2018), plus 36 diverse eukaryotes (24,664 taxa in total). As a query, we used the Tubulin Pfam model (PF00091) spanning the N-terminus as it is the commonly shared domain of this family, for both FtsZ, tubulins, and derivatives. This search also detected FtsZ-like sequences (Makarova and Koonin, 2010), containing the Tubulin_2 Pfam model (PF13809) instead. The scattered and sparse taxonomic distribution of this FtsZ-like (Tubulin_2), together with their irregular domain architecture, that is, they do not have the FtsZ_C nor Tubulin_C Pfam domains but have important N/C-terminus extensions (in Planctomycetes a protein kinase domain is found at the N-terminus). Their sequence divergence (to tubulin/FtsZ family members but also between them) suggested that these proteins had a different evolutionary history to the canonical tubulin/FtsZ protein family (Makarova and Koonin, 2010). Thus, to obtain more accurate multiple sequence alignments (MSAs) and phylogenies, we exclude these sequences from our main reconstructions.

To define the main tubulin/FtsZ subfamilies, we aligned all the sequences (45,929 sequences) with MAFFT v7.310 (Katoh and Standley, 2013), applied gap trimming using trimAl 1.4.22 (–gt 0.7; Capella-Gutiérrez et al., 2009), removed sequences with lower alignment coverage than 60%, built a Fasttree 2.1.10 (Price et al., 2010), and selected the branches of interest to build different datasets (FtsZ, tubulin-like proteins or to exclude FtsZ-like sequences). To avoid the over-representation of some taxa and build taxonomically balanced datasets, we removed redundant sequences (using CD-HIT; Fu et al., 2012) by taxonomic classes (following the GTDB taxonomy) and applied different cut-offs depending on the number of sequences (from 20 sequences, 95% identity, up to >250, 55%).

From these reduced datasets, we built different datasets for specific purposes: that is, focused on prokaryotic FtsZ (Figure 2), all tubulin/FtsZ and homologous proteins with a special focus on archaea (Figures 4, 5), and tubulin/CetZ-like proteins (Figure 6). All trees were visualized and annotated using iTOL (Letunic and Bork, 2019). Raw data for the MSAs and respective phylogenetic trees are provided in Supplementary Data.

For the phylogenies of bacterial FtsZs (Figure 2), we randomly selected one sequence from each taxonomic order from the FtsZ prokaryotic dataset (note that some taxa such as Verrucomicrobiota, CPR were manually increased by selecting representative sequences). Then, we removed spurious sequences based on the visualization of the MSA (trimAl -gt 20) and the respective phylogenetic tree (Fasttree). The final set of sequences was aligned with MAFFT-linsi, and three filtering methods were applied: gap-trimmed alignment (–gt 0.2), BMGE-1.12 (–h 0.55 -m BLOSUM30), using trimAl (–gt 0.7). Due to some miss-aligned regions at the C-terminus in FtsZ from Verrucomicrobia, the manually trimmed alignment was re-aligned as follows: We ran MAFFT-Linsi excluding CPR-Verrucomicrobia divergent sequences, and these later ones were aligned independently; then, both MSAs were treated with trimAl -gt 0.2 and merged with MAFFT (–merge option). Each trimming method had two sub-datasets: one including bacterial and archaeal FtsZs (BA), and the other one containing exclusively bacterial FtsZs (B). From the manually trimmed alignments, B MSA was obtained by removing archaeal sequences from BA MSA, while in BMGE- and trimAl-gt70-trimming MSAs, B, and BA MSAs were aligned independently through MAFFT-Linsi.

For the phylogeny including FtsZ, tubulins, and others (Figures 4, 5), the dataset was constructed by a selection of representative bacterial FtsZs and eukaryotic tubulins, plus an archaeal selection based on redundancy filtering and recovering the discarded paralogs from the remaining taxa. Then, this set of the sequence was aligned with MAFFT-Linsi and trimmed with trimAl (–gt 0.2) and BMGE (the same options as mentioned earlier).

For the phylogeny of tubulins and CetZ(–like) sequences (Figure 6), we manually selected a set of representative sequences branching basally to eukaryotic tubulins extracted from the tree of Figure 3.

Phylogenies were inferred using maximum-likelihood methods with IQ-TREE (Nguyen et al., 2015). We obtained branch support with ultrafast bootstrap (UFBoot2, option -B 1000; Hoang et al., 2018), and the evolutionary models of each set of the sequences were automatically selected using ModelFinder (Kalyaanamoorthy et al., 2017), following the BIC criterion (–m MFP). In addition, we also run UFBoot2 + NNI bootstraps (-B 1000 -bnni options), testing complex evolutionary models such as LG4X (–m LG4X) or C10-30 (–mset LG -madd LG + C10, LG + C20, LG + C30, LG + C10 + R + F, LG + C20 + R + F, LG + C30 + R + F –score-diff all).

To analyze the genome contexts of the tubulin/FtsZ protein family, the DNA sequence corresponding to 10 Kb surrounding (5 kb up/downstream) each of the respective coding genes was extracted. In that fragment, genes were predicted using Prodigal (Hyatt et al., 2010), with default parameters, and the predicted proteins annotated with HMMSCAN (Potter et al., 2018) on the Pfam database (Finn et al., 2011). The abundance of Pfam domains in these proteins was then quantified, counting only the presence and not the repetitions inside the same protein or genome context, that is, multiple copies of the same domain in one gene or multiple gene copies in one genome context was only counted once.

To build the protein models (hidden Markov model, HMM) of the C-termini of the different proteins of interest (CetZ, TubZ, artubulins, Halo.Tub, and Asg.Tub), the sequences were aligned individually, and soft gap trimming was applied with trimAl (–gt 0.1). Then, the C-termini fragments were extracted starting from the end of the Tubulin Pfam domain (N-terminus domain). Protein models were made using HMMBUILD (Potter et al., 2018) and were then used to perform hmmsearches against our dataset of the FtsZ/tubulin protein family without an e-value threshold. Raw data for the construction of protein models and respective searches are provided in Supplementary Data.