Evolutionary Associations of Endosymbiotic Ciliates Shed Light on the Timing of the Marsupial-Placental Split.

Trichostome ciliates are among the most conspicuous protists in the gastrointestinal tract of a large variety of vertebrates. However, little is still known about phylogeny of the trichostome/vertebrate symbiotic systems, evolutionary correlations between trichostome extrinsic traits, and character-dependent diversification of trichostomes. These issues were investigated here, using the relaxed molecular clock technique along with stochastic mapping of character evolution, and binary-state speciation and extinction models. Clock analyses revealed that trichostomes colonized the vertebrate gastrointestinal tract ∼135 Ma, that is, near the paleontological minimum for the split of therian mammals into marsupials and placentals. According to stochastic mapping, the last common ancestor of trichostomes most likely invaded the hindgut of a mammal. Although multiple shifts to fish/amphibian or avian hosts and to the foregut compartments took place during the trichostome phylogeny, only transition to the foregut was recognized as a key innovation responsible for the explosive radiation of ophryoscolecid trichostomes after the Cretaceous/Tertiary boundary, when ungulates began their diversification. Since crown radiations of main trichostome lineages follow those of their mammalian hosts and are in agreement with their historic dispersal routes, the present time-calibrated phylogeny might help to elucidate controversies in the geological and molecular timing of the split between marsupials and placental mammals.

Introduction

The gastrointestinal tract (GIT) of vertebrates is colonized by an extremely diverse, anaerobic microbial community, consisting of archaea, bacteria, fungi, and various protists. Trichostome ciliates are among the most distinctive and largest protists of this complex fermentative association, inhabiting both the foregut and the hindgut compartments of the GIT (Williams and Coleman 1992; Kittelmann et al. 2013; Moon-van der Staay et al. 2014). These ciliates are typically considered as mutualists or harmless commensals, although one species, Neobalantidium coli, can cause severe diarrhea in humans, swines, and rarely also in other domestic and wild animals (Levine 1985; Bradbury 1996; Headley et al. 2008). Mutualistic and commensal trichostomes are very likely not essential for survival of their hosts, as indicated by quite normal life of animals, whose ciliate microbiota has been removed experimentally or has been lost in zoological gardens (Williams and Coleman 1992). In spite of this fact, ciliates obviously contribute to the overall gut function by adding degradative complexity, by their ability to scavenge oxygen, or by their grazing behavior, which helps to shape and regulate prokaryotic populations (Coleman 1986, 1989; Bonhomme 1990). During fermentation of ingested fibre plant material, some trichostome ciliates release large amounts of hydrogen generated in their hydrogenosomes (Yarlett et al. 1984), providing ideal conditions for hydrogenotrophic methanogens.

Trichostomes belong to an ecologically diverse ciliate class, Litostomatea. Their nearest free-living relatives are very likely spathidiids (Strüder-Kypke et al. 2007; Vďačný et al. 2014), that is, predatory ciliates seizing other protists and microscopic animals in a variety of terrestrial or aquatic habitats all over the globe (Foissner and Xu 2007). In spite of the ecological heterogeneity of litostomateans, their monophyletic origin was strongly corroborated by both ultrastructural features and sequence data (Lynn 2008). According to phylogenomic analyses, the Litostomatea are classified within the super-cluster SAL which collates two further ciliate classes, the Spirotrichea and the Armophorea (Gentekaki et al. 2017; Chen et al. 2018).

From an ecological point of view, trichostomes are considered to be rather promiscuous to their hosts, since a number of taxa have been found not only in closely but also comparatively distantly related host species (Williams and Coleman 1992; Dehority 1995; Kittelmann and Janssen 2011). Nonetheless, one distinct association pattern has been recognized, ciliates inhabiting the foregut typically do not occur in the hindgut and vice versa (Moon-van der Staay et al. 2014). Foregut ciliates are typically associated with the forestomach of macropod marsupials (kangaroos and their relatives) (Dehority 1996; Cameron et al. 2001, 2003; Cameron and O’Donoghue 2002a, 2002b, 2002c, 2003a, 2003b, 2003c), with the rumen of even-toed ungulates (Williams and Coleman 1992), or with the crop of birds having a fermentative digestive system (Bardele et al. 2017). On the other hand, hindgut ciliates are associated with fishes, amphibians, and reptiles (Levine 1985; Bradbury 1994), ostriches and rheas (Ponce-Gordo et al. 2002, 2008), primates (Modrý et al. 2009; Tokiwa et al. 2010; Pomajbíková et al. 2012, 2013; Vallo et al. 2012; Schovancová et al. 2013), and especially with elephants and odd-toed ungulates, such as rhinoceroses, horses, and zebras (Ito et al. 2006, 2008, 2010, 2014; Strüder-Kypke et al. 2007; Snelling et al. 2011; Moon-van der Staay et al. 2014).

