Preliminary phylogeny of Tribolium beetles (Coleoptera: Tenebrionidae) resolved by combined analysis of mitochondrial genes

The phylogenetic relationships of the three major species groups of Tribolium (Coleoptera: Tenebrionidae) were inferred using the simultaneous analysis of 642 bp of the most conserved part of mitochondrial DNA (mt DNA) cytochrome oxidase I (COI) and 448–452 bp of mt 16S rDNA. High sequence divergence was observed for both genes even among sibling species. The analysis of the combined segments of COI and 16S rDNA sequences produced a phylogenetic tree with moderate level of confidence. The tree topology showed monophyly of the genus Tribolium whose species were separated into three groups: “brevicornis” group (with T. brevicornis as the only representative), “castaneum” group (with T. castaneum, T. freemani, T. madens and T. audax) and “confusum” group (with T. confusum, T. anaphe and T. destructor). Sibling species pairs T. castaneum – T. freemani and T. madens – T. audax are clearly resolved. The preliminary results presented here give moderate support to the previously proposed phylogeny based on morphological data. 709 * Corresponding author; e-mail: bruvo@irb.hr of the tenebrionid genera Pimelia and Hegeter (Juan et al., 1995, 1996; Pons et al., 2004). However, this part of the COI gene appeared to be too variable for phylogenetic reconstruction within the tenebrionid genus Palorus (Meštrovi et al., 2000), a taxon closely related to the genus Tribolium. For this reason, we decided to use the most conserved part of the COI gene, proposed by Zhang & Hewitt (1997) to be suitable for lowand mid level phylogenetic analysis. The other phylogenetic marker used in this study, 16S rDNA, appears to be at least twice as conserved as COI gene in species of the order Coleoptera (Funk, 1999; Garin et al., 1999). MATERIAL AND METHODS


