Introduction

The classical view emerging from anatomical and archeological studies places Neandertals as a different species from Homo sapiens. This is in agreement with the Out-of-Africa hypothesis (Cann, Stoneking, and Wilson 1987<$REFLINK> ), which predicts that Neandertals coexisted without mating with modern humans who originated in Africa from 100,000 to 200,000 years ago (Stringer and Andrews 1988<$REFLINK> ). Instead, recent anatomical and paleontological research (Wolpoff et al. 2001<$REFLINK> ) supports the multiregional hypothesis, which propounds that some populations of archaic Homo evolved into modern human populations in many regions. Consequently, Neandertals could have contributed to the genetic pool of present-day Europeans.

Studies about modern human genetic diversity (Foley 1998<$REFLINK> ; Jorde et al. 2000<$REFLINK> ) assume that Neandertals were not related to modern humans, although such assumptions have been extensively debated (Nordborg 1998<$REFLINK> ; Hawks and Wolpoff 2001<$REFLINK> ; Relethford 2001<$REFLINK> ).

Five mitochondrial DNA D-loops have been recovered from Neandertal bones (Krings et al. 1997<$REFLINK> ; Krings et al. 1999<$REFLINK> ; Krings et al. 2000<$REFLINK> ; Ovchinnikov et al. 2000<$REFLINK> ). Henceforth, we will refer to these research groups as the Neandertal sequencing groups (NSGs). The phylogenetic analyses of these sequences located Neandertal DNA at the base of modern human diversity (fig. 2C ), suggesting that the Neandertal genes probably did not contribute to the modern human genetic pool.

However, the human D-loop region shows an extreme variation in the substitution rate among sites (Excoffier and Yang 1999<$REFLINK> ; Meyer, Weiss, and von Haeseler 1999<$REFLINK> ) and a high amount of parallel mutations (Tamura and Nei 1993<$REFLINK> ) that complicate phylogenetic reconstructions (Maddison, Ruvolo and Swofford 1992<$REFLINK> ; Ingman et al. 2000<$REFLINK> ). Moreover, from a phylogenetic point of view, the NSG neither selected for the best model of nucleotide substitution nor estimated the model parameters (table 1 ). An abundant literature (Yang, Goldman, and Friday 1994<$REFLINK> ; Swofford et al. 1996<$REFLINK> ; Huelsenbeck and Rannala 1997<$REFLINK> ; Brocchieri 2001<$REFLINK> ; Posada and Crandall 2001<$REFLINK> ) emphasizes the need for an accurate estimation of the nucleotide substitution model and parameters, correction for among-site rate variation, taxa sampling, and selection of the outgroup to attain a correct tree reconstruction.

The purpose of this article is to reanalyze the available Neandertal DNA sequence data. After selecting the appropriate extant human D-loop data sets, we used the maximum-likelihood principle for model selection and parameter estimation to perform a novel phylogenetic reconstruction.

Data and Methods

Sequences and Data Sets

We have used the MOUSE 1.0 database (Burchardt, von Haeseler, and Meyer 1999<$REFLINK> ), an aligned compilation of mtDNA control regions of primate species that contained the Feldhofer Neandertal sequences. An earlier version of this database (HvrBase) was used by the NSG. Other Neandertal sequences were retrieved from GenBank (table 1 ) and manually aligned to our MOUSE-selected data sets.

We carried out two analyses: the first analysis concerns the HVI region and includes three Neandertal sequences. The second analysis concerns the HVI plus HVII regions including two Neandertal sequences. The HVII region of the Mezmaiskaya specimen has not been recovered.

In the first analysis (HVI region) we included human entries containing at least positions 16056 to 16378 of the Anderson et al. (1981)<$REFLINK> reference sequence (fig. 1 ). This segment covers the minimal region common to the three Neandertal entries. Human sequences, duplicated or ambiguous, were deleted. This selection rendered 1,905 human sequences (328 Africans, 471 Asians, 211 Australians-Oceanians, 475 Europeans, and 420 Americans). Three Pan paniscus and three Pan troglodytes sequences were added as outgroup. Also, we discarded alignment columns where one sequence showed an insertion, whereas the rest did not show any insertion. The final alignment contains 3 Neandertal, 1,905 human, and 6 chimpanzee sequences. The range used for analysis covers positions 16023–16400, that correspond to the longest Neandertal sequence (fig. 1 ). In the second analysis we included individuals where both HVI and HVII regions were sequenced at least from positions 16056–16378 and 57–343, respectively (fig. 1 ). These ranges cover the region common to the two Neandertal sequences. Repeated or ambiguous entries were deleted. A total of 377 human entries was selected (79 Africans, 59 Asians, 84 Australians-Oceanians, 123 Europeans, and 32 Americans). The final alignment contains 2 Neandertal, 377 human, and 6 chimpanzee sequences that cover positions 16023–16400 and 57–343. We deleted alignment columns where most sequences showed a gap.