From an evolutionary standpoint, association patterns between trichostomes, their hosts, and parts of their GIT are puzzling. Rumen ciliates are classified in two independent clades separated by several clusters of hindgut ciliates (Snelling et al. 2011; Ito et al. 2014; Moon-van der Staay et al. 2014); some ciliates isolated from birds occur in the crop, while others in the large intestine and both groups do not cluster together in 18S rRNA gene phylogenies (Bardele et al. 2017); ciliates from the colon of poikilotherm vertebrates also do not form a monophyletic cluster (Grim et al. 2015); and ciliates found in the large intestine of primates, including great African apes and humans, also do not group together (Pomajbíková et al. 2012, 2013; Vallo et al. 2012). To cast more light onto this complex phylogenetic picture and to unravel shifts in extrinsic traits of trichostomes during their evolutionary history, the Bayesian approach implemented in the molecular clock theory along with stochastic mapping of character evolution was employed.

The main goal of this study was to use the trichostome phylogenetic framework as a basis to address some of the outstanding questions about formation of their association patterns with vertebrates in the course of evolution. Specifically, 1) when and where the trichostome/vertebrate symbiotic system might have been established; 2) which part of the vertebrate GIT compartments might have been first colonized by ciliates; 3) what was the key innovation responsible for the explosive radiation of some trichostome lineages in mammals; 4) which extrinsic traits of trichostomes tend to evolve in a correlated fashion; and finally 5) whether the trichostome 18S rRNA gene evolves at the same rate as in their nearest free-living relatives. Addressing these issues might also help to elucidate controversies in the geological and molecular timing of the split between marsupials and placental mammals.

Results and Discussion

Phylogenetic Relationships among Trichostome Ciliates

Results of Bayesian and maximum likelihood (ML) phylogenetic analyses were congruent with previous studies (Strüder-Kypke et al. 2007; Snelling et al. 2011; Ito et al. 2014; Moon-van der Staay et al. 2014; Bardele et al. 2017). The trichostome order Vestibuliferida was revealed as paraphyletic, encompassing the monophyletic order Macropodiniida and the biphyletic order Entodiniomorphida. The order Macropodiniida was fully statistically supported in both phylogenetic analyses and occurrence of its members is restricted to the forestomach of Australian macropod marsupials (Cameron and O’Donoghue 2004). This endemic Australian clade was consistently classified in a sister-group position to the vestibuliferid Balantidium ctenopharyngodoni–B. polyvacuolum cluster with very high statistical support (Bayesian posterior probability 1.00/ML 98% bootstrap value) (supplementary fig. S1, Supplementary Material online).

The order Entodiniomorphida was split into two independent and fully statistically supported lineages, the hindgut family Buetschliidae and a highly diverse cluster comprising all other entodiniomorphid families, whose members inhabit both the foregut and the hindgut compartments of the GIT. Buetschliids grouped inconsistently with very poor statistical support either with the vestibuliferid Balantidium duodeni–B. entozoon cluster in Bayesian analyses, or with the vestibuliferid family Paraisotrichidae in ML analyses. On the other hand, the remaining entodiniomorphid families were depicted as a sister group of the vestibuliferid family Isotrichidae with full support in Bayesian analyses and very strong support in ML analyses (98% bootstrap) (supplementary fig. S1, Supplementary Material online). Relationships among members of the remaining entodioniomorphid families were complex and consistent with those in the study of Ito et al. (2014).

The paraphyletic picture of the order Vestibuliferida indicates that it might represent the stem lineage of the whole trichostome assemblage. Already Strüder-Kypke et al. (2007) argued that vestibuliferids might be the most ancestral trichostomes, because of their plesiomorphic holotrichous ciliation and occurrence in a comparatively large variety of vertebrate higher groups.

Origin and Evolution of the Trichostome-Vertebrate Symbiotic Systems

Since free-living relatives of trichostomes are common in aquatic environments worldwide (Lynn 2008), water was very likely the source from which vertebrates acquired the progenitor of trichostomes, as already suggested by Dogiel (1947). Indeed, experimental studies confirmed that hindgut trichostomes can survive passage through the acidic conditions of the stomach (Adam 1953) and are transmitted via the fecal-oral route and contaminated water (Levine 1985; Bradbury 1996; Egan et al. 2010).