INTRODUCTION
The genus Tribolium Macleay, 1825 belongs to the family Tenebrionidae (order Coleoptera) and comprises 33 described species (Sokoloff, 1972), several of which are important stored product pests.Some of the species are easily cultured in the laboratory and the genus has been extensively used in genetic analyses (Sokoloff, 1972;Beeman et al., 1996;Bennett et al., 1999;Brown et al., 2000).Thanks to its convenient genetic set-up, Tribolium castaneum is predisposed to become one of the three most important insect model systems beside Drosophila and Bombyx mori (Schmidt & Tautz, 1999;Brown et al., 2003;Nagaraju & Goldsmith, 2002).
Several studies have attempted to establish evolutionary trends in these beetles.According to morphological characters, Hinton (1948) divided Tribolium species into five groups: "brevicornis" group (considered as the most primitive), "confusum", "alcine", "castaneum" and "myrmecophilum" group.According to Hinton, the first four groups separated in the Cretaceous period, while the "myrmecophilum" group separated from the "castaneum" group in the Quaternary.The "brevicornis", "castaneum" and "confusum" species-groups are predominant in the number of species and their distribution, and can be further split into sections.The remaining two groups comprise only five species that are geographically restricted to Australia and Madagascar.Within Tribolium, even the sibling species T. madens -T.audax and T. castaneum -T.freemani, able to hybridize, are morphologically distinguishable by numerous characters pointing to a high species divergence in the group (Halstead, 1969;Wade & Johnson 1994).
Cladistic relationships among some Tribolium species inferred from cytological studies of chromosome numbers, revealed the "castaneum" group as more "primitive" than the "confusum" group (Smith, 1952).Evolutionary relationships of the most common Tribolium species were also studied by electrophoretic analysis of isozymes (Wool, 1982).According to this study, T. castaneum and T. confusum are more closely related to T. brevicornis than they are to each other, which is in disagreement with cytological and morphological data.A comparison of the defensive secretions from Tribolium species indicated that "brevicornis" and "castaneum" are sister species groups, and that T. brevicornis and T. audax are members of the same species group (Howard, 1987).
To our knowledge, there has been no study of this genus so far that would be based on phylogenetically informative molecular data.At the DNA level, relationships between some Tribolium species were inferred using characteristics of highly repetitive DNAs, their sequence homology and abundance (Juan et al., 1993;Ugarkovi et al., 1996).The cladogram obtained from satellite DNA hybridization data, karyological data and chemical compounds of secretion (Juan et al., 1993) supported the "castaneum" species-group, although T. audax was placed outside of this group, while "confusum" and "brevicornis" species-groups appeared as sister clades, in contrast to the phylogeny proposed by Hinton (1948) andHoward (1987).
In this work we address the phylogenetic relationships among the eight most studied species of Tribolium: T. castaneum, T. confusum, T. destructor, T. madens, T. audax, T. anaphe, T. freemani and T. brevicornis, which are also the most widespread and economically important.These eight species represent, according to Hinton (1948), three major species groups: the "castaneum" group (T.castaneum, T. audax, T. freemani and T. madens), the "confusum" group (T.confusum, T. anaphe, T. destructor) and the "brevicornis" group (T.brevicornis).Despite our efforts, we were unable to get hold of any of the additional "brevicornis" group species, nor any of the species belonging to "alcine" and "myrmecophilum" groups, and therefore we regard this study as an initial and preliminary analysis of Tribolium phylogeny.In order to resolve the relationships of these eight species we used two genes that are regarded as the most conserved among mitochondrial markers routinely used for phylogenetic inference: cytochrome oxidase I (COI) and 16S rDNA.COI has been used for phylogenetic inference of many insects (Lunt et al., 1996) -the variable part of this gene has been used to deduce phylogenetic relationships among congeneric species of the tenebrionid genera Pimelia and Hegeter (Juan et al., 1995(Juan et al., , 1996;;Pons et al., 2004).However, this part of the COI gene appeared to be too variable for phylogenetic reconstruction within the tenebrionid genus Palorus (Meštrovi et al., 2000), a taxon closely related to the genus Tribolium.For this reason, we decided to use the most conserved part of the COI gene, proposed by Zhang & Hewitt (1997) to be suitable for low-and mid level phylogenetic analysis.The other phylogenetic marker used in this study, 16S rDNA, appears to be at least twice as conserved as COI gene in species of the order Coleoptera (Funk, 1999;Garin et al., 1999).

Beetle material
Insect cultures (laboratory strains) used in this study were obtained from the Central Science Laboratory (Sand Hutton, York, UK).A natural population of the species T. madens was collected in a flourmill (Virovitica, northern Croatia).Total DNA was extracted from a single individual according to Steiner et al. (1995), or from pooled individuals following the standard protocols.

Phylogenetic analyses
The alignment of the COI and 16S rDNA sequences was done with CLUSTAL X version 1.81 (Thompson et al., 1997).Bio-Edit 5.09 (Hall, 1999) was used for manual inspection and correction of 16S rRNA alignments according to secondary structure of Drosophila yakuba 16S rRNA (Gutell & Fox, 1988).Aligned sequences were reformatted appropriately and entered into MEGA 2.1 (Kumar et al., 2001) for calculation of sequence statistics.In order to downweight the gaps present in 16S ribosomal gene (24 gaped positions altogether), we applied the "simple method" for recoding gaps described in Simmons & Ochoterena (2000) and Simmons et al. (2001) using the program GapCoder (Young & Healy, 2002).The matrices obtained this way were then used in both parsimony and Bayesian searches.The alignments used in phylogenetic analyses are available on request from the authors.
COI and 16S rDNA data sets were tested using the partition homogeneity test (Farris et al., 1994) implemented within PAUP to assess whether combining of COI and 16S sequences is appropriate for phylogenetic analysis.Phylogenetic analyses were performed on separate and combined datasets by parsimony (MP) and maximum likelihood (ML) searches using PAUP 4.0b10 (Swofford, 1998), and by Bayesian analysis with Markov Chain Monte Carlo (MCMC) sampling in MrBayes 3.1 (Huelsenbeck & Ronquist, 2001).
For MP analyses, we used the heuristic search with treebisection-reconnection (TBR) branch-swapping option.The extent of support for internal nodes was estimated using the bootstrap method (1000 pseudoreplicates) (Felsenstein, 1985).Bremer support and partitioned Bremer support values (Baker & DeSalle, 1997) were calculated based on strict consensus MP trees using the program TreeRot 2 (Sorenson, 1999).
The model of DNA substitution was determined for each gene with the program MODELTEST 3.04 (Posada & Crandall, 1998).The estimated parameters (GTR+G+I, being the model preferred by both Akaike information criterion (AIC) and hierarchical likelihood ratio test (hLRT) for both partitions) were then invoked for ML heuristic searches with simple sequence addition and TBR branch-swapping option.
Bayesian MCMC analyses with four parallel Markov chains were run under the GTR+G+I model (unlinked parameters for each partition, overall rate allowed to be different across partitions) for 1000000 generations, with sampling frequency of one in every hundred trees.The consensus tree (with posterior probabilities for each node as a measure for support) was constructed based on the trees sampled after the convergence of likelihood scores (burn-in).

Sequence variation of COI gene
Sequence analysis of a 642 bp fragment of the most conserved part of mitochondrial COI gene was performed in eight Tribolium species as well as in two outgroup species, Pimelia elevata and Tenebrio molitor.Intraspecific variation, assessed by screening four to eight specimens of each species, ranges from 0 to 1%, with the exception of T. freemani where intraspecific variability is 2.3% mainly due to the slippage in one of the four analysed clones.In addition, one of T. brevicornis sequences appeared shorter due to several internal deletions leading to frame-shift mutations.It can be suspected that such a sequence represents a nuclear paralogue of a part of COI gene, and for this reason it was not considered for further analyses.Other COI fragments did not present any peculiarities that would point to their nuclear provenience (review on mitochondrial pseudogenes in Bensasson et al., 2001).In the case of T. madens, the specimens from two diferent populations, laboratory and natural, were analyzed.Differences between these populations (0.06-0.50%) are within the range of intraspecies divergences of other examined species.
Among the 642 examined nucleotide positions of tested Tribolium sequences, 223 positions were polymorphic (34.7%) and 153 (23.8%) phylogenetically informative.Nucleotide composition, number of variable sites and Ti/Tv ratios are presented in Table 1.The sequences have an average frequency of A+T nucleotides of 63.0%, increasing to 74.3% at third codon positions.Similar frequency of A+T has been observed in two other tenebrionid genera -Pimelia (Juan et al., 1995) and Hegeter (Juan et al., 1996).The observed A+T bias is lower than in other insects, what is suggested to be a general characteristic of beetle COI/COII genes (Liu & Beckenbach, 1992;Funk et al., 1995).The average transition/transversion ratio (calculated from the alignment of the entire dataset) in examined sequences is 1.3, being the highest at first codon positions, while dropping close to 1 at third codon positions.Of the 223 variable positions, 80.7% are at third codon positions, 15.7% at first, while the second codon positions are the most conserved, with only 3.6% of nucleotide changes (Table 1).
Comparison of inferred amino acid (aa) sequences of COI in Tribolium using mtDNA code for D. melanogaster (De Bruijin, 1983) showed that the 214 aa long polypeptide is highly conserved with only 20 changes (10%) of which 12 are conservative.Most of these changes are located in the variable part of the protein between positions 182-212, what is in accordance with the suggested aa variability of insect COI (Lunt et al., 1996).Maximum aa divergence between pairs of sequences within Tribolium amounts to 7%.
Uncorrected p-distances were calculated for each species pair.The minimum distance of 12% is calculated for T. confusum-T.anaphe species pair, while all other sequence pairs are within the range from 14% to the maximum divergence of 19% for T. destructor-T.audax.The genetic distances between Tribolium and the outgroups (Tenebrio molitor and Pimelia monticola) range between 16% and 20%.The fact that some ingroup divergences are higher than the ingroup-outgroup distance, together with low Ts/Tv ratio, and high number of Tv A C as well as variable sites in the 3 rd codon position (80.7%),suggests that there is a high level of saturation due to multiple changes in the COI data set.

Sequence variation of 16S rDNA
Sequence analysis of a 448-452 bp fragment of the 16S rDNA has been performed in the same species used to obtain COI data.Intraspecific variation was in the range from 0 to 0.89%.The fragment length was conserved within the species, but in the alignment of T. freemani clones two gaps are created, although the clones have the same length.Since analyzed 16S rDNA sequences showed unexpectedly high interspecific divergence (see below), they could not be easily aligned with a significant confidence.The computer-generated alignment was therefore further adjusted considering published secondary structures for Drosophila yakuba 16S rRNA (Gutell & Fox, 1988).Some stems show noncanonical G-T pairs, which also help in maintaining stem stability (Gutell et al., 1994).In the highly variable region between positions 190-249, which corresponds to G3 helix (De Rijk et al., 1996), it was not possible to predict the stem and loop regions.
Sequence analysis reveals significantly higher overall A+T content in 16S than in COI (Table 1).The base composition is clearly biased towards AT, especially in the loops.The percent of variable sites is similar between stems and loops, consistent with previous observations that both regions contain sites allowed to evolve with heterogeneous rates regardless of their secondary structure (Simon et al., 1994).If G3 helix is excluded from the analysis, the percentage of variable sites is lower in 16S than in COI, but when G3 is considered, the divergences are highly similar.Since it is expected that insect 16S rDNA evolves at least twice slower than COI (Funk, 1999), this suggests saturation of 16S rRNA gene as well.The overall Ts/Tv ratio of 16S is twice as low as that of COI, as expected for 16S and COI sequences showing similar levels of divergence (Funk et al., 1995).The 16S Ts/Tv ratio was twice higher in stems (similar to that of COI) than in loops, in which Ts may be completely saturated.
Within the examined Tribolium species, genetic distances reveal the minimum value of 7% for the T. madens -T.audax pair and 9% for the T. anaphe -T.confusum pair, and the maximum value of 18% for the T. freemani -T.destructor pair.Genetic distances between Tribolium and the Tenebrio / Pimelia outgroups are in the range from 16 to 22%.

Phylogenetic analysis COI
COI trees obtained by all three methods revealed the same topology (not shown).Each of the species is clearly monophyletic, species cluster in species-groups as proposed by morphology, and within each of the species groups the relations among the species are well supported.In addition, COI supports the hypothesis of the sister-group relationship between "castaneum" and "confusum" groups.
We also analyzed the effects of several weighting schemes; all produced the same tree topology, so we proceeded with the analyses by using simultaneously two weighting methodologies, substitution weighting 2 : 1 for Tv/Ts and character weighting 2 : 5 : 1 for the first vs.second vs. third codon positions, being one of the most commonly used weighting schemes.Ignoring transitions and substitutions on third codon positions resulted in low phylogenetic resolution, indicating that such an approach results in a loss of information.

16S rDNA
The trees obtained by all three methods for 16S rDNA gene revealed fairly similar topologies (not shown).The major uncertainty concerned the position of T. brevicornis, which was separated as a species group or joined to the "confusum" clade.16S is phylogenetically most suitable for resolving relationships among taxa exhibiting intermediate levels of COI divergence (~5-15%) (Funk et al., 1995;Funk, 1999).However, COI divergence between species of Tribolium was 12-19%, and the 16S marker alone turned out to be insufficient to resolve unequivocally the phylogeny of the genus.

Combined COI and 16S rDNA
Since separate analysis of mitochondrial COI and 16S genes could not completely resolve the phylogenetic relationships among Tribolium species under study, further improvement was obtained by the analysis of combined 642 bp long COI fragment plus 448-452 bp long 16S fragment.The ILD test (Farris et al., 1994)  gesting that the two data sets are combinable.Partitioned Bremer support values (Fig. 1a) indicate that the COI partition makes a slightly higher contribution to the overall support for the nodes than does the 16S partition.
The trees obtained in all three MP (Fig. 1a), ML (not shown) and Bayesian (Fig. 1b) analyses of combined sequences clearly showed monophyly of the genus Tribolium and also of each of the species."Castaneum" and "confusum" clades are well supported, and within these two groups the relationships between the species as well as the two sibling species pairs (T.freemani -T.castaneum and T. madens -T.audax) are clearly resolved.Most nodes have high confidence values, except for the position of T. brevicornis.This species is positioned as a sister taxon to all other Tribolium species with a moderate support in MP and ML trees, while in the Bayesian tree it appears as a sister taxon to "castaneum" clade (though with low support).

DISCUSSION
Despite its great economical importance, the Tenebrionid genus Tribolium has not been extensively studied so far on the level of informational macromolecules, except for highly repetitive DNA traits of the several stored-product pest species (Juan et al., 1993;Ugarkovi et al., 1996;Mravinac et al., 2004Mravinac et al., , 2005)).The reason for this neglection most probably resides in fact that the great majority of Tribolium species are free-living, many of them are very rare and/or inhabit remote and not easily reachable areas, and, to our knowledge, very few people might have (or, have had) the access to these species.This is especially so for the species belonging to "myrmecophilum" and "alcine" groups (restricted to Australia and Madagascar, respectively), but also for the "brevicornis" group species with the exception of T. brevicornis.Our original intention was to collect and include in our phylogenetic analysis the representatives of all Tribolium species-groups, but, unfortunately, this task turned out to be impossible.However, despite of this initial failure, we 712  (Hinton, 1948;Halstead, 1969) are indicated with asterisks, and proposed species groups with their geographic regions of origin (Hinton, 1948)  decided to conduct this study in order to establish a preliminary phylogeny that could in turn be expanded and completed by the addition of the missing taxa.
In our study we included eight Tribolium species belonging to three species groups -"confusum", "castaneum" and "brevicornis".In order to resolve the phylogenetic relationships of these species we analyzed the most conserved part of the COI gene and examined part of 16S rDNA which evolves at least twice as slow as COI and is suitable for analysis of divergent species within a genus, differentiating among genera and families (Funk, 1999).Separate analyses of these two mitochondrial genes resulted in similar classification of Tribolium species, but, taken independently, they were not sufficient to resolve the phylogenetic relationships of the eight species under study.The phylogenetic relationships among examined Tribolium species were therefore analyzed using the combined segments of COI and 16S rDNA genes.Such an analysis showed monophyly of the genus Tribolium and clearly resolved species in three groups (Fig. 1), mostly in the same manner as they have been clustered according to the examination of morphological characters (Hinton, 1948).In the absence of the "myrmecophilum" and "alcine"-group species it is obviously not possible to confirm the sister-group relationships which would involve these species-groups, but the results of MP and ML suggest, with a moderate support, that the "brevicornis" group is probably the sister group to the other Tribolium species clade.However, since the Bayesian analysis suggested a slight possibility of a closer relationship of T. brevicornis with the "castaneum" group, we regard that the exact position of T. brevicornis could not be unequivocally resolved.We believe that the addition of species belonging to "brevicornis" group would greatly enhance its resolution.
Phylogenetic analysis of thirty-seven coleopteran species based on the mt COI gene revealed very long branches for species belonging to the families Scarabaeidae and Tenebrionidae, indicating high divergence among members of these families (Howland & Hewitt, 1995).Although the present study involved relatively conserved segments of COI and 16S rDNA mitochondrial genes, the sequence divergence was surprisingly high among examined Tribolium species, in a range from 12 to 19% for COI and from 7 to 18% for 16S.Low intraspecific variability is observed for both tested genes in all analyzed Tribolium species, including laboratory and wild populations, which can be almost exclusively found in storehouses (Wool, 1982).
Phylogenetic analysis of the collembolan genus Orchesella (Frati et al., 2000) and coleopteran genera Neochlamisus, Ophraella, Timarcha and Chrysolina (Funk et al., 1995;Funk, 1999;Garin et al., 1999;Gómez-Zurita et al., 2000) revealed intrageneric divergences in a high range from approx.0.7% to 21% for COI and from 0.2% to 12% for 16S rDNA.In contrast to the insect species mentioned above, observed divergences are dramatically pronounced even among very close relatives within Tribolium species: siblings T. madens -T.audax and T. freemani -T.castaneum.For these species pairs, observed divergence is 14% and 15% for COI, and 7% and 13.5% for 16S rDNA, respectively, being within the level of divergence that has been usually observed between species groups in coleopteran genera (Funk et al., 1999).In species of the genus Pissodes (Coleoptera), which can hybridize producing fertile hybrids, the divergence of COI ranges from 0.5% to 7.5% (Langor & Sperling, 1997).Although COI in Tribolium showed higher sequence divergence than observed in other coleopteran species (e.g., Funk et al., 1995;Funk, 1999;Garin et al., 1999;Gómez-Zurita et al., 2000), the maximum divergence of COI at the amino acid level within Tribolium amounts to 7%, being within the divergence of 11% predicted for species within cole-opteran genera.This is in agreement with the notion that COI is the most slowly evolving among mitochondrial protein coding genes at the amino acid level (Funk, 1999).
Based on the facts presented here, we assume that COI and 16S rRNA genes are not sufficient to completely resolve the phylogeny of the genus Tribolium.Future work on this genus should therefore include additional molecular markers, preferably some of the nuclear genes or genomic regions that proved useful for investigations of sub-familial interrelationships of insects, such as ribosomal RNA genes, elongation factors or ITS regions (e.g.Danforth et al., 2005;Gómez-Zurita & Galián, 2005;Gómez-Zurita et al., 2005;Becerra, 2004;Whitfield et al., 2002;Rokas et al., 2002;Cryan et al., 2001;Clark et al., 2001;Krzywinski et al., 2001;Mardulyn & Whitfield, 1999).
In conclusion, our study gives moderate support to the topology of the phylogenetic tree based on morphological traits proposed by Hinton (1948).We strongly hope that this preliminary analysis of phylogenetic relationships within this genus will serve as a stimulus for the Tribolium-oriented community to initiate a more thorough taxon sampling of wild-living species of Tribolium in order to fill the currently existing gaps and to reveal the "true" phylogeny of this important beetle-genus.