Model Testing

In order to select the best model of nucleotide substitution fitting each sequence set, we have followed the procedure described by Posada and Crandall (1998)<$REFLINK> , implemented in the MODELTEST program. This approach first obtains a neighbor-joining tree (Saitou and Nei 1987<$REFLINK> ), using the Jukes-Cantor model of nucleotide substitution (Jukes and Cantor 1969<$REFLINK> ) and then uses a likelihood ratio test statistic to select the best model and estimate its parameters, keeping the same tree topology.

Pairwise Distances

Uncorrected and maximum-likelihood pairwise distances have been computed using PAUP* 4.0 (Swofford 1998<$REFLINK> ). Maximum-likelihood distances were calculated using the best model fitting the data.

Phylogenetic Analysis and Bootstrapping

Phylogenetic analyses were carried out using PAUP* 4.0 (Swofford 1998<$REFLINK> ). For each data set the procedure was: (1) the sequence order was randomized, (2) the maximum-likelihood distances (using the selected model and the best parameters) were computed, and (3) a neighbor-joining tree (Saitou and Nei 1987<$REFLINK> ) was constructed for 100 bootstrap replicates (Felsenstein 1985<$REFLINK> ). Interior branch tests (Sitnikova, Rzhetsky, and Nei 1995<$REFLINK> ) were carried out with the program PHYLTEST 2.0 (Kumar 1996<$REFLINK> ) using the most complex model available in this program: the Kimura two-parameter model (Kimura 1980<$REFLINK> ) with gamma correction.

Results and Discussion

Selection of Best Model of Nucleotide Substitution, and Maximum-Likelihood Estimation of the Model Parameters

The best model fitting the HVI and HVI + HVII sets is the Tamura-Nei with a gamma distribution for rate heterogeneity and a proportion of invariable sites, out of the 56 different models implemented in the MODELTEST program. Actually, this model was developed for the primate D-loop region (Tamura and Nei 1993<$REFLINK> ). The advantage of applying a correction for rate variation among sites in this region has also been reported (Tamura and Nei 1993<$REFLINK> ; Wakeley 1993<$REFLINK> ; Excoffier and Yang 1999<$REFLINK> ; Meyer, Weiss, and von Haeseler 1999<$REFLINK> ). It has been shown that correction for rate variation among sites and consideration of invariable sites improve phylogenetic reconstruction in cases where a long branch attraction (LBA) artifact exists (Philippe and Laurent 1998<$REFLINK> ). As we will show later, chimpanzee sequences can attract the Neandertal ones.

The parameters estimated in our work (table 2 ) are similar to those obtained by other authors (Excoffier and Yang 1999<$REFLINK> ; Meyer, Weiss, and von Haeseler 1999<$REFLINK> ). There is a strong difference between transitions and transversions, and a low value for the gamma distribution shape parameter, which indicates a strong heterogeneity in the substitution rate among sites. The value of the shape parameter estimated here is greater than other estimations (Excoffier and Yang 1999<$REFLINK> ; Meyer, Weiss, and von Haeseler 1999<$REFLINK> ) because we have introduced a parameter for the proportion of invariable sites.

Phylogenetic Reconstruction and Bootstrap

Schematic representations of the bootstrap tree topologies derived from our HVI and HVI + HVII data sets are shown in figure 2A and B . The detailed topologies, including 1,914 and 385 sequences, respectively, are not shown here, but are available upon request. Bootstrap values were not significant for most nodes. Although the three Neandertal sequences clustered together, we obtained no support for a branch separating the Neandertal cluster from the human sequences. Our HVI tree (fig. 2A ) is similar to the NSG one, but branch support is not significant. The HVI + HVII tree (fig. 2B ) places 10 Africans sequences as outgroup of Neandertals and other humans, but again, the bootstrap value is not significant. Therefore, our results suggest that the Neandertal sequences cannot be considered an outgroup for modern humans. This result is in contrast to the tree obtained by the NSG, where the Neandertal sequences appear basal to modern humans when using chimpanzee sequences as outgroup (fig. 2C ).