Reconstructions of ancestral traits indicated that the last common progenitor of trichostomes inhabited the hindgut of a mammal (figs. 1 and2). Very likely there were two independent shifts to poikilotherm vertebrates within the vestibuliferid assemblage, one in the Balantidium duodeni–B. entozoon cluster and the other one in the Balantidium ctenopharyngodoni–B. polyvacuolum cluster (supplementary fig. S2, Supplementary Material online). Furthermore, there were three shifts into the foregut (supplementary fig. S3, Supplementary Material online): the macropodiniid clade colonized the forestomach of kangaroos and wallabies, while the isotrichid and ophryoscolecid clades independently invaded the rumen of even-toed ungulates. Colonization of birds and primates also happened independently within the vestibuliferid and entodioniomorphid clades. However, the large intestine of elephants/odd-toed ungulates was invaded only once and very likely by the last common ancestor of vestibuliferids and entodiniomorphids (figs. 1 and2). All these host-site specific associations were also recognized as statistically significant correlates during the trichostome evolution (table 1).

. 1.

SIMMAP reconstruction of ancestral groups of hosts of trichostomes, based on a set of 100 randomly selected post burn-in trees from each run of the Bayesian analysis. Relative proportions of characters states were mapped onto the best scoring ML tree whose branching pattern is shown in detail in supplementary figure S1, Supplementary Material online. Coding of characters is summarized in supplementary table S1, Supplementary Material online.

. 2.

SIMMAP reconstruction of ancestral sites of the gastrointestinal tract of hosts of trichostomes, based on a set of 100 randomly selected post burn-in trees from each run of the Bayesian analysis. Relative proportions of characters states were mapped onto the best scoring ML tree whose branching pattern is shown in detail in supplementary figure S1, Supplementary Material online. Coding of characters is summarized in supplementary table S1, Supplementary Material online.

Table 1.

Evolutionary State-By-State Associations of Trichostome Ciliates with Their Host Taxa and Particular Sites of Their Gastrointestinal Tract.

State-By-State Associations	D Value	P Value	M Value	P Value
Macropod marsupials vs. forestomach	0.0592	0.0431	0.1622	<0.001
Macropod marsupials vs. crop	−0.0005	0.0863	0.0001	0.1255
Macropod marsupials vs. rumen	−0.0219	0.1765	−0.0028	0.4000
Macropod marsupials vs. large intestine	−0.0368	0.0706	−0.0112	0.1647
Fish/amphibians vs. forestomach	−0.0033	0.2078	0.0000	0.7765
Fish/amphibians vs. crop	−0.0003	0.0588	0.0003	0.0533
Fish/amphibians vs. rumen	−0.0183	0.0235	−0.0029	0.0706
Fish/amphibians vs. large intestine	0.0218	0.0275	0.0270	0.0353
Birds vs. forestomach	−0.0010	0.3059	0.0001	0.7569
Birds vs. crop	0.0035	0.0078	0.0142	0.0078
Birds vs. rumen	−0.0013	0.1882	0.0010	0.3569
Birds vs. large intestine	−0.0012	0.2471	0.0001	0.8392
Elephants/odd-toed ungulates vs. forestomach	−0.0299	0.1412	−0.0026	0.5294
Elephants/odd-toed ungulates vs. crop	−0.0023	0.0431	0.0000	0.8627
Elephants/odd-toed ungulates vs. rumen	−0.1206	0.0000	−0.0375	0.0078
Elephants/odd-toed ungulates vs. large intestine	0.1527	0.0000	0.1885	0.0000
Ruminants vs. forestomach	−0.0233	0.2235	−0.0018	0.5020
Ruminants vs. crop	−0.0004	0.0706	0.0010	0.0706
Ruminants vs. rumen	0.1694	<0.001	0.2615	0.0000
Ruminants vs. intestine	−0.1457	<0.001	−0.0748	<0.001
Great apes vs. forestomach	−0.0004	0.4196	0.0001	0.8353
Great apes vs. crop	0.0000	0.0941	0.0001	0.1098
Great apes vs. rumen	−0.0016	0.1961	0.0001	0.8824
Great apes vs. large intestine	0.0021	0.1412	0.0027	0.1608
Pigs vs. forestomach	−0.0014	0.3608	0.0000	0.7529
Pigs vs. crop	−0.0001	0.0627	0.0001	0.0667
Pigs vs. rumen	−0.0056	0.0824	−0.0005	0.1137
Pigs vs. large intestine	0.0071	0.0902	0.0089	0.1098

Note.—D- and M values express the strength of character state correlations, while P values express their statistical significance.

In bold are significant P values (< 0.05 ) for the D Values and M Values

