BACKGROUND
progress of molecular techniques has eased the use of genomic data for phylogenetic analyses.
nevertheless, whole genomes are currently available for relatively few metazoans.
molecular studies of phylogenetic relationships within higher taxonomic groups, e.g.
at the intra-ordinal level, therefore still rely on individual genes, among which the nuclear and mitochondrial ribosomal rna genes are the most frequently sequenced.
a pattern of highly variable positions, nested within conserved, slowly substituting sites across the alignment, yields a valuable resource for studying phylogenetic relationships of both recent and ancient splits  <cit> .
this, combined with the ease of amplification, has lead to a widespread use of rrna genes in phylogenetics and furthermore uncovered several specific properties of these genes, which should be considered, using these sequences as phylogenetic markers.
paired regions in rrna sequences evolve via selectively neutral substitutions in the form of compensatory mutations  <cit>  to maintain energetically stable secondary structures.
additionally, a strong bias of nucleotide composition between paired and unpaired areas has been observed  <cit> .
ribosomal rna loop and stem regions are therefore subject to very different selectional regimes, which can hamper the use of rrna genes for phylogenetic purposes  <cit> .
in particular, correlated variation of nucleotides in stem regions, has been suggested to corrupt phylogenetic analyses as these covariation patterns of paired sites do not display independent phylogenetic signal.
ignoring this correlation results in an overestimation of phylogenetic information of these sites, which can lead to inflated measurements of tree robustness  <cit> .
as a solution, rrna secondary structural information as independent set of characters have been advocated to aid tree reconstruction by the use of specific rna substitution models  <cit> .
application of rna substitution models in phylogenetics is still confined to a few studies  <cit> , most of them emphasing improvement of the analyses.
in contrast, a recent study on hexapod phylogeny found mixed rna/dna model setups leading to a higher sensitivity to systematic problems   <cit> .
this has been suggested as a result of potential homoplasy in loop positions.
rna models virtually downweight stem partitions, leading to an increased impact of loops.
if these loop positions are saturated and/or misaligned, this "noisy" signal might dominate the phylogenetic signal of unsaturated stem positions and lead to inaccurate tree reconstruction.
in the present study, we want to test this hypothesis by comparing the performance of mixed rna/dna model setups in the tree reconstruction of different ribosomal rna data sets with the level of relative homoplasy in loop and stem positions.
current studies on the topic of modelling rrna data in tree reconstruction have utilised simulation analyses  <cit> , which can generally be seen to be a sophisticated complement to empirical studies in order to test hypotheses in algorithmically rooted phylogenetics.
however, in letsch et al.
<cit> , tree reconstructions on simulated data were not able to reveal a potential correlation between homoplasy and data modelling.
consequently, the present analyses were based on case studies.
eight ribosomal rna data sets were initially compiled, using from one to three mitochondrial and/or nuclear ribosomal rna gene partitions, covering a broad spectrum of phylogenetic levels , tunicata , heterobranchia  chilopoda , hexapoda , mammalia , primates  and anisoptera ).
all data sets were aligned with the rnasalsa alignment software  <cit> , considering rrna secondary structures.
ambiguously aligned positions where identified and excluded prior to the tree reconstruction.
based on the complete aligned data sets, we further conducted maximum likelihood  tree reconstructions with the raxml v <dig> .
<dig> software package  <cit>  with  a standard dna model setups and   <dig> mixed rna/dna model setups.
in the latter, loop positions are covered by a standard dna model and stem positions are covered by a specific rna model.
performance of different model setups was compared according to recent morphological and molecular expectations of taxonomy.
to test the relative homoplasy between stem and loop regions, all alignments were divided into unpaired  and paired  positions according to a consensus secondary structure.
both partitions were then separately tested for homoplasy by estimating the level of substitutional saturation.
additionally, ml analyses were conducted on loop and stem partitions separately and the results were compared to the trees from the combined data sets.
the analyses setup is depicted in figure  <dig>
RESULTS
phylogenetic analyses
in the following, we represent and discuss the results of echinodermata, tunicata and mammalia data sets as examples.
discussion on the results of all other data sets are provided in additional file  <dig>  to investigate the impact of homoplasy in loop regions on the behaviour of mixed rna/dna models setups in the tree reconstruction, the tree reconstruction results of all  <dig> rna models were compared to the trees relying on the dna model setup.
trees were evaluated in comparisons with recent morphological and molecular understanding of the accordant group, where we mainly focus on "benchmark clades" to check the reliability of each model setup.
clades were defined as "benchmark clades", if they have repeatedly received support in previous studies, based on independent morphological and/or molecular data.
phylogeny of echinoderm classes
echinodermata is divided into five extant classes, the crinoidea , ophiuroidea , asteroidea , holothuroidea  and echinoidea .
monophyly in these five classes is well founded  <cit> , whereas the relationships among them are still debated.
nevertheless, there is some consensus regarding major aspects of echinoderm phylogeny  <cit> .
crinoids are seen as the most basal split within echinodermata, forming the sister group to the four remaining classes .
furthermore, there is strong support for a sister group relationship of echinoids and holothurians .
debates on the phylogenetic position of the stellate forms  revolve around two competing hypotheses: are the ophiurids alone sister group to echinozoa  <cit>  or do asteroids and ophiuroids form a clade , which is then the sister taxon to echinozoa  <cit> ?
the above outlined hypotheses are only reflected by the results of the gtr and the rna6a model setups.
these trees all show basal crinoidea and eleutherozoa divided into asterozoa and echinozoa.
in contrast, all other mixed rna/dna model setups show either holothuroidea or a clade of holothuroidea + echinoidea the sister taxon to the rest of echinodermata , thus clearly contradicting current expectations of echinoderm phylogeny.
the position of appendicularia within tunicata
molecular approaches to the phylogeny of tunicata are generally hampered by the base composition biases and elevated substitution rates in aplousobranchia, appendicularia and mogulidae .
appendicularia retain larval characters throughout their lifespan, which made an understanding of their phylogeny crucial for understanding the evolution of body plans and developmental modes in tunicata  <cit> .
recent molecular studies using phylogenomic or rrna data to target the phylogeny of tunicata usually recover appendicularia as sister group to all other tunicate groups  <cit> .
however, this position is suspected to be a result of a "long branch attraction" artefact, due to genome-wide elevated substitution rates in this group  <cit> .
as an alternative, appendicularia as sister to stolidobranchia has been recovered through analyses of 18s rrna genes  <cit> .
however, this position was generally weakly supported and has been discussed as a possible result of base composition bias in appendicularia and mogulidae, a family of stolidobranchia  <cit> .
these problems are reflected by the results of our study on tunicates , which either show appendicularia as first split within tunicata  or as sister group to stolidobranchia .
according to the currently unresolved position of appendicularia, none of these alternatives can be chosen as superior.
relationships within mammalia
the first molecular analyses on the phylogeny of mammalia, using mitochondrial genes have remarkably challenged previous morphological hypotheses on the relationship among mammalian groups  <cit> .
however, subsequent studies on nuclear markers and more sophisticated analyses of mitochondrial genomes led to more consistent hypotheses of mammalian relationships, which are in several aspects congruent to morphological studies  <cit> .
general congruence among these independent markers has resulted in a well resolved and strongly corroborated backbone tree of mammalian groups, representing the four "superorders" xenarthra, afrotheria, euarchontoglires and laurasiatheria and several subgroups, e.g.
monophyletic theria, the paenungulata , tetytheria , and euarchonta .
consequently, the evaluation of tree reconstructions targeting mammalian phylogeny has been formalised by defining theses groups as "benchmark clades"  <cit> , whose appearance is used to evaluate the performance of the method that was used.
in the present study, most analyses are highly congruent in their results of mammalian relationships and display many of the proposed benchmark clades with sufficient support.
however, in the rna6a, rna <dig> and rna16a analyses, a clade combining afrotheria + xenarthra is sistergroup to laurasiatheria and rodentia appear paraphyletic to the remaining eutherian groups, with muroidea + anomalurus as first split within eutheria.
this reflects the suggestions of several previous analyses based on mt genes, but must been interpreted as a result of model misspecification ignoring among-site rate variation  <cit>  and compositional bias  <cit> .
in contrast, all other analyses show afrotheria + xenarthra as first split within eutheria and paraphyletic rodentia, but the latter are nested within monophyletic laurasiatheria.
it is notable, that bootstrap support values for potentially correct groupings increased, if mixed rna/dna model setups are applied .
separate analyses of loops and stems
homoplasy due to multiple substitutions was tested with the index of substitution saturation   <cit> , which assumes a critical index of substitution saturation  that defines a threshold for significant saturation in the data.
the issc is compared with the observed iss of the data.
if the iss value is larger than the critical issc values, saturation is assumed.
to contribute to different tree shapes, the issc is estimated, using a symmetrical  and an asymmetrical  tree topology.
the test for homoplasy reveals striking differences between paired and unpaired positions.
in the stem portions of all data sets, the asymmetrical and the symmetrical issc is always larger than the observed iss.
the differences are significant, thus indicating that the paired partitions are not saturated.
in contrast, we detected potential saturation of substitution in the unpaired positions of several data sets.
for echinodermata, hexapoda, tunicata, chilopoda and heterobranchia, the iss of was notably larger than the asymmetrical issc, suggesting substantial saturation in these alignments.
the complete results of the saturation tests are summarised in table  <dig>  the comparison matrix in figure  <dig> further depicts the results of all tree reconstruction analyses in relation of substitution saturation in the loop partitions.
the saturation test results of the loop and stem partitions were additionally compared to the saturation test results of the combined data sets of the groups exhibiting saturation in the loop regions.
as displayed in table  <dig> saturation vanishes in all of the combined data sets.
characteristics of the applied test data sets including the results of the test for substitution saturation in loop and stem partitions.
iss: estimated index of substitution saturation for the data set.
iss.c s and iss.c a: critical values for the index of substitution saturation  if the true tree is symmetrical  or asymmetrical .
iss > iss.c indicates saturation.
*original
results of the index of substitution saturation test for combined data sets, where homoplasy due to multiple substitution had been previously been detected in the loop partitions.
iss.c s and iss.c a: critical values for the index of substitution saturation  if the true tree is symmetrical  or asymmetrical .
iss > iss.c indicates saturation.
subsequently, ml tree reconstructions were conducted on the separated loop and stem partitions, using a dna model setup.
to characterise the phylogenetic signal in both partitions, we checked whether the trees from the paired or the unpaired partition were more congruent to the combined data results.
trees resulting from all three setups  were compared with the robinson-foulds  <cit>   distance score, which accounts for topology differences.
this indicates a closer similarity of trees based on combined and unpaired data.
comparisons between the combined data set, analysed under different mixed model schemes, usually strengthen this effect.
with the exception of chilopda and hexapoda, most rf distances between the combined data and the unpaired data diminished, whereas rf distances between the combined data and the paired data often increased .
discussion
relative homoplasy and rna modelling
to our knowledge, this is the first work on a separate characterization of homoplasy in paired and unpaired regions of rrna sequences.
we examined relative homoplasy separately, as due to their distinct physiological function in protein biosynthesis, paired and unpaired positions can be expected to evolve differentially and might therefore also differ in their sensitivity to numerous substitutions occurring at the same position, thus hiding or completely erasing phylogenetic signal.
our results indicate separate exploratory analyses of loops and stems as crucial, because homoplasy due to multiple substitutions in loop positions could not be detected if the combined data sets  are tested for excessive homoplasy .
pooling loop portions may be subject to the same kind of underestimation of homoplasy, if rates among loop regions are highly heterogeneous, as is apparent in van de peer et al.
<cit> .
thus, it may be advisable to estimate homoplasy in each loop region separately.
these homoplastic substitution patterns have generally been addressed as "substitutional saturation"  <cit> .
however, "saturation" is a concept that is only relevant to distance based analyses, where "saturation" refers to saturation curves, in which increasing phylogenetic depth does not increase pairwise distances.
phylogenies have nodes at many levels, from the tips to the root, and character based analyses can insulate homoplasy and mediate errors due to homoplastic sites in ways that distance based analyses cannot  <cit> , for example by increasing the taxon sampling to break up long branches  <cit> .
in this context, it can further be stated that "saturation" is not an inherent character of an aligned sequence position in a given alignment, it rather depends on the considered phylogenetic level.
additional measurements of the iss in the tunicate data set, applied to the loop partitions of distinct monophyletic subgroups within tunicata, shows an iss significantly smaller than the issc for both symmetrical and asymmetrical topologies, thus indicating a decrease of homoplasy in these groups .
the iss method applied here accounts for position specific nucleotide frequency pattern in a given alignment, which are supposed to reflect the occurrence of multiple substitutions in the data set  <cit> , thus implying relative homoplasy.
therefore, the use of the iss might be a reasonable heuristic to estimate the level of homoplasy according to the deep splits discussed in the considered data sets.
results of the index of substitution saturation test for different tunicate subgroups.
iss.c s and iss.c a: critical values for the index of substitution saturation  if the true tree is symmetrical  or asymmetrical .
iss > iss.c indicates saturation.
as depicted in figure  <dig>  most of the rna models only lead to reasonable tree hypotheses, if the loop regions are found to contain meaningful phylogenetic data.
in data sets identified as saturated, nearly all rna models failed to recover an expected hypothesis.
in contrast, the standard dna model setup usually led to trees congruent with recent views on the relationships in the accordant groups.
in this context, especially the results of the hexapod data are interesting, as these analyses did not provide superior trees by the dna models setups, although saturation is detected.
in the case of hexapods, no model setup led to the expected tree hypothesis.
this indicates a generally insufficient phylogenetic signal of the 18s and 28s rna data to resolve the shortest of the internodes among the deep hexapod splits.
however, in  <cit> , bayesian inference analyses of the identical data set led to wrong tree hypotheses in the mixed rna/dna model setup, but not in the standard dna model setup.
kjer  <cit>  found that topologies were virtually identical between the standard and rna models, but support values varied: in some cases, expected nodes were more strongly supported with gtr models  other cases favoured the doublet model .
similarly, support for a paraphyletic anisoptera  went down with the doublet model.
nevertheless it is noteworthy, that mixed rna/dna model setups frequently led to an increased support for probably correct clades in data sets without homoplastic loops , thus supporting previous studies on rna models in phylogenetics  <cit> .
our analyses setup also allows comparisons between the different rna model setups.
current rna substitution models can be divided into two distinct classes, rooted in population genetics  <cit> .
models of the first class assume a one-step process of compensatory substitution in paired positions, thus allowing double substitutions : a mutation in a base pair  may led to slightly deleterious ug or gu pairs.
if selection against these intermediates is strong, these are kept in low frequency in the population.
if a second mutation occurs at the corresponding site , drift in gene frequency may lead to a domination of this new base pairing in the population .
in contrast, models of the second class assume a two-step process of compensatory substitution in paired positions, considering one substitution in base-pairs, with a probability of zero for all double substitutions.
this approach considers a fixation of intermediate states in the population at a high frequency, after a mutation in a base pair and before a second mutation at the corresponding site .
rna models can be further discriminated by their treatment of mismatches: 6-state models completely ignore these pairings, 7-state models lump all mismatches in one category, whereas 16-state models apply distinct frequency and substitution rate parameters to the individual mismatches.
previous studies on rna models in phylogenetics have predicted the superiority of models allowing double substitutions and the superiority of the most general models   <cit> .
the latter is corroborated by the aicc modeltest of the present analyses, which frequently show higher likelihoods and aicc values for the most general models .
furthermore, the rna6a model led to the expected topologies in two data sets  showing relative homoplasy in loop partitions, whereas all other rna models fail to display presumably correct trees, if significant homoplasy was identified .
in this context, the rna6d and rna16b models are performing worst, as both are only able to display one potentially correct tree hypotheses.
additionally, congruencies between the performance of rna models and the results of the aicc modeltest can be drawn from our results.
according to the aicc modeltest, in the class of the rna <dig> models, the most general rna6a model is always superior to all other rna <dig> models , which is further congruent to its performance in tree reconstruction analyses.
this is not reflected by the rna <dig> and rna <dig> models, where the models with the highest aicc scores .
potential pitfalls of rna modelling
consequences of different evolutionary constraints in stem and loop regions of rrna sequences for phylogenetic analyses has long been suspected and led to different recommendations for weighting stem positions in parsimony analyses  <cit> .
beside suggestions for simple one-half weighting of paired positions  <cit> , empirical investigation of compensatory substitution rates in stem positions  <cit>  reveals a rate of about 40% of that expected under a hypothesis of perfect compensation.
therefore, the weighting of stem characters is suggested to be reduced by no more than 20%.
consequently, in model based tree reconstruction methods, like bayesian inference and maximum likelihood, it should be reasonable to use specific rna models  as simply applying a standard dna model to data from one part of the helical regions.
application of these rna models has frequently been justified by a consistent phylogenetic signal of coevolved paired sites, decreasing the information content in the data  <cit> .
analyses ignoring this interdependence should tend to overestimate the support for dubious or even wrong nodes in a tree  <cit> .
due to a reduced number of effective sites, the application of specific rna models, which take interdependencies into account, reduces tree confidence, but is more reliable in the light of the information content in the data  <cit> .
our results actually imply a reduced impact of stem positions in the combined data set, if mixed rna/dna model setup are used.
this is depicted by the tree distance results of the separate analyses of the stem and loop partitions.
in most data sets, the distances between the trees based on only the loop partition and the combined data are reduced, if rna models are applied for the combined data, whereas the distances between the stem partition and the combined data are mostly enlarged .
this could have been expected, if coevolution in paired sites is assumed and thus these positions do not provide independent phylogenetic information.
however, for several of the currently tested data sets, the substitution saturation test reveals that the unpaired positions clearly experience excessive homoplasy, which indicates a loss of phylogenetic information, as these positions are no longer informative  <cit> .
in this context, stem positions in the current data sets should contain more reliable signal, compared to loop regions, as they exhibit a much lesser grade of homoplasy due to multiple substitutions.
consequently, the application of rna models increases the relative impact of noisy positions in the data set and reduces the influence of more informative portions.
thus, the results of the current analyses corroborate the hypothesis proposed by letsch et al.
<cit> .
rna models doubtlessly provide a better depiction of the phylogenetic information content of rrna data sets, but this might be a trap, if homoplasy is far greater in loop positions.
in this case, the informative phylogenetic content is obscured by noise.
the situation might probably be depicted best, if we thought of a weighting scheme for rrna data sets: in standard dna model setup, the signal of paired positions is virtually weighted twice, as both positions are linked and signal of pairs can be seen redundant.
as outlined above, previous studies have mostly targeted this as a problem  <cit> , but the current analyses showed relative homoplasy as delimiting the confidence of the phylogenetic signal provided by loop regions, revealing a more or less hidden coherence between two factors - covariation and homoplasy - contributing to the phylogenetic signal of the rrna data sets.
for this scenario it can be stated, that in contrast to a previously proposed overestimation of wrong support by ignoring site interdependencies  <cit> , the application of rna models will tend to overestimate the support for dubious or wrong nodes in a tree.
as depicted above, loops and stems can be expected to experience different selectional regimes, which has resulted in the development of the specific rna models.
nevertheless, it has been noted as early as  <dig>  <cit>  that substitution rates do not fit neatly into stem-loop partitions, and thus weighting according to stems vs. loops might be problematic, which was later demonstrated by van de peer  <cit> .
consequently, selectional constraints on rrna may not only differ between paired and unpaired regions, but also among the individual loop or stem regions, which would depend on the individual function and their relative location in the 3d rrna molecule.
binding sites of ribosomal proteins, for example protein l11-binding domain  within the lsu rrna domain ii and the sarcin-ricin loop within domain vi, constituting the gtpase-associated center  <cit>  or the lsu rrna domain v, which contains the peptidyl transferase center   <cit> , are highly conserved throughout metazoa.
furthermore, many of the rrna regions of domain iv that are involved in trna and inter-subunit interactions are also preserved  <cit> .
in contrast, the domain i of the mt lsu is highly variable on sequence level and until now, no conserved secondary structures could be detected  <cit> .
consequently, the partitioning into loops and stems must be seen as only an relative coarse approximation to model rrna sequences.
in future phylogenetic studies on rrna, more sophisticated partitioning schemes, depending on the function, base composition and relative location of rrna regions, would be able to enhance model based tree reconstruction analyses
CONCLUSIONS
the results of the present study can be interpreted as a trade-off between using specific rna models for a hopefully more accurate estimation of covariation in paired sites and the risk of augmenting relatively homoplastic unpaired positions in the tree reconstruction.
for future phylogenetic studies based on rrna sequences, we would therefore highly recommend a separate test for saturation of substitution in loop and stem partitions of the aligned data set.
the use of a mixed rna/dna model setup should be avoided if saturation occurs in the loop partitions, as otherwise the valuable phylogenetic signal of the stem partitions might be masked by potentially noisy signal provided by the loops.
in contrast, if no substantial homoplasy is detected in the data, the use of mixed rna/dna models can be highly recommend, as these lead to an increased support to probably correct clades.
based on the presents results, we cannot advocate an general exclusion of the potentially noisy loop positions: first, noise is not an inherent character of a certain nucleotide position, but depends on the considered phylogenetic level.
and second, differences among loop  regions can be expected and excluding these regions as a whole reduces the phylogenetic signal of the data set.
consequently, we rather recommend to think about enhanced partitioning strategies, which would allow a more careful modelling of rrna sequences and provide a first approach to detect noisy signal among loop  partitions.
however, covariation and substitution saturation are only two parameters of the evolutionary inherent pattern displayed  in the data.
other phenomena, like non-stationarity of substitution rates across sites and branches as well as base composition heterogeneity, might also maul the signal content of the data set.
a previous study  <cit>  based on nuclear rrna genes, identified deviation of base composition in certain clades as probably misleading tree reconstruction analyses, rather than the covariation pattern in stem regions.
a sophisticated estimation of evolutionary pattern in rrna sequence data is therefore principally desirable and newly developed methods should be applied, which are able to consider background knowledge as covariation, non-stationary processes or heterogeneity in the data  <cit> .

 progress of molecular techniques has eased the use of genomic data for phylogenetic analyses. nevertheless, whole genomes are currently available for relatively few metazoans. molecular studies of phylogenetic relationships within higher taxonomic groups, e.g. at the intra-ordinal level, therefore still rely on individual genes, among which the nuclear and mitochondrial ribosomal rna genes are the most frequently sequenced. a pattern of highly variable positions, nested within conserved, slowly substituting sites across the alignment, yields a valuable resource for studying phylogenetic relationships of both recent and ancient splits. this, combined with the ease of amplification, has lead to a widespread use of rrna genes in phylogenetics and furthermore uncovered several specific properties of these genes, which should be considered, using these sequences as phylogenetic markers. paired regions in rrna sequences evolve via selectively neutral substitutions in the form of compensatory mutations to maintain energetically stable secondary structures. additionally, a strong bias of nucleotide composition between paired and unpaired areas has been observed. ribosomal rna loop and stem regions are therefore subject to very different selectional regimes, which can hamper the use of rrna genes for phylogenetic purposes. in particular, correlated variation of nucleotides in stem regions, has been suggested to corrupt phylogenetic analyses as these covariation patterns of paired sites do not display independent phylogenetic signal. ignoring this correlation results in an overestimation of phylogenetic information of these sites, which can lead to inflated measurements of tree robustness. as a solution, rrna secondary structural information as independent set of characters have been advocated to aid tree reconstruction by the use of specific rna substitution models. application of rna substitution models in phylogenetics is still confined to a few studies, most of them emphasing improvement of the analyses. in contrast, a recent study on hexapod phylogeny found mixed rna/dna model setups leading to a higher sensitivity to systematic problems. this has been suggested as a result of potential homoplasy in loop positions. rna models virtually downweight stem partitions, leading to an increased impact of loops. if these loop positions are saturated and/or misaligned, this "noisy" signal might dominate the phylogenetic signal of unsaturated stem positions and lead to inaccurate tree reconstruction. in the present study, we want to test this hypothesis by comparing the performance of mixed rna/dna model setups in the tree reconstruction of different ribosomal rna data sets with the level of relative homoplasy in loop and stem positions. current studies on the topic of modelling rrna data in tree reconstruction have utilised simulation analyses, which can generally be seen to be a sophisticated complement to empirical studies in order to test hypotheses in algorithmically rooted phylogenetics. however, in letsch et al., tree reconstructions on simulated data were not able to reveal a potential correlation between homoplasy and data modelling. consequently, the present analyses were based on case studies. eight ribosomal rna data sets were initially compiled, using from one to three mitochondrial and/or nuclear ribosomal rna gene partitions, covering a broad spectrum of phylogenetic levels, tunicata, heterobranchia chilopoda, hexapoda, mammalia, primates and anisoptera ). all data sets were aligned with the rnasalsa alignment software, considering rrna secondary structures. ambiguously aligned positions where identified and excluded prior to the tree reconstruction. based on the complete aligned data sets, we further conducted maximum likelihood tree reconstructions with the raxml v. software package with a standard dna model setups and mixed rna/dna model setups. in the latter, loop positions are covered by a standard dna model and stem positions are covered by a specific rna model. performance of different model setups was compared according to recent morphological and molecular expectations of taxonomy. to test the relative homoplasy between stem and loop regions, all alignments were divided into unpaired and paired positions according to a consensus secondary structure. both partitions were then separately tested for homoplasy by estimating the level of substitutional saturation. additionally, ml analyses were conducted on loop and stem partitions separately and the results were compared to the trees from the combined data sets. the analyses setup is depicted in figure  phylogenetic analyses in the following, we represent and discuss the results of echinodermata, tunicata and mammalia data sets as examples. discussion on the results of all other data sets are provided in additional file to investigate the impact of homoplasy in loop regions on the behaviour of mixed rna/dna models setups in the tree reconstruction, the tree reconstruction results of all rna models were compared to the trees relying on the dna model setup. trees were evaluated in comparisons with recent morphological and molecular understanding of the accordant group, where we mainly focus on "benchmark clades" to check the reliability of each model setup. clades were defined as "benchmark clades", if they have repeatedly received support in previous studies, based on independent morphological and/or molecular data. phylogeny of echinoderm classes echinodermata is divided into five extant classes, the crinoidea, ophiuroidea, asteroidea, holothuroidea and echinoidea. monophyly in these five classes is well founded, whereas the relationships among them are still debated. nevertheless, there is some consensus regarding major aspects of echinoderm phylogeny. crinoids are seen as the most basal split within echinodermata, forming the sister group to the four remaining classes. furthermore, there is strong support for a sister group relationship of echinoids and holothurians. debates on the phylogenetic position of the stellate forms revolve around two competing hypotheses: are the ophiurids alone sister group to echinozoa or do asteroids and ophiuroids form a clade, which is then the sister taxon to echinozoa? the above outlined hypotheses are only reflected by the results of the gtr and the rna6a model setups. these trees all show basal crinoidea and eleutherozoa divided into asterozoa and echinozoa. in contrast, all other mixed rna/dna model setups show either holothuroidea or a clade of holothuroidea + echinoidea the sister taxon to the rest of echinodermata, thus clearly contradicting current expectations of echinoderm phylogeny. the position of appendicularia within tunicata molecular approaches to the phylogeny of tunicata are generally hampered by the base composition biases and elevated substitution rates in aplousobranchia, appendicularia and mogulidae. appendicularia retain larval characters throughout their lifespan, which made an understanding of their phylogeny crucial for understanding the evolution of body plans and developmental modes in tunicata. recent molecular studies using phylogenomic or rrna data to target the phylogeny of tunicata usually recover appendicularia as sister group to all other tunicate groups. however, this position is suspected to be a result of a "long branch attraction" artefact, due to genome-wide elevated substitution rates in this group. as an alternative, appendicularia as sister to stolidobranchia has been recovered through analyses of 18s rrna genes. however, this position was generally weakly supported and has been discussed as a possible result of base composition bias in appendicularia and mogulidae, a family of stolidobranchia. these problems are reflected by the results of our study on tunicates, which either show appendicularia as first split within tunicata or as sister group to stolidobranchia. according to the currently unresolved position of appendicularia, none of these alternatives can be chosen as superior. relationships within mammalia the first molecular analyses on the phylogeny of mammalia, using mitochondrial genes have remarkably challenged previous morphological hypotheses on the relationship among mammalian groups. however, subsequent studies on nuclear markers and more sophisticated analyses of mitochondrial genomes led to more consistent hypotheses of mammalian relationships, which are in several aspects congruent to morphological studies. general congruence among these independent markers has resulted in a well resolved and strongly corroborated backbone tree of mammalian groups, representing the four "superorders" xenarthra, afrotheria, euarchontoglires and laurasiatheria and several subgroups, e.g. monophyletic theria, the paenungulata, tetytheria, and euarchonta. consequently, the evaluation of tree reconstructions targeting mammalian phylogeny has been formalised by defining theses groups as "benchmark clades", whose appearance is used to evaluate the performance of the method that was used. in the present study, most analyses are highly congruent in their results of mammalian relationships and display many of the proposed benchmark clades with sufficient support. however, in the rna6a, rna and rna16a analyses, a clade combining afrotheria + xenarthra is sistergroup to laurasiatheria and rodentia appear paraphyletic to the remaining eutherian groups, with muroidea + anomalurus as first split within eutheria. this reflects the suggestions of several previous analyses based on mt genes, but must been interpreted as a result of model misspecification ignoring among-site rate variation and compositional bias. in contrast, all other analyses show afrotheria + xenarthra as first split within eutheria and paraphyletic rodentia, but the latter are nested within monophyletic laurasiatheria. it is notable, that bootstrap support values for potentially correct groupings increased, if mixed rna/dna model setups are applied. separate analyses of loops and stems homoplasy due to multiple substitutions was tested with the index of substitution saturation, which assumes a critical index of substitution saturation that defines a threshold for significant saturation in the data. the issc is compared with the observed iss of the data. if the iss value is larger than the critical issc values, saturation is assumed. to contribute to different tree shapes, the issc is estimated, using a symmetrical and an asymmetrical tree topology. the test for homoplasy reveals striking differences between paired and unpaired positions. in the stem portions of all data sets, the asymmetrical and the symmetrical issc is always larger than the observed iss. the differences are significant, thus indicating that the paired partitions are not saturated. in contrast, we detected potential saturation of substitution in the unpaired positions of several data sets. for echinodermata, hexapoda, tunicata, chilopoda and heterobranchia, the iss of was notably larger than the asymmetrical issc, suggesting substantial saturation in these alignments. the complete results of the saturation tests are summarised in table the comparison matrix in figure further depicts the results of all tree reconstruction analyses in relation of substitution saturation in the loop partitions. the saturation test results of the loop and stem partitions were additionally compared to the saturation test results of the combined data sets of the groups exhibiting saturation in the loop regions. as displayed in table saturation vanishes in all of the combined data sets. characteristics of the applied test data sets including the results of the test for substitution saturation in loop and stem partitions. iss: estimated index of substitution saturation for the data set. iss.c s and iss.c a: critical values for the index of substitution saturation if the true tree is symmetrical or asymmetrical. iss > iss.c indicates saturation. *original results of the index of substitution saturation test for combined data sets, where homoplasy due to multiple substitution had been previously been detected in the loop partitions. iss.c s and iss.c a: critical values for the index of substitution saturation if the true tree is symmetrical or asymmetrical. iss > iss.c indicates saturation. subsequently, ml tree reconstructions were conducted on the separated loop and stem partitions, using a dna model setup. to characterise the phylogenetic signal in both partitions, we checked whether the trees from the paired or the unpaired partition were more congruent to the combined data results. trees resulting from all three setups were compared with the robinson-foulds distance score, which accounts for topology differences. this indicates a closer similarity of trees based on combined and unpaired data. comparisons between the combined data set, analysed under different mixed model schemes, usually strengthen this effect. with the exception of chilopda and hexapoda, most rf distances between the combined data and the unpaired data diminished, whereas rf distances between the combined data and the paired data often increased. discussion relative homoplasy and rna modelling to our knowledge, this is the first work on a separate characterization of homoplasy in paired and unpaired regions of rrna sequences. we examined relative homoplasy separately, as due to their distinct physiological function in protein biosynthesis, paired and unpaired positions can be expected to evolve differentially and might therefore also differ in their sensitivity to numerous substitutions occurring at the same position, thus hiding or completely erasing phylogenetic signal. our results indicate separate exploratory analyses of loops and stems as crucial, because homoplasy due to multiple substitutions in loop positions could not be detected if the combined data sets are tested for excessive homoplasy. pooling loop portions may be subject to the same kind of underestimation of homoplasy, if rates among loop regions are highly heterogeneous, as is apparent in van de peer et al.. thus, it may be advisable to estimate homoplasy in each loop region separately. these homoplastic substitution patterns have generally been addressed as "substitutional saturation". however, "saturation" is a concept that is only relevant to distance based analyses, where "saturation" refers to saturation curves, in which increasing phylogenetic depth does not increase pairwise distances. phylogenies have nodes at many levels, from the tips to the root, and character based analyses can insulate homoplasy and mediate errors due to homoplastic sites in ways that distance based analyses cannot, for example by increasing the taxon sampling to break up long branches. in this context, it can further be stated that "saturation" is not an inherent character of an aligned sequence position in a given alignment, it rather depends on the considered phylogenetic level. additional measurements of the iss in the tunicate data set, applied to the loop partitions of distinct monophyletic subgroups within tunicata, shows an iss significantly smaller than the issc for both symmetrical and asymmetrical topologies, thus indicating a decrease of homoplasy in these groups. the iss method applied here accounts for position specific nucleotide frequency pattern in a given alignment, which are supposed to reflect the occurrence of multiple substitutions in the data set, thus implying relative homoplasy. therefore, the use of the iss might be a reasonable heuristic to estimate the level of homoplasy according to the deep splits discussed in the considered data sets. results of the index of substitution saturation test for different tunicate subgroups. iss.c s and iss.c a: critical values for the index of substitution saturation if the true tree is symmetrical or asymmetrical. iss > iss.c indicates saturation. as depicted in figure most of the rna models only lead to reasonable tree hypotheses, if the loop regions are found to contain meaningful phylogenetic data. in data sets identified as saturated, nearly all rna models failed to recover an expected hypothesis. in contrast, the standard dna model setup usually led to trees congruent with recent views on the relationships in the accordant groups. in this context, especially the results of the hexapod data are interesting, as these analyses did not provide superior trees by the dna models setups, although saturation is detected. in the case of hexapods, no model setup led to the expected tree hypothesis. this indicates a generally insufficient phylogenetic signal of the 18s and 28s rna data to resolve the shortest of the internodes among the deep hexapod splits. however, in, bayesian inference analyses of the identical data set led to wrong tree hypotheses in the mixed rna/dna model setup, but not in the standard dna model setup. kjer found that topologies were virtually identical between the standard and rna models, but support values varied: in some cases, expected nodes were more strongly supported with gtr models other cases favoured the doublet model. similarly, support for a paraphyletic anisoptera went down with the doublet model. nevertheless it is noteworthy, that mixed rna/dna model setups frequently led to an increased support for probably correct clades in data sets without homoplastic loops, thus supporting previous studies on rna models in phylogenetics. our analyses setup also allows comparisons between the different rna model setups. current rna substitution models can be divided into two distinct classes, rooted in population genetics. models of the first class assume a one-step process of compensatory substitution in paired positions, thus allowing double substitutions: a mutation in a base pair may led to slightly deleterious ug or gu pairs. if selection against these intermediates is strong, these are kept in low frequency in the population. if a second mutation occurs at the corresponding site, drift in gene frequency may lead to a domination of this new base pairing in the population. in contrast, models of the second class assume a two-step process of compensatory substitution in paired positions, considering one substitution in base-pairs, with a probability of zero for all double substitutions. this approach considers a fixation of intermediate states in the population at a high frequency, after a mutation in a base pair and before a second mutation at the corresponding site. rna models can be further discriminated by their treatment of mismatches: 6-state models completely ignore these pairings, 7-state models lump all mismatches in one category, whereas 16-state models apply distinct frequency and substitution rate parameters to the individual mismatches. previous studies on rna models in phylogenetics have predicted the superiority of models allowing double substitutions and the superiority of the most general models. the latter is corroborated by the aicc modeltest of the present analyses, which frequently show higher likelihoods and aicc values for the most general models. furthermore, the rna6a model led to the expected topologies in two data sets showing relative homoplasy in loop partitions, whereas all other rna models fail to display presumably correct trees, if significant homoplasy was identified. in this context, the rna6d and rna16b models are performing worst, as both are only able to display one potentially correct tree hypotheses. additionally, congruencies between the performance of rna models and the results of the aicc modeltest can be drawn from our results. according to the aicc modeltest, in the class of the rna models, the most general rna6a model is always superior to all other rna models, which is further congruent to its performance in tree reconstruction analyses. this is not reflected by the rna and rna models, where the models with the highest aicc scores. potential pitfalls of rna modelling consequences of different evolutionary constraints in stem and loop regions of rrna sequences for phylogenetic analyses has long been suspected and led to different recommendations for weighting stem positions in parsimony analyses. beside suggestions for simple one-half weighting of paired positions, empirical investigation of compensatory substitution rates in stem positions reveals a rate of about 40% of that expected under a hypothesis of perfect compensation. therefore, the weighting of stem characters is suggested to be reduced by no more than 20%. consequently, in model based tree reconstruction methods, like bayesian inference and maximum likelihood, it should be reasonable to use specific rna models as simply applying a standard dna model to data from one part of the helical regions. application of these rna models has frequently been justified by a consistent phylogenetic signal of coevolved paired sites, decreasing the information content in the data. analyses ignoring this interdependence should tend to overestimate the support for dubious or even wrong nodes in a tree. due to a reduced number of effective sites, the application of specific rna models, which take interdependencies into account, reduces tree confidence, but is more reliable in the light of the information content in the data. our results actually imply a reduced impact of stem positions in the combined data set, if mixed rna/dna model setup are used. this is depicted by the tree distance results of the separate analyses of the stem and loop partitions. in most data sets, the distances between the trees based on only the loop partition and the combined data are reduced, if rna models are applied for the combined data, whereas the distances between the stem partition and the combined data are mostly enlarged. this could have been expected, if coevolution in paired sites is assumed and thus these positions do not provide independent phylogenetic information. however, for several of the currently tested data sets, the substitution saturation test reveals that the unpaired positions clearly experience excessive homoplasy, which indicates a loss of phylogenetic information, as these positions are no longer informative. in this context, stem positions in the current data sets should contain more reliable signal, compared to loop regions, as they exhibit a much lesser grade of homoplasy due to multiple substitutions. consequently, the application of rna models increases the relative impact of noisy positions in the data set and reduces the influence of more informative portions. thus, the results of the current analyses corroborate the hypothesis proposed by letsch et al.. rna models doubtlessly provide a better depiction of the phylogenetic information content of rrna data sets, but this might be a trap, if homoplasy is far greater in loop positions. in this case, the informative phylogenetic content is obscured by noise. the situation might probably be depicted best, if we thought of a weighting scheme for rrna data sets: in standard dna model setup, the signal of paired positions is virtually weighted twice, as both positions are linked and signal of pairs can be seen redundant. as outlined above, previous studies have mostly targeted this as a problem, but the current analyses showed relative homoplasy as delimiting the confidence of the phylogenetic signal provided by loop regions, revealing a more or less hidden coherence between two factors - covariation and homoplasy - contributing to the phylogenetic signal of the rrna data sets. for this scenario it can be stated, that in contrast to a previously proposed overestimation of wrong support by ignoring site interdependencies, the application of rna models will tend to overestimate the support for dubious or wrong nodes in a tree. as depicted above, loops and stems can be expected to experience different selectional regimes, which has resulted in the development of the specific rna models. nevertheless, it has been noted as early as that substitution rates do not fit neatly into stem-loop partitions, and thus weighting according to stems vs. loops might be problematic, which was later demonstrated by van de peer. consequently, selectional constraints on rrna may not only differ between paired and unpaired regions, but also among the individual loop or stem regions, which would depend on the individual function and their relative location in the 3d rrna molecule. binding sites of ribosomal proteins, for example protein l11-binding domain within the lsu rrna domain ii and the sarcin-ricin loop within domain vi, constituting the gtpase-associated center or the lsu rrna domain v, which contains the peptidyl transferase center, are highly conserved throughout metazoa. furthermore, many of the rrna regions of domain iv that are involved in trna and inter-subunit interactions are also preserved. in contrast, the domain i of the mt lsu is highly variable on sequence level and until now, no conserved secondary structures could be detected. consequently, the partitioning into loops and stems must be seen as only an relative coarse approximation to model rrna sequences. in future phylogenetic studies on rrna, more sophisticated partitioning schemes, depending on the function, base composition and relative location of rrna regions, would be able to enhance model based tree reconstruction analyses  the results of the present study can be interpreted as a trade-off between using specific rna models for a hopefully more accurate estimation of covariation in paired sites and the risk of augmenting relatively homoplastic unpaired positions in the tree reconstruction. for future phylogenetic studies based on rrna sequences, we would therefore highly recommend a separate test for saturation of substitution in loop and stem partitions of the aligned data set. the use of a mixed rna/dna model setup should be avoided if saturation occurs in the loop partitions, as otherwise the valuable phylogenetic signal of the stem partitions might be masked by potentially noisy signal provided by the loops. in contrast, if no substantial homoplasy is detected in the data, the use of mixed rna/dna models can be highly recommend, as these lead to an increased support to probably correct clades. based on the presents results, we cannot advocate an general exclusion of the potentially noisy loop positions: first, noise is not an inherent character of a certain nucleotide position, but depends on the considered phylogenetic level. and second, differences among loop regions can be expected and excluding these regions as a whole reduces the phylogenetic signal of the data set. consequently, we rather recommend to think about enhanced partitioning strategies, which would allow a more careful modelling of rrna sequences and provide a first approach to detect noisy signal among loop partitions. however, covariation and substitution saturation are only two parameters of the evolutionary inherent pattern displayed in the data. other phenomena, like non-stationarity of substitution rates across sites and branches as well as base composition heterogeneity, might also maul the signal content of the data set. a previous study based on nuclear rrna genes, identified deviation of base composition in certain clades as probably misleading tree reconstruction analyses, rather than the covariation pattern in stem regions. a sophisticated estimation of evolutionary pattern in rrna sequence data is therefore principally desirable and newly developed methods should be applied, which are able to consider background knowledge as covariation, non-stationary processes or heterogeneity in the data. 