Given that neither the bootstrap method (Felsenstein 1985<$REFLINK> ) nor the likelihood mapping method used by Krings et al. (1997<$REFLINK> , 1999<$REFLINK> , 2000)<$REFLINK> have statistical meaning, we performed the interior branch test. It has been shown that the confidence value of this test is the complement of the P values in the standard statistical test (Sitnikova, Rzhetsky, and Nei 1995<$REFLINK> ). We divided our data sets into four clusters: (1) chimpanzees, (2) Neandertals, (3) extant human sequences located at the base of the HVI + HVII tree (fig. 2B ) that correspond to 10 !Kung entries of the MOUSE database, and (4) the rest of humans. We then determined whether the lengths of the interior branches of the three possible trees for the four sequence clusters are significantly different from 0. The result of this test is not significant for either the HVI or HVI + HVII regions. Thus, like the bootstrap analyses, the interior branch test suggests that a polytomy is the best representation for the evolutionary tree relating Neandertal and extant humans.

Pairwise Distances

The NSG reported that the pairwise comparisons between the Neandertal and human sequences demonstrate that Neandertals are outside of modern human D-loop variability. In particular, Krings et al. (1997)<$REFLINK> stated that ‘a total of 0.002% of the pairwise comparisons between human mtDNA sequences were larger than the smallest difference between the Neandertal and the humans.’ We think that this point merits further analysis. The current database is biased because of the overrepresentation of some populations and the underrepresentation of others. For instance, the MOUSE database contains 6,012 entries for the HVI region, but 31% of the entries belong to only 20 populations out of 206 populations represented (10% of the total of populations). The extreme cases are 306 Koreans, 126 Yaps, 120 Cayapa Amerindians, 119 Mandeka, 115 Palau, and 100 white British. There are also 1,417 entries of undetermined population (40% of them are from North America and 23% European, but only 9% are from Africa). Thus, African populations containing the most ancient lineages and the highest variation are underrepresented in the database.

Because of the database overrepresentation of some human populations, the distribution of pairwise distances is biased. A large part of pairwise comparisons are made between individuals belonging to the same population. Likewise, it is expected that most individuals from a single population will show similar distances to a given outgroup (Neandertal, in this case). To overcome this problem, we considered another sample of the human variation. We first sorted the HVI sequences in our data set according to its uncorrected distance to the reference sequence (Anderson et al. 1981<$REFLINK> ). Then we grouped them into 171 classes, containing equidistant sequences (considering four decimals), and chose one sequence at random from each class. The computation of pairwise distances between the 171 randomly selected sequences and the Neandertals rendered 1.6% of human-human comparisons larger than the smallest difference between Neandertals and humans. Likewise, 27% of the comparisons are lower than the largest human-human difference. This result suggests that Neandertals sequences are not so different from those of extant humans, in contrast to the NSG claims.

Final Considerations

A main conclusion can be extracted from our analyses: the phylogenetic position of the ancient DNA sequences recovered from Neandertal bones is sensitive to the phylogenetic methods employed. It depends on the model of nucleotide substitution, the branch support method, and the set of data used. Adcock et al. (2001)<$REFLINK> recovered HVI sequences of archaic human bones from Australia, and their phylogenetic analysis showed that two of the specimens were outgroups even for the most ancient African lineages. They concluded that this is an evidence supporting the multiregional hypothesis. However, a second analysis carried out by Cooper et al. (2001)<$REFLINK> that took into account the heterogeneity of rates between sites and a large sample of modern humans, showed that both HVI sequences are located among extant humans. This case illustrates the influence of the nucleotide substitution model on the phylogenetic reconstruction of the human D-loop region.

The NSG studies used poor parameter models of nucleotide substitution for their analyses, whereas we opted for complex (parameter rich) models following the likelihood ratio test. Both alternatives have pros and cons. Yang (1997)<$REFLINK> and Posada and Crandall (2001)<$REFLINK> encourage the use of the best-fit model, although they reported some examples where simple models can recover the true phylogenetic tree better than complex models. Several authors highlight the importance of using the best nucleotide substitution model for a given data set (Huelsenbeck and Hillis 1993<$REFLINK> ; Sullivan and Swofford 1997<$REFLINK> ; Cunningham, Zhu, and Hillis 1998<$REFLINK> ). The same considerations apply to the choice of the branch support method. The bootstrap method has some caveats (Cummings, Otto, and Wakeley 1995<$REFLINK> ), and the quartet puzzling (likelihood mapping) has also been criticized for overestimating branch support values (Cao, Adachi, and Hasegawa 1998<$REFLINK> ).