BiSSE analyses (Maddison et al. 2007; FitzJohn 2012) indicated that the foregut trichostome lineages have a significantly higher speciation rate, causing their net-diversification rate to be on an average much higher than in the hindgut lineages (fig. 3). Likelihood ratio test also favored asymmetrical character state transition rates, that is, transition rate from hindgut to foregut is statistically significantly higher (q10 > q01) than vice versa (table 2 and fig. 3). Since colonization of the foregut happened three times independently and no backward shifts were detected by stochastic character mapping (supplementary fig. S3, Supplementary Material online), the association of trichostomes with the foregut may be considered as a key innovation. Apparently, cladogenic success of the foregut trichostomes is coupled with rumination, a peculiar behavior responsible for much better mechanical maceration of the dietary plant material and hence also for new niches available for bacteria and ciliates (Kittelmann et al. 2013; Moon-van der Staay et al. 2014).

. 3.

Posterior densities of differences in speciation, extinction, transition, and net diversification rates between trichostome lineages associated with the foregut (subscript 0) and the hindgut (subscript 1). Differences in rates between foregut and hindgut lineages are shifted toward the positive side of the x-axis, suggesting significant difference in their rates. The average rate of molecular evolution in trichostomes is almost two times higher than in their nearest free-living relatives, spathidiids. Arrows denote median values. Ninety-five percent intervals are indicated below posterior densities.

Table 2.

Parametrization and Fitting of Five Binary State Diversification Models for Trichostome Ciliates Associated with Foregut (subscript 0) and Hindgut (subscript 1) of Their Hosts.

Model	Free Pars	λ0	λ1	μ0	μ1	q01	q10	lnLik	AIC	χ2	P
Unconstrained	6	0.266	0.086	0.234	0.065	0.000	0.004	–600.89	1213.8	–	–
λ0=λ1	5	0.156	0.156	0.114	0.143	0.000	0.001	–606.98	1224.0	12.164	4.87e–04
μ0=μ1	5	0.163	0.129	0.115	0.115	0.000	0.001	–605.05	1220.1	8.3038	3.96e–03
q01=q10	5	0.314	0.073	0.295	0.051	0.002	0.002	–607.05	1224.1	12.306	4.52e–04
λ0=λ1, μ0=μ1, q01=q10	3	0.169	0.169	0.150	0.150	0.001	0.001	–612.61	1231.2	23.434	3.28e–05

Note.—Free pars, free parameters; λ, speciation rate; µ, extinction rate; q, transition rate; lnLik, log likelihood of the model; AIC, Akaike information criterion; χ2, chi–square for likelihood ratio test; P, probability for likelihood ratio test.

Time-Calibrated Evolution of Trichostomes in Light of Origin and Dispersal of Their Hosts

According to molecular clock analyses, the radiation of trichostomes very likely occurred after the beginning of the Cretaceous ∼135 Ma, which is near the paleontological minimum for the split of therian mammals into marsupials and placentals (fig. 4). The known fossil record indicates that this happened in the Early Cretaceous of China ∼125 Ma (Cifelli and Davis 2003; Luo et al. 2003). During the Late Cretaceous, marsupials reached North America, where they further radiated and spread in two directions. One dispersal route was to Australia through South America and Antarctica (Woodburne and Case 1996). The other route was during the Eocene to Eurasia and North Africa, where they became extinct later on. On the other hand, placental mammals spread during the Late Cretaceous from Laurasia to Africa and to South America via North America (Kardong 2015) (fig. 5). After considering all these issues, it could be anticipated that the trichostome/vertebrate symbiotic system might have been established during the Early Cretaceous in Laurasia or North America in the hindgut of an early therian mammal.

. 4.

Chronogram showing posterior estimates of divergence times of trichostome ciliates obtained with the Bayesian relaxed molecular clock technique. The 95% credibility intervals are shown as bars. Horizontal axis represents the time scale in million years.

. 5.

Paleogeography of the Late Cretaceous and dispersal routes of marsupials and placentals (based on Kardong 2015). Marsupials radiated from North America in two directions: 1) to Australia via South America and Antarctica and 2) to Eurasia and North Africa, where they became extinct during the Eocene (dashed arrows). Placentals most likely originated in Laurasia and from there spread to Africa and Americas. Hoatzin, the ancient neotropical bird that is sister to all other landbirds, was very likely infected early after its emergence after the K/T boundary, by vestibuliferid trichostomes originated from marsupials or placentals migrating to South America from North America.