furthermore, the extremes of length variation among taxa, combined with regional substitution rate variation can mislead the alignment of rrna sequences and thus distort subsequent tree reconstructions. the application of specific rna models in tree reconstructions is hampered by interaction between the appropriate modelling of covarying sites in stem regions, and excessive homoplasy in some loop regions. our results depict a trade-off between increased accuracy in estimation of interdependencies in helical regions with the risk of magnifying positions lacking phylogenetic signal. however, recent developments in phylogenetic methodology now allow a comprehensive integration of secondary structures in alignment and tree reconstruction analyses based on rrna sequences, which has been shown to correct some of these problems. we found an explicit impact of rna substitution models on tree reconstruction analyses. in this context, the rna6a model outperformed all other models, including the more parametrized rna and rna models. rna models often failed to recover reasonable trees when single-stranded regions are excessively homoplastic, because these regions contribute a greater proportion of the data when covarying sites are essentially downweighted. failure to account for covariation patterns in helical regions of ribosomal rna genes has the potential to misdirect the estimation of the phylogenetic signal of the data. here, we explore the potentials of rna substitution models and the interactions of specific model setups with the inherent pattern of covariation in rrna stems and substitution rate variation among loop regions. we can therefore conclude that caution is warranted when applying rrna covariation models, and suggest that loop regions be independently screened for phylogenetic signal, and eliminated when they are indistinguishable from random noise. in addition to covariation and homoplasy, other factors, like non-stationarity of substitution rates and base compositional heterogeneity, can disrupt the signal of ribosomal rna data. all these factors dictate sophisticated estimation of evolutionary pattern in rrna data, just as other molecular data require similarly complicated corrections.
 failure to account for covariation patterns in helical regions of ribosomal rna genes has the potential to misdirect the estimation of the phylogenetic signal of the data. furthermore, the extremes of length variation among taxa, combined with regional substitution rate variation can mislead the alignment of rrna sequences and thus distort subsequent tree reconstructions. however, recent developments in phylogenetic methodology now allow a comprehensive integration of secondary structures in alignment and tree reconstruction analyses based on rrna sequences, which has been shown to correct some of these problems. here, we explore the potentials of rna substitution models and the interactions of specific model setups with the inherent pattern of covariation in rrna stems and substitution rate variation among loop regions.  we found an explicit impact of rna substitution models on tree reconstruction analyses. the application of specific rna models in tree reconstructions is hampered by interaction between the appropriate modelling of covarying sites in stem regions, and excessive homoplasy in some loop regions. rna models often failed to recover reasonable trees when single-stranded regions are excessively homoplastic, because these regions contribute a greater proportion of the data when covarying sites are essentially downweighted. in this context, the rna6a model outperformed all other models, including the more parametrized rna and rna models.  our results depict a trade-off between increased accuracy in estimation of interdependencies in helical regions with the risk of magnifying positions lacking phylogenetic signal. we can therefore conclude that caution is warranted when applying rrna covariation models, and suggest that loop regions be independently screened for phylogenetic signal, and eliminated when they are indistinguishable from random noise. in addition to covariation and homoplasy, other factors, like non-stationarity of substitution rates and base compositional heterogeneity, can disrupt the signal of ribosomal rna data. all these factors dictate sophisticated estimation of evolutionary pattern in rrna data, just as other molecular data require similarly complicated corrections. 