Phylogenetic study of Class Armophorea (Alveolata, Ciliophora) based on 18S-rDNA data

The 18S rDNA phylogeny of Class Armophorea, a group of anaerobic ciliates, is proposed based on an analysis of 44 sequences (out of 195) retrieved from the NCBI/GenBank database. Emphasis was placed on the use of two nucleotide alignment criteria that involved variation in the gap-opening and gap-extension parameters and the use of rRNA secondary structure to orientate multiple-alignment. A sensitivity analysis of 76 data sets was run to assess the effect of variations in indel parameters on tree topologies. Bayesian inference, maximum likelihood and maximum parsimony phylogenetic analyses were used to explore how different analytic frameworks influenced the resulting hypotheses. A sensitivity analysis revealed that the relationships among higher taxa of the Intramacronucleata were dependent upon how indels were determined during multiple-alignment of nucleotides. The phylogenetic analyses rejected the monophyly of the Armophorea most of the time and consistently indicated that the Metopidae and Nyctotheridae were related to the Litostomatea. There was no consensus on the placement of the Caenomorphidae, which could be a sister group of the Metopidae + Nyctorheridae, or could have diverged at the base of the Spirotrichea branch or the Intramacronucleata tree.


Sequence acquisition
Since the monophyly of the Armophorea is questionable and the phylogenetic affinities are variable (Shin et al., 2000;Miao et al., 2009a,b;Vd'acny et al., 2010), we broadly sampled ciliophoran 18S sequences to include representatives of all recognized ciliate classes (Lynn, 2008). One hundred and ninety-five ciliate 18S rDNA sequences, including 44 armophorean sequences, were downloaded from the NCBI database (Table S1) and assembled for multiple alignment. To avoid confusion when discussing the systematics of armophoreans, only sequences that were identified at least to the level of family or order were sampled. The armophorean sequences used represented the Families Nyctotheridae Affa 'a, 1987(Order Clevelandellida Puytorac & Grain, 1976), Caenomorphidae Poche, 1913and Metopidae Kahl, 1927(Order Armophorida Jankowski, 1964. Many of the armorphorean 18S sequences available from the NCBI were partial so that the missing homologous regions were treated as absent.

Sequence alignments
There is generally little agreement on how to treat 'ambiguously alignable' regions. Some authors recommend elimination of the nucleotide positions of ambiguous alignments as a means of improving phylogenetic hypotheses (Olsen and Woese, 1993;Swofford et al., 1996;Talavera and Castresana, 2007), while others consider that such positions contain information that is potentially useful for phylogenetic reconstructions (Lutzoni et al., 2000;Aagesen, 2004;Redelings and Suchard, 2009). In this study, we opted to preserve this information and to explore different alignments (Wheeler, 1995;Doyle and Davis, 1998). To assess how different alignment criteria might influence phylogenetic hypotheses obtained from the ciliate 18S data, the sequences were multiple-aligned using the ClustalW algorithm and the 18S rRNA secondary structure. The resulting alignment files were inspected in BioEdit v7.0.5 (Hall, 1999) to code leading and trailing gaps as missing data. The overall and mean p-distances displayed in Tables 1 and 2 were calculated with the program MEGA   572 Paiva et al.  (Tamura et al., 2011), using pairwise deletion as a treatment for gaps and missing data.

ClustalW alignment
The sequences were aligned with the program ClustalW 1.81 through the CIPRES Science Gateway (Miller et al., 2010). Gap-opening (GOP) and gap-extension (GEP) penalties for multiple-alignment were optimized based on the averaged Q-scores, calculated with the program TuneClustalX (Hall, 2004) and used as a benchmark for alignment accuracy (Hall, 2005). The range explored for the parameters included GOP values of 10, 20, 30 and 40, each combined with integer GEP values that varied evenly from 1 to GOP/2. For comparison, we also included the ClustalW default values, i.e., GOP = 10 and GEP = 0.2, to yield 76 combinations ( Figure 1). The alignment of the highest averaged Q-score was inspected and refined by eye in BioEdit as a means of improving the averaged Q-score. This alignment is referred to as the Q-score optimized alignment (QOA).

Secondary structure alignment
As an alternative approach, the sequences were aligned based on the eukaryotic 18S rRNA secondary structure using the SINA web aligner (Pruesse et al., 2007) with its default settings. The alignment was inspected in BioEdit to remove gap-only columns followed by further refinement by eye that took into account the structural similarity among the sequences. This alignment is referred to as the secondary structure alignment (SSA).

Phylogenetic analyses Sensitivity analysis
To evaluate whether differences in the GOP-GEP choices in the ClustalW alignment affected the hypotheses for ciliate 18S phylogeny, the BioNJ algorithm (Gascuel, 1997) implemented in the program PAUP* 4b10 (Swofford, 2003) was used to build neighbor-joining (NJ) trees from p-distance matrices of each alignment. For this purpose, the alignments were not refined manually in order to prevent altering the decisions made by ClustalW for each GOP-GEP combination. The resulting trees were gathered with PAUP* and a strict consensus tree was built to show the insensitive branches, with emphasis on ciliate higher taxa ( Figure 2).

Analyses of the QOA and SSA data sets
Both data sets were independently analyzed using Bayesian inference (BI), maximum likelihood (ML) and maximum parsimony (MP) frameworks. In all resulting trees, the root was placed a posteriori at the Intramacronucleata-Postciliodesmatophora split (Lynn, 2008). The BI and ML routines employed a GTR + I (= 0.19) + G (= 0.6) model, selected based on the Akaike information criterion (AIC) (Akaike, 1974;Bos and Posada, 2005) in MODELTEST 3.7 (Posada and Crandall, 1998). To test whether the optimal trees were significantly differ-Expanded phylogeny of Armophorea 573  ent from suboptimal ones, all of the trees were statistically compared based on the data sets and optimality criteria by which they had been generated, with emphasis on the monophyly/non-monophyly of the Armophorea; this comparison was done using the approximately unbiased (AU) test for maximum likelihood and the Templeton test for maximum parsimony (Templeton, 1983;Shimodaira, 2002). The tests were done using the package CONSEL v0.1i (Shimodaira and Hasegawa, 2001;Shimodaira, 2004) and PAUP* 4b10, respectively. The taxonomy of higher taxa displayed in the phylogenetic trees (Figures 2-6) is mostly according to Lynn (2008), although the taxonomy of the Ventrata and Lamellicorticata agrees with Vd'acny et al. (2010).
Bayesian inference was determined with MrBayes 3.1.2 (Ronquist and Huelsenbeck, 2003) and was based on two independent Markov chain Monte Carlo (MCMC) simulations that were run with four chains of 1,000,000 generations. The trees were sampled every 200 generations (temperature of heat chains = 0.2), and the first 100,000 574 Paiva et al. generations were discarded as burn-in. A 50% majority rule consensus of the trees remaining after burn-in was used to calculate the Bayesian posterior probabilities of the recovered kinships; these probabilities were used as node support measures for BI (Schneider, 2007).
For ML, the data sets were analyzed with the program PhyML 3.0 (Guindon and Gascuel, 2003), starting with a BioNJ tree for which the likelihood was improved via SPR branch-swapping to generate the ML tree. Node stability was evaluated via the SH-like aLRT branch support Expanded phylogeny of Armophorea 575 Figure 4 -BI/ML tree hypothesized from the secondary structure alignment (-lnL = 68901.83179). Arrows indicate sequences that supposedly belong to metopids; values associated with nodes are Bayesian posterior probabilities/aLRT support; * = full support; -= support < 50%. Scale bar = 2 substitutions per 10 nucleotide positions. (Anisimova and Gascuel, 2006;Guindon et al., 2010). The MP analyses were done with the program TNT 1.1 (Goloboff et al., 2003(Goloboff et al., , 2008, using a strategy that combined routines of parsimony-ratchet (Nixon, 1999) with tree-drifting, tree-fusing and sectorial searches (Goloboff, 1999) in order to find optimal cladograms. Gaps were coded as a "fifth base" to accommodate their phylogenetic information (Giribet and Wheeler, 1999;Schneider, 2007), and only parsimony-informative characters were analyzed.
We agree with Goloboff (1993) and Platnick et al. (1991Platnick et al. ( , 1996 concerning the incompleteness of equal-weighted cladistics and thus applied the implied-weighting approach of Goloboff (1993) to the data. For this, Goloboff's parameter K was explored at integer intervals varying evenly from 1 to 10. The resulting trees were summarized via strict consensus. The common synapomorphies of the nodes of interest were assessed by optimizing the characters of each data set onto the trees using TNT. Node support was mea- 576 Paiva et al. sured via 1,000 symmetric resampling (Goloboff et al., 2003(Goloboff et al., , 2008 replicates, assuming K = 10.

Sequences and alignments
The global optimum for the averaged Q-scores distribution was found at GOP = 30, GEP = 14 (Figure 1), and the corresponding alignment was selected as optimal among the GOP-GEP combinations that were explored. However, regions with higher averaged Q-scores may exist for wider GOP-GEP intervals. Manual refinement of the optimal alignment improved the averaged Q-score by0 .75% (56.5533). The ClustalW default parameters produced a suboptimal result (55.5774) within the GOP-GEP space explored (Figure 1). Expanded phylogeny of Armophorea 577 The alignment based on the corresponding 18S rRNA secondary structure generated a data set for which the averaged Q-score was 49.7616. Although this value appeared to be considerably suboptimal, it should not be interpreted as strictly in SSA since the Q-score is a measure of alignment quality that reflects the minimization of changes among nucleotides, whereas alignments based on secondary structure tend to minimize the changes among RNA structures (Thompson et al., 1994, Hall, 2005Kjer et al., 2007).
The G + C content of the ciliate 18S data was 44.5 mol%. QOA yielded 2,099 characters, of which 1,312 (62.5%) were parsimony-informative, 408 were constant and 379 variable, but parsimony-uninformative. On the other hand, SSA yielded 2,424 characters, of which 1,412 (i.e., 58.3%) were parsimony-informative, 455 were constant, and 557 variable, but parsimony-uninformative. Thus, although QOA provided fewer characters it contained slightly more information for parsimony analysis than SSA. The quantitative difference in the number of characters between the two data sets is explained by SSA having more gap entries than QOA (25.6% vs. 15.3%).
The overall mean p-distance of the 18S data was 0.187 for QOA and 0.171 for SSA; the within-and among-group distances also varied according to the alignment criteria. For the armophorean families examined here, the lowest within-group mean distances were found within the Caenomorphidae in SSA whereas this same family had the largest mean distance in QOA (Table 1). Inter-group distance comparisons of the 18S sequences indicated that the Metopidae and Nyctotheridae were genetically closer to each other than to the Caenomorphidae (Table 2).

Sensitivity to GOP-GEP variation
An analysis of sensitivity to GOP-GEP variation ( Figure 2) indicated that most higher taxa relationships in the Intramacronucleata (especially those of spirotrich clusters) depended upon how indels were estimated during multiple alignments and the influence of such estimates on the calculation of distance matrices in NJ methods. This situation is aggravated by differences in the placement of 'difficult' positions for each combination of parameters. These findings not only corroborate previous observations based on cladistic analyses of implied-aligned matrices of ciliate 18S rRNA (Kivimaki et al., 2009), but also emphasize the need for thorough indel parameter exploration during automated nucleotide alignments (Morrison and Ellis, 1997;Carroll et al., 2006;Smythe et al., 2006;Kjer et al., 2007).

Bayesian inference, Maximum likelihood and Maximum parsimony results
Although a considerable number of sequences from all major ciliate groups was analyzed in this work, the following discussion focuses on the kinships of the Armophorea. BI and ML yielded two topologies (one for each alignment), with the log likelihood of the ML tree ob-tained with SSA being slightly higher than that for QOA (Figures 3 and 4). In MP analyses, the strict consensus of cladograms resulting from SSA provided more resolution than that from QOA. This finding suggested that the former was slightly more robust to variation in Goloboff's K parameter than the latter as it dealt with the relationships among higher taxa and affinities within Metopidae and Nyctotheridae (Figures 5 and 6). The matrices resulting from both QOA and SSA showed considerable character incongruence, as indicated by the relatively low ensemble consistency index of their resulting cladograms. On the other hand, the ensemble retention index was relatively high, indicating that most nucleotide primary homologies contributed to synapomorphy. The cladograms from QOA were rather more consistent and showed slightly more homology than those from SSA. Clades representing the main divergences of armophorean lineages hypothesized from the QOA matrix were united by more synapomorphies than those from SSA, except for the Caenomorphidae (Figures 5 and 6).
In all analyses (BI, ML and MP), the Litostomatea was adelphotaxon of some armophoreans, with high support (> 80%) most of the times. However, a completely monophyletic but relatively weakly supported Armophorea was only hypothesized by the BI and ML trees from QOA (Figure 3), which were significantly different (AU test; p < 0.05) from those in which the armophoreans were not monophyletic. For all other trees, the Armophorea were polyphyletic, and Litostomatea was sister to Metopidae + Nyctotheridae. In these trees, the Caenomorphidae branched outside the Lamellicorticata and diverged at the base of Intramacronucleata or near Licnophora spp. and the remaining spirotrich branch (Table 3). These topologies also differed significantly from those in which the armophoreans were monophyletic (AU test, Templeton test; p < 0.05).
The Caenomorphidae were distributed in a fully supported symmetric group, dichotomized in branches containing four terminals, the position of which varied slightly depending on the alignment criteria and analytic framework. Three Caenomorphid terminals that were classified to family level in NCBI/GenBank, namely AJ009658, AJ009661 and AJ009662, unambiguously branched within the Metopidae (Figures 3-6). Moving these sequences into the main Caenomorphidae branch consistently augmented the mean p-distance of this group but had little effect on that of the Metopidae (Table 1). We therefore suppose that such sequences might belong to actual metopids. They were originally mentioned in a paper by van Hoek et al. (1999) as Caenomorpha-"like" species, so their identity is unknown.
The Metopidae comprised a pectinate line of branches, paraphyletic in relation to the monophyletic Nyctotheridae, and showed the least stable phylogenetic pattern among armophoreans in MP analyses; the latter was seen as polytomies in the consensus trees (Figures 5 and 6), whereas the BI-ML trees of both alignments yielded little inconsistency (Figures 3 and 4). Remarkably, Metopus contortus (Quennerstedt, 1867) Kahl, 1932, always diverged at the base of the Metopidae + Nyctotheridae, and Metopus palaeformis Kahl, 1927, was consistently monophyletic, even though the genus Metopus was not. This situation not only reflects the position of these species in relation to nyctotherids, but also the finding that the branch containing Brachonella arose from within Metopus. The Nyctotheridae were hypothesized to be monophyletic in all analyses, with the affinities of Nyctotherus ovalis Leidy, 1950, sequences varying according to the alignment criteria and phylogenetic framework. The monophyly of Nyctotherus depended on the position of its type species N. velox, which was unstable, although Nyctotheroides was always monophyletic.
Regarding their internal kinships, the monophyly of Armophorea was rejected in all but two analyses (Figures 4-6; Table 3) in which it was weakly supported by the data (Figure 3). Miao et al. (2009b) also found armophoreans to be not monophyletic, but in a different scenario than that hypothesized here. Thus, these authors found Caenomorpha uniserialis to branch off the base of the Litostomatea, while the Metopidae + Nyctotheroidae were a sister group of the protohypotrichs. This contrasts with the study by Shin et al. (2000), in which C. uniserialis branched off the base of a Metopus + Nyctotherus group with moderate to high support (74-98) in distance-based, ML and MP trees.
In his recent classification of the Ciliophora, Lynn (2008) considered the Armophorea to contain two orders (Armophorida and Clevelandellida). The former included the Families Metopidae and Caenomorphidae, while the latter included the Family Nyctotheridae plus five other families that unfortunately were not represented in our study (see Lynn, 2008). The affinity of the Metopidae to the Nyctorheridae contradicts the Order Armophorida proposed by Lynn (2008), who considered metopids to be closely related to caenomorphids. On the other hand, our results seem to fit the system of Jankowski (2007) better, with the caenomorphids, metopids and nyctotherids placed in three separate orders, viz. Armophorida, Metopida Jankowski, 1980, and Clevelandellida, respectively. The sequence of Epalxella antiquorum (Penard, 1922) Corliss, 1960, representative of the Order Odontostomatida Sawaya, 1940, a group traditionally associated with armophoreans (Jankowski, 1964a(Jankowski, ,b, 2007, clustered consistently with Class Plagiopylea Small & Lynn, 1985 (Ventrata) in all analyses (not shown). These findings corroborate a previous study by Stoeck et al. (2007) who found E. antiquorum to be related to trimyemids and plagiopylids. Lynn (2008) thus tentatively transferred the odontostomatids from the Armophorea to the Plagiopylea, but considered that phylogenetic analyses of further representatives and of other markers were necessary in order to decide on their affinity. The placement of Brachonella within the pectinate assemblage of Metopus terminals casts some doubt on the validity of the former, as it involves their morphological separation (see Esteban et al. (1995) and Foissner and Agatha (1999)).
Although the non-monophyly of armophoreans has been reported in the literature (e.g. Miao et al., 2009a,b), it has never been discussed in detail. The unambiguous proximity of the Metopidae and Nyctotheridae is frequent (e.g., Riley and Katz, 2001;Affa'a et al., 2004;Gong et al., 2009), and their position as an adelphotaxon of Litostomatea agrees with the recent study by Vd'acny et al. (2010), who proposed the name Lamellicorticata for the taxon formed by Armophorea and Litostomatea. Accordingly, one putative morphological synapomorphy of this group is the plate-like organization of the postciliary microtubules that form a layer right and between the ciliary rows Expanded phylogeny of Armophorea 579  BI -Bayesian inference, ML -Maximum likelihood, MP -Maximum parsimony, QOA -Q-score optimized alignment; SSA -secondary structure alignment. (Foissner and Agatha, 1999;Lynn, 2008;Vd'acny et al., 2010).
In a detailed fine structure study of Caenomorpha medusula Perty, 1852 (Figures 7 and 8), Santa-Rosa (1975) found that postciliary microtubules were not developed in the somatic kinetids, thus precluding the organization mentioned above (Figures 7C,D). Consequently, if the Armophorea are to be considered monophyletic, then a loss of the plate-like arrangement of postciliary microtubules is assumed to have occurred after the Caenomorphidae lineage diverged at the base of the armophorean cluster (Figure 3). On the other hand, assuming that caenomorphids are distantly related to armophoreans and that no further traditional armophoreans are found to lack the plate-like ar-rangement of postciliary microtubules, the presence of such features can be assumed to be a feasibly consistent synapomorphy of the Lamellicorticata ex Caenomorphidae.
The classification of metopids alongside with caenomorphids is generally based in the assumption of homology of the perizonal ciliary stripe by Small and Lynn (1985) and Puytorac (1994), as discussed by Foissner and Agatha (1999). Foissner and Agatha (1999) described and compared the morphogenetic process in Metopus hasei Sondheim, 1929 andM. inversus (Jankowski, 1964) Foissner and Agatha, 1999 to that described for C. medusula by Martin-Gonzalez et al. (1987) and concluded that they have different morphogenetic origin and function (Foissner and Agatha, 1999). Accordingly, the metopid perizonal stripe 580 Paiva et al.  generates only the paroral for the opisthe, whereas the caenomorphid stripe generates the paroral plus the adoral membranelles and the opisthe's juvenile perizonal stripe (Martin-Gonzalez et al., 1987;Foissner and Agatha, 1999). Additionally, the kinetome organization of metopids differs from that of caenomorphids (Santa-Rosa, 1975;Sola et al., 1990;Silva-Neto, 1993;Decamp and Warren, 1997;Foissner and Agatha, 1999). In M. hasei and M. inversus the somatic dikinetids have a barren anterior kinetosome, while in perizonal dikinetids both kinetosomes were ciliated (Foissner and Agatha, 1999). On the other hand, in C. medusula all of the kinetosomes in the meridian (bell) kineties are ciliated whereas the posterior kinetosome of perizonal dikinetids is barren (Santa-Rosa, 1975) ( Figures  7D and 8C). The three lowermost perizonal dikinetids in M. hasei and M. inversus are not positioned equidistantly, compared to the equidistant placement in C. medusula. While the foregoing features can be used to support hypotheses that the Caenomorphidae are not closely related to the Metopidae + Nyctotheridae, there is morphologic evidence to support the monophyly of Armophorea, although sometimes ambiguously. The most obvious morphological characteristic is the body torsion present in metopids and caenomorphids (Jankowski, 1964b;Corliss, 1979). Among the metopids, this torsion is conspicuous in the campanulate-shaped representatives of Brachonella Jankowski, 1964. Furthermore, Brachonella darwini (Kahl, 1927) Jankowski, 1964b, exhibits a thorn-like posterior projection resembling those of the Caenomorpha (Jankowski, 1964b). The presence of interkinetosomal connectives jointing adoral membranelles of C. medusula (see Figure 8A) characterizes them as heteromembranelles (Puytorac and Grain, 1976;Lynn, 2008) that also occur in clevelandellids (Lynn, 2008).
Such evidence might be considered support for a close relationship between caenomorphids and clevelandellids. However, the presence of heteromembranelles best fits the phylogenetic trees as two independent gains, viz. one in the Caenomorphidae branch and another in the Nyctotheridae (representing the clevelandellids). This also applies even to trees in which Armophorea is monophyletic, given the paraphyly of Metopidae (Figure 3). A diplostichomonad paroral, in which a ridge separates two rows of kinetosomes ( Figure 8B) was found in C. medusula by Santa-Rosa (1975), thus matching this structure's configuration in clevelandellids (Paulin, 1967;Puytorac and Grain, 1976;Takahashi and Imai, 1989;Grim, 1998), but also in the metopid Parametopidium circumlabens (Biggar & Wenrich, 1932) Aescht, 2001(Silva-Neto, 1993. However, this possibly differs from the seemingly linearly arranged oral dikinetids in Metopus (Esteban et al., 1995;Foissner and Agatha, 1999;Lynn, 2008).
The BI-ML trees inferred from the secondary structure alignment and the MP cladograms (Figures 5 and 6; Table 3) also show the possibility of caenomorphids having either diverged at the base of the Intramacronucleata tree, as suggested by Lynn and Wright (2013) or to be related to Spirotrichea. The phylogeny of ciliates based on other markers exhibit different branching patterns (Lynn, 2008) in which metopids and nyctotherids are closely related to spirotrichs. Based on a-tubulin amino acids, Israel et al. (2002) hypothesized a neighbor-joining cluster formed by M. palaeformis + N. ovalis with the spirotrich Euplotes spp. distantly placed from the litostomatean cluster. Moreover, based on histone H4 data, a neighbor-joining tree was hypothesized by Katz et al. (2004) in which M. palaeformis and N. ovalis branched off a trichotomy formed by Spirotrichea and the remaining ciliate clusters (the rooting method was not specified), except for litostomateans, which were not included (Katz et al., 2004). In any case, a-tubulin and H4 phylogenies must be interpreted cautiously because of paralogy (Israel et al., 2002;Katz et al., 2004). These results corroborate the close affinity of metopids to nyctotherids shown by our analyses. However, further data, especially on caenomorphids, is still required to improve our understanding of armophorean kinships based on a-tubulin and H4 phylogenies.

Concluding Remarks
The present study has shown that different nucleotide alignment criteria and the use of different phylogenetic frameworks provided different hypotheses to explain the evolutionary affinities of armophoreans based on the 18S marker. This and the sensitivity of some basal branching patterns of the Intramacronucleata to GOP-GEP variation highlight the importance of explicitness in nucleotide alignment criteria. Whereas the 18S phylogeny of the Armophorea results in an apparently stable placement of metopids and nyctotherids near the litostomateans, the same cannot be said for caenomorphids. Moreover, assumptions regarding the evolution of morphological features based on 18S phylogeny are quite general and ambiguous and must not be over-interpreted since various aspects of the life cycle (which includes morphogenesis) and fine structure of most representatives of Armophorea (Foissner and Agatha, 1999;Lynn, 2008) remain unknown. Improvements in taxon sampling for phylogenetic analyses, the use of additional molecular markers, and advances in our knowledge of the life cycle and fine structure of armophoreans should shed some light on the natural history of these organisms.