Marsupials most likely brought trichostomes to Australia in their GITs, which is well corroborated by a combination of the following events: 1) marsupials radiated in North America during the Late Cretaceous, when macropodiniid trichostome clade branched off from vestibuliferids (fig. 4); 2) marsupials existed in Australia from the early Paleocene ∼65 Ma (Woodburne and Case 1996) and the marsupial order Diprotodontia, which includes the dominant Australian terrestrial herbivores, radiated in the Paleogene ∼55 Ma (Meredith et al. 2009), similarly as did the endemic Australian macropodiniid trichostome clade (fig. 4); and 3) there were no placentals which could transfer trichostomes to Australia before ∼15 Ma, when Australia moved closer to Indonesia (Rich 1991).

Concerning odd- and even-toed ungulates, both groups originated during the Upper Cretaceous in the Northern Hemisphere. Their crown radiations followed after the Cretaceous-Tertiary (K/T) boundary ∼65 Ma and the Miocene (5–23 Ma) was an important epoch for their further diversification and dispersals between Africa, Eurasia, and North America (Hassanin and Douzery 2003; Springer et al. 2003; Steiner and Ryder 2011; Zhou et al. 2011; Vilstrup et al. 2013). Crown radiations of vestibuliferid and entodiniomorphid lineages inhabiting ungulates and protoungulates basically copy their evolution. Thus, all main vestibuliferid and entodiniomorphid lineages emerged before 65 Ma and the entodiniomorphid family Ophryoscolecidae, which inhabits the rumen of even-toed ungulates, began to intensively radiate slightly after the K/T boundary ∼63 Ma (fig. 4). The explosive burst of new lineages at the base of the ophryoscolecid clade was also recognized by MEDUSA (Alfaro et al. 2009) which indicated that ophryoscolecids radiated without apparent extinction, that is, under a Yule model with a net diversification rate of r = 0.067 lineages per 1 My (lnLik = –577.17, AICc = 1162.49). On the other hand, the rest of the trichostome tree of life diversified circa 3.5 times slower under a birth–death model, with a net diversification rate of r = 0.019 lineages per 1 My and an extinction fraction of ɛ = 0.51 (lnLik = –584.44, AICc = 1172.93).

Invasion of vestibuliferids into fishes and amphibians happened two times independently within the vestibuliferid clade (fig. 1). First, in the early trichostome phylogeny, when the Balantidium duodeni–B. entozoon cluster diverged from buetschliids ∼115 Ma. The second invasion occurred ∼75 Ma, when the Balantidium ctenopharyngodoni–B. polyvacuolum cluster detached from the lineage leading to macropodiniids (fig. 4). Colonization of the bird GIT took place also two times independently. About 75 Ma, when Aviisotricha branched off from the Isotricha–Dasytricha cluster to inhabit the crop of Opisthocomus hoazin. This ancient bird lives in South America and emerged ∼65 Ma (Prum et al. 2015), which is very near the minimum posterior estimate of the Aviisotricha origin. The second invasion into birds happened comparatively recently also within the vestibuliferid cluster, when Neobalantidium coli invaded the large intestine of ostriches and rheas (figs. 1 and4). However, this ciliate is not restricted to birds (Ponce-Gordo et al. 2002, 2008) and its reports from them are only indicative of low host specificity (Levine 1985; Bradbury 1996; Headley et al. 2008) (fig. 1 and supplementary table S1, Supplementary Material online).

The association of trichostomes with primates is well represented by the Troglodytella clade, whose members live only in the large intestine of great African apes (Modrý et al. 2009; Tokiwa et al. 2010; Pomajbíková et al. 2012, 2013; Vallo et al. 2012). Great apes radiated 5.7–7.0 Ma, which almost perfectly matches the present posterior estimate of 5.9 Ma for the Troglodytella diversification (fig. 4). Already Vallo et al. (2012) recognized that T. abrassarti coevolves with chimpanzees and its molecular diversity may cast light also onto problematic parts of their phylogeny.

Trichostomes Elucidate Timing of Marsupial–Placental Split