Secondly, we believe that the likelihood mapping values supporting Neandertals as a different species might be artifactually increased. To illustrate this point, we took two human HVI sequences showing maximal divergence between them, a chimpanzee, and the Feldhofer Neandertal. The three possible quartets were evaluated with the program TREE-PUZZLE (Strimmer and von Haeseler 1996<$REFLINK> ) with the same model and parameter used by Krings et al. (1997)<$REFLINK> . It can be seen (fig. 3 ) that the chimpanzee branch is very long, and that the Neandertal branch is longer than the two humans, opening the possibility of a LBA artifact. It has been reported that quartet methods are very prone to LBA (Huelsenbeck 1998<$REFLINK> ). Thus, the best-scored topology is usually the one joining long branches in exhaustive quartet evaluations. Another issue relates to the significance of the log-likelihood values. In the trees drawn in figure 3 , the best score is for the tree that joins Neandertal and chimpanzee sequences. However, the log-likelihoods are very similar between the alternative topologies. The Kishino-Hasegawa test (Kishino and Hasegawa 1989<$REFLINK> ) yielded no significant differences between them at the 5% level. This example illustrates that likelihood mapping can select as the best topology one that is not significantly better than the other two. A third argument against the use of likelihood mapping in this study has to do with the database population bias mentioned earlier. It is not advisable to perform likelihood mapping with a database where most human sequences are very similar among themselves.

Ovchinnikov et al. (2000)<$REFLINK> also obtained a high bootstrap value (95%) for the Neandertal clade, segregating them from modern humans. When we used the same phylogenetic procedures, i.e., corrected Tamura-Nei distances with a gamma distribution of 0.4, we obtained a lower bootstrap value (62%). An explanation for this discrepancy could be that Ovchinikov et al. did not select human sequences by length and did not eliminate sequences with ambiguities. Missing data are a common nuisance factor in phylogenetic analysis (Wiens 1998<$REFLINK> ). Although we obtained a lower bootstrap value, we think that 62% is still a high bootstrap value compared with the ones obtained when maximum-likelihood distances are used (fig. 2A and B ). It has been shown that maximum-likelihood distances have some advantages over corrected distances, reducing sampling variances of the transition-transversion ratio estimations (Swofford et al. 1996<$REFLINK> ). Meyer, Weiss, and von Haeseler (1999)<$REFLINK> estimated a transition-transversion ratio of 15.55 for the HVI region. This points to an underestimation by Ovchinnikov et al. (2000)<$REFLINK> of the number of substitutions along the tree branches.

The low bootstrap values obtained in our work indicate the lack of resolution of the HVI and HVII sequences to determine the phylogenetic position of Neandertals because no tree topology is favored. In part, this is caused by short sequences with a high proportion of invariable sites. Ingman et al. (2000)<$REFLINK> noted that the D-loop region, in spite of its wide use for human phylogenetic analysis, is not an appropriate region of the mitochondrial genome for answering phylogenetic questions.

Another interesting issue is the possibility that the Neandertal sequences were artifactual. Caldararo and Gabow (2000)<$REFLINK> noticed that some of the nucleotide substitutions found in the Feldhofer Neandertal sequence matched those found in other human ancient HVI sequences. They attributed these similarities to diagenetic changes that occurred in fossil DNA. Hansen et al. (2001)<$REFLINK> studied different clones of ancient DNA amplifications and concluded that PCR can introduce errors caused by the miscoding lesions of ancient DNA. In the absence of diagenetic changes, the mean distance between an ancient sequence and an extant one ought to be larger than the distance between a more ancient sequence and the extant one. So, the more ancient a sequence, the shorter its distance to the extant sequences (fig. 4 ). We tested this prediction by computing the average distances between each Neandertal HVI sequence and our data set of extant HVI human sequences (table 3 ). It can be seen that the most recent specimen (Mezmaiskaya) shows a shorter distance than the oldest one (Feldhofer). This result can be explained by a different nucleotide substitution rate among Neandertal lineages (populations) or by miscoding DNA lesions in Neandertal fossils. Had DNA damage increased the differences between Neandertal and modern humans, Neandertals would be more akin to modern humans than what recent claims suggest.