Fig. 1 .
Fig. 1.Phylogenetic relationships among eight examined Tribolium species proposed on the basis of combined analysis of COI and 16S data sets.Acronyms "p" and "h" in the names of clones denote population and haplotype numbers, respectively.(a) The strict consensus MP tree (length 1178, CI = 0.610, RI = 0.852); confidence levels of tree bifurcation were assessed by 1000 bootstrap replicates.Bootstrap values are indicated above the nodes, while the partitioned Bremer support values are given below the nodes (COI/16S).Two pairs of sibling species(Hinton, 1948;Halstead, 1969) are indicated with asterisks, and proposed species groups with their geographic regions of origin(Hinton, 1948) are indicated with vertical bars.(b) Majority-rule consensus tree resulting from Bayesian MCMC analysis.Numbers on the nodes are posterior probability values.Branch lengths are proportional to number of changes.
Fig. 1.Phylogenetic relationships among eight examined Tribolium species proposed on the basis of combined analysis of COI and 16S data sets.Acronyms "p" and "h" in the names of clones denote population and haplotype numbers, respectively.(a) The strict consensus MP tree (length 1178, CI = 0.610, RI = 0.852); confidence levels of tree bifurcation were assessed by 1000 bootstrap replicates.Bootstrap values are indicated above the nodes, while the partitioned Bremer support values are given below the nodes (COI/16S).Two pairs of sibling species(Hinton, 1948;Halstead, 1969) are indicated with asterisks, and proposed species groups with their geographic regions of origin(Hinton, 1948) are indicated with vertical bars.(b) Majority-rule consensus tree resulting from Bayesian MCMC analysis.Numbers on the nodes are posterior probability values.Branch lengths are proportional to number of changes.

TABLE 1 .
rendered a P value of 0.168, which is nonsignificant, sug-Summary of nucleotide composition, number of variable sites and Ti/Tv ratio for COI and 16S rDNA data of examined Tribolium beetles.