Molecular data have yielded conflicting results for the timing of the marsupial–placental split and multiple molecular markers suggest divergence times by ∼40–65 Ma older (Kumar and Hedges 1998; Messer et al. 1998; Woodburne et al. 2003) than the available paleontological minimum (Cifelli and Davis 2003; Luo et al. 2003). Interestingly, the present time-calibrated phylogeny of trichostomes indicates that these obligate endosymbionts invaded the GIT of mammals ∼135 Ma (fig. 4), which is comparatively near the fossil appearance of early therians (Luo et al. 2003). Older constraints on the root of the trichostome tree of life failed to produce a joint time prior, documenting that they do not represent proper calibration points (Barba-Montoya et al. 2017). Accordingly, the last common ancestor of trichostomes most likely lived in the hindgut of early therian mammals before they segregated into marsupials and placentals. Without this supposition it is not possible to explain the divergence times and presence of trichostomes in endemic Australian marsupials (fig. 5) and in a variety of placental orders. Because there were no other mammals apart from marsupials that could bring trichostomes to Australia before the Eocene (Rich 1991; Woodburne and Case 1996), marsupials had to be colonized by macropodiniid trichostomes prior to radiating from Laurasia or North America in the Early or Late Cretaceous. The most parsimonious solution seems to be that trichostomes were already present in the most recent common ancestor of marsupials and placentals. In this light, the present time-calibrated phylogeny of trichostomes suggests that mammalian molecular data might have overestimated the timing of the marsupial–placental split. Discrepancies between clock- and fossil-based estimates for divergence dates of mammals might be, in fact, an artifact attributable to relatively small changes in molecular evolutionary rates through time (Beck and Lee 2014).

Symbiosis Might Increase Rate of 18S rRNA Gene Evolution

The present Bayesian dating analysis using MCMCTree (Yang 2007) estimated the clock rate for the trichostome 18S rRNA gene to be 3.06 × 10−4 nucleotide substitutions per site per 1 My on an average, with the 95% credibility interval spanning a range from 2.62 × 10−4 to 3.55 × 10−4. This speed of 18S rRNA gene evolution is at the same order of magnitude as in other groups of ciliates investigated so far, that is, 1.76 × 10−4 on an average (1.13–2.54 × 10−4) in free-living litostomateans (calculated from original data; Rajter and Vďačný 2016); 1.25–3.32 × 10−4 on an average in oligohymenophoreans (Wright and Lynn 1997; Rataj and Vďačný 2018), and 3.96 × 10−4 on an average (1.63–6.25 × 10−4) in armophoreans (Vďačný P, unpublished data).

However, the average rate of molecular evolution in trichostomes is almost two times higher than in their free-living litostomatean relatives (see above and fig. 3). A range of correlates with the increased rate of molecular evolution has been analyzed, including life history traits (e.g., metabolic rate, body size, life spans, generation times, population sizes), environmental factors (temperature, UV radiation, geological complexity), and evolutionary processes (speciation rate, biogeographic history) (Davies et al. 2004; Thomas et al. 2006; Bromham 2011; Yang 2014). In case of trichostomes, the comparatively high and stable body temperature of their predominantly homoiotherm hosts might be responsible for their increased clock rate with respect to free-living litostomateans. Strüder-Kypke et al. (2007) proposed another trigger of faster evolution in trichostomes, host switching. This seems to be a plausible explanation, since armophoreans have the highest substitution rate among ciliates and clevelandellid armophoreans are rather promiscuous to their hosts, occurring in a variety of invertebrates (annelids, mollusks, sea cucumbers, millipedes, and cockroaches) and poikilotherm vertebrates (Earl 1972; Lynn 2008).

Materials and Methods

Alignment and Tree-Building Methods

Operation taxonomic units (OTUs) were selected to represent all main trichostome lineages available in GenBank. Their 18S rRNA gene sequences were aligned on the GUIDANCE2 server (http://guidance.tau.ac.il/ver2/; last accessed January 2018), using the MAFFT algorithm and 100 bootstrap repeats (Sela et al. 2015). Since the score of the resulting alignment was very high (0.9732), no masking strategy was employed. The best fitting evolutionary substitution model of the alignment was estimated and selected under the Akaike Information Criterion in jModelTest ver. 0.1.1 (Guindon and Gascuel 2003; Posada 2008).

Maximum likelihood analyses were performed under the best fitting GTR + I + Γ model in PHYML ver. 3.0 with SPR tree-rearrangement and 1,000 nonparametric bootstrap replicates on the South of France bioinformatics platform (http://www.atgc-montpellier.fr/phyml/; last accessed January 2018) (Guindon et al. 2010). Bayesian analyses were also conducted under the GTR + I + Γ evolutionary model on the CIPRES portal ver. 3.1 (http://www.phylo.org/; last accessed January 2018), using the program MrBayes (Ronquist and Huelsenbeck 2003) on XSEDE ver. 3.2.6 (Miller et al. 2010), with two independent runs each having four chains 5,000,000 generation long. Every 100th tree was sampled and the first 25% of trees were discarded as burn-in. A maximum clade credibility tree was constructed and calculation of posterior probabilities of its nodes was based on the remaining 75% of sampled trees. Maximum likelihood (ML) and Bayesian trees were rooted a posteriori according to the studies of Moon-van der Staay et al. (2014) and Bardele et al. (2017) in FigTree ver. 1.2.3 (Andrew Rambaut, http://tree.bio.ed.ac.uk/software/figtree/; last accessed January 2018).

Calibration Points and Molecular Clock Analyses

Since trichostomes are obligate endosymbionts of the digestive tract and occurrence of some trichostome clades is restricted to certain vertebrate groups, calibration points were based on their fossil appearance. Specifically, 1) a maximum bound of 145 Ma was used on the root, corresponding to fossil remnants of early therian mammals (Cifelli and Davis 2003); 2) a minimum bound of 53 Ma was assigned to the endemic trichostome Australian clade, reflecting the radiation of the marsupial order Diprotodontia (Meredith et al. 2009); 3) a minimum bound of 65 Ma was associated with the node leading to Aviisotricha hoazini which lives exclusively in the crop of Opisthocomus hoazin, representing the most ancient bird lineage (Prum et al. 2015); 4) minimum and maximum bounds of 5.0–7.7 Ma were used for the Troglodytella node, whose members occurs exclusively in great African apes (Pomajbíková et al. 2012, 2013); 5) minimum and maximum bounds of 10–15 Ma were used for the Ostracodium node, whose members live in cattle but not in deer (Williams and Coleman 1992); and 6) a minimum bound of 2 Ma was associated with the Ophryoscolex node, since O. caudatus occurs in sheep and goats, whereas O. purkynjei in cattle (Williams and Coleman 1992). All minimum calibration bounds were assigned the truncated Cauchy distribution with P = 0.1 and c = 0.2, so that the density falls off rapidly away from the mode (Yang and Rannala 2006; Inoue et al. 2010). Whether these calibration points represent a reasonable and joint time prior for Bayesian dating analyses was tested without sequence data in the MCMCTree program (Yang 2007), as recommended by Barba-Montoya et al. (2017).

Before molecular dating analyses, a likelihood ratio test of the molecular clock was performed using BASEML (Yang 2007) on the best scoring ML tree, without fossil calibrations. The log likelihood under the clock was ℓ0 = –16,496.6893, while it was ℓ1 = –16,153.1617 without the clock. Comparison of twice the log likelihood difference 2Δℓ = 2(ℓ1−ℓ0) = 687.06 with a χ2 distribution with 134 df (n−2, where n is the number of OTUs) indicated significant difference with P < 0.001. Consequently, the strict clock could be rejected and divergence times were calculated with a Bayesian relaxed molecular clock technique.

Divergence times were estimated over the best scoring ML tree in MCMCTree (Yang and Rannala 2006; Rannala and Yang 2007; Yang 2007), with the following settings: 1) time unit in million years; 2) the HKY85 + Γ5 evolutionary model; 3) a gamma prior G(6, 2) for the transition/transversion rate ratio κ and a gamma prior G(1, 1) for the gamma shape parameter α; 4) birth and death rates at λ = µ = 2 and sampling fraction at ρ = 0.1, producing a nearly flat kernel density (Yang and Rannala 2006); 5) the mean rate of nucleotide substitution µ at 1.75 × 10−4 per site per 1 My implemented as a diffuse gamma prior G(1, 5,710) (Vďačný 2015); and 6) the variance of the logarithm of the rate σ2 incorporated as a diffuse gamma prior G(1, 145), whose mean represents the reciprocal of prior mean of the root age (see above), as recommended by Thorne et al. (1998). Markov Chain Monte Carlo (MCMC) analyses were run for 40,000 generations, sampling every two iterations, after a burn-in of 2,000 iterations. Two independent MCMC analyses were performed to confirm the consistency of the results obtained. The convergence to the stationary distribution and the effective sample size of divergence time estimates were inspected using the program Tracer ver. 1.6 (Rambaut and Drummond 2007).

Reconstruction of Ancestral Extrinsic Traits

Reconstruction of evolutionary history of four extrinsic traits of trichostomes was conducted in SIMMAP ver. 1.5.2 (Bollback 2006): groups of host organisms, thermoregulation modes of hosts as well as parts and sites of the GIT of hosts. Data were obtained mostly from Williams and Coleman (1992), Cameron and O’Donoghue (2004), Ponce-Gordo et al. (2008), Pomajbíková et al. (2013), Ito et al. (2014), Moon-van der Staay et al. (2014), Bardele et al. (2017) and works cited therein. Characters, their states and coding are summarized in supplementary table S1, Supplementary Material online. All character states were treated as unordered.

Priors for the beta distribution of state frequencies and for the gamma distribution of the overall rate of character change were estimated for each character separately on the best scoring ML tree, with an MCMC analysis implemented in SIMMAP. Samples from the posterior distributions of the MCMC analyses were used to find the best fitting distributions in the R statistical environment, using the R script provided with the SIMMAP package. Priors used for reconstruction analyses are summarized in table 3.

Table 3.

Priors for the Beta Distribution of Character State Frequencies and for the Gamma Distribution of the Overall Rate of Character Change Used in Reconstruction of Ancestral Extrinsic Traits of Trichostome Ciliates.

Extrinsic Trait	Beta Distribution	Gamma Distribution
Hosts of trichostome ciliates	1/k	α=14.02, β=0.66, k=60
Thermoregulation modes of hosts	α=3.00, k=31	α=2.18, β=0.36, k=60
Parts of gastrointestinal tract of hosts	α=4.57, k=31	α=3.75, β=0.66, k=60
Sites of gastrointestinal tract of hosts	1/k	α=6.66, β=0.75, k=60

Note.—For further details on extrinsic traits, see supplementary table S1, Supplementary Material online.

A set of 100 randomly selected post burn-in trees from each run of the Bayesian analysis served to incorporate phylogenetic uncertainty. Ten samples were analyzed per each tree with 20 priors drawn from the prior distributions. Branch lengths were rescaled that the overall length of trees is one. Results were plotted as pie charts and mapped onto the best scoring ML tree using the R script “PlotSimMap.R” (https://github.com/nylander/PlotSimMap; last accessed January 2018).

Correlation of Extrinsic Traits

State-by-state associations of trichostomes to host taxa and sites of their GIT during the evolutionary history was analyzed with the correlation method of Huelsenbeck et al. (2003), as implemented in SIMMAP (Bollback 2006). Like in reconstruction of extrinsic ancestral states, phylogenetic uncertainty was incorporated by a set of 100 randomly selected trees from the posterior distribution of each run of the Bayesian analysis. Ten samples were analyzed per each tree with ten priors drawn from the prior distribution. The number of predictive samples was set to ten, whereby predictive sampling served to determine the posterior P values for correlation testing. To assess the strength of the state-by-state associations, both D- and M values were calculated. The D statistic is a measure of association between one state and another, expressed as frequency of occurrence of states on the phylogeny. On the other hand, the M statistic evaluates the correlation between character histories along the phylogeny for two characters and their states, taking into account the fraction of time one state is associated with another in a character history.

Diversification Analyses of Trichostomes

Changes in tempo of trichostome cladogenesis were assessed by MEDUSA, using a stepwise approach based upon the Akaike information criterion for detecting multiple shifts in birth and death rates (Alfaro et al. 2009). Diversification-rate shifts were identified with the maximum-likelihood technique, applied to the time tree estimated with the Bayesian relaxed molecular clock technique.

The time tree also served as a scaffold for testing the effect of particular parts of the GIT on diversification of trichostome lineages over their evolutionary history. To this end, the binary-state speciation and extinction (BiSSE) approach of Maddison et al. (2007), as implemented in the R-package diversitree (FitzJohn 2012), was employed. The BiSSE model has six free parameters: speciation (λ0, λ1) and extinction (µ0, µ1) rate for both character states of a binary trait as well as a transition rate of change between the two states (q01, q10). These six parameters were variously constrained to specify nested diversification models. Statistically significant differences between constrained models and the full model with six free parameters were assessed with the likelihood ratio test. We minimized the potential of accepting the null hypothesis when the alternate hypothesis is true (type II error) by specifying the proportion of the taxa sampled for each character state, as recommended by Silvestro et al.

Posterior distributions of BiSSE estimates of diversification parameters were also analyzed with MCMC, using an exponential prior. The step argument w for the MCMC sampler was determined by running a 100 step long chain. The range of observed samples was then used as a measure of the “step size” in a 10,000 step long MCMC chain. The 95% credibility intervals for marginal distributions of diversification parameters were calculated and plotted using the R profiles.plot function.

Conclusions

Although trichostomes appear promiscuous to their hosts at first glance, they exhibit a clustering specific for higher taxa of their hosts and individual gastrointestinal compartments. Moreover, crown radiations of main trichostome lineages follow those of their mammalian hosts. This interesting finding might help to elucidate controversies in the geological and molecular timing of the split between marsupials and placental mammals. Whether microbe–host interactions might also explain shifts in the molecular clock rate of the 18S rRNA gene, needs further research. Armophorean and oligohymenophorean ciliates could serve as good model objects for further studies, since they comprise both symbiotic and free-living lineages.

Supplementary Material

Supplementary data are available at Molecular Biology and Evolution online.

Click here for additional data file.