The Neighbor-joining Method: A New Method for
Reconstructing Phylogenetic Trees’
Naruya Saitou2 and Masatoshi Nei
Center for Demographic and Population Genetics, The University of Texas
Health Science Center at Houston

A new method called the neighbor-joining method is proposed for reconstructing
phylogenetic trees from evolutionary distance data. The principle of this method
is to find pairs of operational taxonomic units (OTUs [ =neighbors]) that minimize
the total branch length at each stage of clustering of OTUs starting with a starlike
tree. The branch lengths as well as the topology of a parsimonious tree can quickly
be obtained by using this method. Using computer simulation, we studied theefficiency of this method in obtaining the correct unrooted tree in comparison withthat of five other tree-making methods: the unweighted pair group method of analysis, Far-r-is’smethod, Sattath and Tversky’s method, Li’s method, and Tateno etal.‘s modified Fan-is method. The new, neighbor-joining method and Sattath andTversky’s method are shown to be generally better than the other methods.

Introduction
In the construction of phylogenetic trees, the principle of minimum evolution
or maximum parsimony is often used. The standard algorithm of the tree-making
methods based on this principle is to examine all possible topologies (branching patterns) or a certain number of topologies that are likely to be close to the true tree andto choose one that shows the smallest amount of total evolutionary change as the final
tree. This method is quite time consuming, and, when the number of operationaltaxonomic units (OTUs) is large, only a small proportion of all possible topologies is
examined. However, there are methods in which the process of searching for theminimum evolution tree is built into the algorithm, so that a unique final topology is
obtained automatically. Some examples are the distance Wagner (DW) method (Farris1972), modified Farris (MF) methods (Tateno et al. 1982; Faith 1985), and the neighborliness methods of Sattath and Tversky (ST method; 1977) and Fitch ( 198 1). Thesemethods are not guaranteed to produce the minimum-evolution tree, but their efficiency in obtaining the correct tree is often better than that of the standard maximumparsimony algorithm (Saitou and Nei 1986). In the following we would like to present
a new method (the neighbor-joining [NJ] method) that produces a unique final tree
under the principle of minimum evolution. This method also does not necessarily
produce the minimum-evolution
tree, but computer simulations have shown that it
I. Key words: phylogenetic tree, neighbor-joining method, minimum-evolution

tree, parsimonioustree.
2. Current address: Department of Anthropology, University of Tokyo, Tokyo 113, Japan.
Address for corespondence and reprints: Dr. Masatoshi Nei, Center for Demographic and Population
Genetics, Graduate School of Biomedical Sciences, The University of Texas Health Science Center at Houston,
P.O. Box 20334, Houston, Texas 77225.
Mol. Biol. Evol. 4(4):406-425. 1987.
0 1987 by The University of Chicago. All rights reserved.
0737-4038/87/0404-0007$02.00

406

Neighbor-joining Method407

is quite efficient in obtaining the correct tree topology. It is applicable to any type of
evolutionary distance data.
Algorithm

The algorithm of the NJ method is similar to that of the ST method, whose
objective is to construct the topology of a tree. Unlike this method, however, the NJ
method provides not only the topology but also the branch lengths of the final tree.
Before discussing the algorithm of the present method, let us first define the term
“neighbors.” A pair of neighbors is a pair of OTUs connected through a single interior
node in an unrooted, bifurcating tree. Thus, OTUs 1 and 2 in figure 1 are a pair ofneighbors because they are connected through one interior node, A. There are twoother pairs of neighbors in this tree (viz., [5, 61 and [7, 81). The number of pairs ofneighbors in a tree depends on the tree topology. For a tree with N (24) OTUs, theminimum number is always two, whereas the maximum number is N/2 when N isan even number and (N - 1)/2 when N is an odd number.
If we combine OTUs 1 and 2 in figure 1, this combined OTU ( l-2) and OTU 3
become a new pair of neighbors. It is possible to define the topology of a tree by
successively joining pairs of neighbors and producing new pairs of neighbors. For
example, the topology of the tree in figure 1 can be described by the following pairsof neighbors: [l, 21, [5, 61, 17, 81, [l-2, 31, and [l-2-3, 41. Note that there is anotherpair of neighbors, [5-6,7-81, that is complementary to [l-2-3,4] in defining the topology.
In general, N - 2 pairs of neighbors can be produced from a bifurcating tree of N
OTUs. By finding these pairs of neighbors successively, we can obtain the tree topology.
Our method of constructing a tree starts with a starlike tree, as given in figure2(a), which is produced under the assumption that there is no clustering of OTUs. In

2
6

FIG. 1.-An
lengths.

unrooted tree of eight OTUs, l-8. A-F are interior nodes, and italic numbers are branch

408

Saitou and Nei

practice, some pairs of OTUs are more closely related to each other than other pairs
are. Consider a tree that is of the form given in figure 2(b). In this tree there is only
one interior branch, XY, which connects the paired OTUs (1 and 2) and the others(3, 4, . . . , N) that are connected by a single node, Y. Any pair of OTUs can take thepositions of 1 and 2 in the tree, and there are N(N - 1)/2 ways of choosing them.
Among these possible pairs of OTUs, we choose the one that gives the smallest sumof branch lengths. This pair of OTUs is then regarded as a single OTU, and the next
pair of OTUs that gives the smallest sum of branch lengths is again chosen. This
procedure is continued until all N - 3 interior branches are found.
The sum of the branch lengths is computed as follows: Let us define Do and Lab
as the distance between OTUs i and j and the branch length between nodes a and b,
respectively. The sum of the branch lengths for the tree of figure 2(a) is then given by(1)
since each branch is counted N - 1 times when all distances are added. On the other
hand, the branch length between nodes X and Y (Lxr) in the tree of figure 2(b) is
given by1=-[
Lxy 2(N-2)

5 (Dlk+D2k)-(N_2)(L1x+L2~)-2
k=3

5

&I-

(2)

i=3

The first term within the brackets of equation (2) is the sum of all distances that
include Lxy, and the other two terms are for excluding irrelevant branch lengths. If
we eliminate the interior branch (XY) from figure 2(b), two starlike topologies (onefor OTUs 1 and 2 and the other for the remaining N - 2 OTUs) appear. Thus, Llx+ L2x and c: 3 Liy can be obtained by applying equation (1):
Lx+L2x=D12,

5LiY=&

i=3

8

FIG. 2.-(a),clustered.

W

2 4.

(3b)

3si-cj8

A starlike tree with no hierarchical structure; and (b), a tree in which OTUs 1 and 2 areNeighbor-joining Method

409

Adding these branch lengths, we find that the sum (SIz) of all branch lengths of thetree in figure 2(b) becomes

It can be shown that equation (4) is the sum of the least-squares estimates of branch
lengths (see Appendix A).
In general, we do not know which pairs of OTUs are true neighbors. Therefore,
the sum of branch lengths (S,) is computed for all pairs of OTUs, and the pair thatshows the smallest value of Sii is chosen (inferred) as a pair of neighbors. In practice,even this pair may not be a pair of true neighbors; but, for a purely additive tree withno backward and parallel substitutions, this method is known to choose pairs of trueneighbors (see the following section- Criterion for Minimum-Evolution
Tree-fordetail). At any rate, if S12 is found to be smallest among all Sij values, OTUs 1 and 2
are designated as a pair of neighbors, and these are joined to make a combined OTU(l-2). The distance between this combined OTU and another OTU j is given by
Dc,-2)j=(Dlj(3 5jIN).

+ D2jU(5)

Thus, the number of OTUs is reduced by one, and, for the new distance matrix, theabove procedure is again applied to find the next pair of neighbors. This cycle isrepeated until the number of OTUs becomes three, where there is only one unrooted tree.
The branch lengths of a tree can be estimated by using Fitch and Margoliash’s(1967) method. Suppose that OTUs 1 and 2 are the first pair to be joined in the treeof figure 1. Llx and L2X are then estimated byLx= (012

+

Dlz-

DzzW,

(64

where D,z = (Cy= 3 Dli)/(N - 2) and D 2z = (2;’ 3 Dzi)/(N - 2). Here, 2 representsa group of OTUs including all but 1 and 2, and Dlz and Dzz are the distances between
1 and 2 and 2 and 2, respectively (see Nei 1987, pp. 298-302, for an elementary
exposition of this method). LIX and L2X are the least-squares estimates for the tree of
figure 2(b) (see Appendix A), and they are estimates of LIA and Lu, respectively, infigure 1. Once L lA and Lu are estimated, OTUs 1 and 2 are combined as a single
OTU (l-2), and the next neighbors are searched for. Suppose that (l-2) and 3 are thenext neighbors to be joined, as in figure 1. Branch lengths Lt1_2jBand L3* are obtained by applying equations (6a) and (6b). Furthermore, LAB is estimated by
L, 1_2)B- (D12)/2. The above procedure is applied repeatedly until all branch lengths
are estimated. If a tree is purely additive, this method gives the correct branch lengths
for all branches (see Appendix B).
The principle of the NJ method can be extended to character-state data such as
nucleotide or amino acid differences. In this case, one can use the total number of4 10

Saitou and Nei

Table 1Distance Matrix for the Tree in Figure 1
OTU
OTU

1

2

3

4

5

6

7

.
.
.
5 .

7
8
11
13
16
13
17

5
8
10
13
10
14

5
7
10
7
11

8
11
8
12

5
6
10

9
13

8

2
3
4
6
7
8

.

substitutions in place of the sum of branch lengths ($), though the actual procedure
is a little more complicated than that given above (Saitou 1986, pp. 90-98). However,since the algorithm turns out to be very similar to that of Hartigan (1973), we shallnot present it here. Note also that most character-state data can be converted into
distance data so that the above simpler algorithm applies.
An example: consider the distance matrix given in table 1. The distance D, inthis matrix is obtained by adding all relevant branch lengths between OTUs i and j
in figure 1 under the assumption that there is no backward and parallel substitution.
The result of application of the NJ method is presented in table 2 and figure 3. In the
Table 2SC Matrices for Two Cycles of the NJ Method for the Data in Table 1= [ 1, 21

A. Cycle 1: Neighbors

OTU
OTU
2
3
4
5
6
7
8

.
.
.
.

1

2

3

4

5

6

7

36.67
38.33
39.00
40.33
40.33
40.17
40.17

38.33
39.00
40.33
40.33
40.17
40.17

38.67
40.00
40.00
39.83
39.83

39.67
39.67
39.50
39.50

37.00
38.83
38.83

38.83
38.83

37.67

B. Cycle 2: Neighbors= [5, 61

OTU
OTU

1-2

3

4

5

6

7

.

31.50
32.30
33.90
33.90
33.70
33.70

32.30
33.90
33.90
33.70
33.70

33.70
33.70
33.50
33.50

31.30
33.10
33.10

33.10
33.10

31.90

3
4
5
6
7
8

.
.
.

Neighbor-joining Method

4 11

6
5

(b)

(a>

(c)

1
2

8

4

5

Cd)

(e)

(f)

FIG. 3.-Application of the neighbor-joining method to the distance matrix of table 1. Italic numbers
branch lengths, and branches with thicker lines indicate that their lengths have been determined.

search for the first pair of neighbors (cycle l), OTUs 1 and 2 are chosen because Si2
is smallest among 28 &‘s (see table 2). Si2 (= 36.67) is smaller than the sum (SO= 39.28) of branch lengths of the starting starlike topology, but, interestingly, some&‘s are larger than S o. Diz and Dzz in equations (6a) and (6b) become 13 and 10,
respectively. Thus, the branch lengths LIA and Lu are obtained to be (7 + 13 - lo)/
2 = 5 and (7 + 10 - 13)/2 = 2, respectively, which are identical with those of thetrue tree in figure 1 (fig. 3[a]). OTUs 1 and 2 are then combined, and the average
distances (Dt1_2Jj;
j = 3, . . . , 8) are computed by equation (5). In the next step (cycle
2 in table 2), OTUs 5 and 6 are found to be a pair of neighbors, and L5E and LeE areestimated to be 1 and 4, respectively, which are again identical with those of the true
tree (fig. 3[b]). In cycle 3, OTUs (l-2) and 3 are chosen as a pair of neighbors, andthe branch lengths for LsBand LcI_2)B
become 1 and 5.5, respectively. Thus, the branch
length LABis estimated to be 5.5 - y2 = 2. These are again the correct values (fig.
3[c]). In cycle 4, [l-2-3, 41 is identified as a pair of neighbors (fig. 3[d]), and in cycle
5 [l-2-3-4, 5-61 is chosen. The choice of the latter pair of neighbors automaticallyleads to the identification of the final pair of neighbors [7, 81. The Sij for [l-2-3-4, 56] is identical with that for [7, 81. The topology of the reconstructed tree is therefore
given by figure 3(e), which is identical with that of figure 1. The branch lengths LTF(=LTz = 2) and LgF(=LgZ = 6) are obtained by using equations (6a) and (6b), whereas- Do_2_3_4M5_6j/2= 2 (fig. 3[f]). It is thus clear that all branch
LDFbecomes L, 1_2_3_4_5_6Jz
lengths as well as the topology are correctly reconstructed in the present case.
Criterion for the Minimum-Evolution TreeIn this section, we first show that the algorithm developed above produces thecorrect tree for a purely additive tree. We shall then discuss a criterion for the minimumevolution tree.

412

Saitou and Nei

Consider a tree for N (24) OTUs and assume that OTUs 1 and 2 are a pair oftrue neighbors. For an additive tree, we obviously have the following inequalities.
D,,+Do<Dli+D2j,

012 + D, < D,j + Dzi,

(7)

where i and j are any OTUs (3 I i < j I N). Under this condition, it can be shownthat S12 is smaller than Slj or ,S’zj(3 I j I N). To show this, let us consider the pairingof OTUs 1 and 3 as an example. The total length of the tree with this pairing can be
written as-&;

Uhk+D3k,+&+3.

k-2

k23

In a similar manner,

-&

i

UW

(D,*+D2k)+&D,2.

k-3

Hence,$3~&2=&+,3-D,,)+&$@2k-D3k)

k-4

=&

t

(9)

[(D13+Dzk)-(D12+D3&

k-4

If we note the inequalities in formula (7), Di2 + D3k < Di3 + D2k (4 I k I N).
Therefore, S13 > S 12. The same inequality also holds for any other pairs involving
OTUs 1 and 2: Sy > Si2 and S2j > Si2 (3 I j I N). Furthermore, in our algorithm
we search for a pair of OTUs that shows the smallest Sti. Therefore, if OTUs 1 and 2
are such a pair, Si2 must be smallest among all $‘s. However, this is not what weneed in our algorithm. Our algorithm requires that if S12 is smallest among all &‘s,
OTUs 1 and 2 are neighbors. Proof of this theorem is somewhat complicated, but it
can be done (see Appendix C). Therefore, our algorithm produces the correct unrooted
tree for a purely additive tree.
Of course, actual data usually involve backward and parallel substitutions, sothat there is no guarantee that the correct topology is obtained by the NJ method.
However, computer simulations, which will be discussed below, have shown that,compared with other methods, the NJ method is efficient in obtaining the correct
topology.
In constructing the topology of a tree, Sattath and Tversky (1977) and Fitch(198 1) used the inequalities in formula (7). Their method is to count the number(neighborliness) of cases satisfying formula (7) for each pair of OTUs and choose thepair showing the largest number as neighbors. Since Sattath and Tversky’s (1977)
algorithm uses equation (5) for making the new distance matrix, their method is
expected to give a result similar to ours. Fitch ( 198 1) uses interior-distance matrices
for constructing the topology, so that his algorithm is different from ours. Nevertheless,these three methods as well as some other tree-making methods require the same4 14

Saitou and Nei

An example: applying his neighborlinessmethod to Case’s (1978) data on immunologicaldistance, Fitch ( 198 1) constructed a phylogenetic tree of nine frog (Rana)
species. If we use the NJ method, a slightly different tree is obtained (fig. 5); that is,while the closest species to the R. aurora-R. boyhi group is R. cascadae in Fitch’s tree,
it is R. muscosa in our tree. The latter topology is also obtained by the ST method.
We can apply the minimality
test in formula ( 10) to see which topology is morereasonable. The test can be done if we consider the four OTU groups, i.e., the aurora
and boylii group, muscosa, cascadae, and the remaining
five species. Applicationofthe test supports the topology presented in figure 5 rather than Fitch’s. Comparisonof the sum of branch lengths between the two topologies also supports the topology
in figure 5. (This particular comparison
was conducted under the condition that all
branch lengths are nonnegative
and that each estimated [patristic] distance is greaterthan or equal to the corresponding
observed distance, because Fitch’s tree was constructed under this condition.)
We also note that the branch lengths estimated by the
NJ method are close to those estimatedby a linear programming
method (seeFitch 198 1).
Efficiency of the NJ Method in Recovering

the Correct Topology

Since the exact evolutionary
pathways of extant organisms are usually unknown,
it is not suitable to use real data for examining the efficiency of a tree-making
method.
Therefore, we employed a computer simulation,comparing reconstructed
trees with
their model trees. In this study we compared the efficiency of the NJ method withthat of five other methods: UPGMA (Sokal and Sneath 1963), the DW method, theST method, Li’s (LI; 198 1) method, and the MF method. The LI method is a transformed distance method (see Nei [ 1987, pp. 302-3051 for the explanationof the transformed distance method), and the MF method is a modificationof Farris’s (1972)
method. All these methods produce a unique parsimonious
tree from distance data.
We considered both cases of constant and varying (expected) rates of nucleotide substitution.
R. aurora
R. boylii
muscosacascadaetemporaria
24.0

pretiosa
R. catesbeiana
R. pipiens
R. tarahwnarae

FIG. 5.-Treeobtainedby the NJ methodfrom immunologicaldistance

data of Case (1978)

Neighbor-joining Method

4 15Constant Rate of Nucleotide Substitution
To examine the effect of topological differences, we considered two different
model trees (trees [A] and [B] of fig. 6), both of which consist of eight OTUs. Modeltree (A) has two neighboring pairs ([ 1, 21 and [7, S]), whereas model tree (B) has four([l, 21,[3,41, 15,619 and [7, 81). To make the effect of branch lengths comparable forthe two model trees, we assumed that the interior branch length (a) is the same forboth trees. We also tried to make the average @) of all pairwise distances (0;s) nearly
the same for the two trees. Hence, we set c = b + 3a or c = b + 3a, where a, b, andc are the expected branch lengths (expected numbers of nucleotide substitutions persite) given in figure 6. In a computer simulation conducted with the same topology
as that of model tree (A), Tateno et al. (1982) set a = b. In the present study, we set
a Q b in model tree (A) so that the differences between different Du’s were relatively
smaller. This makes it more difficult to reconstruct the correct tree than in the caseof Tateno et al.‘s simulation.
The scheme of the computer simulation used is as follows: The ancestral sequenceof a given number of nucleotides was generated by using pseudorandom numbers,
and this sequence was assumed to evolve according to the predetermined branchingpattern of the model tree. Random nucleotide substitutions were introduced in eachbranch of the tree following a Poisson distribution with the mean equal to the expectedbranch length. Although the expected rate of nucleotide substitution was the same forall lineages, the actual number of substitutions varied considerably with lineage becauseof stochastic errors. After the nucleotide sequences for eight OTUs were produced,
nucleotide differences were counted for all pairs of sequences, and the evolutionary
distance (Jukes and Cantor 1969) was computed for each pair of OTUs. With the six
tree-making methods mentioned above, tree topologies were determined from data
either on the proportion of different nucleotides between the two sequences compared(p) or on the Jukes-Cantor distance (d). Note that p is a metric, whereas d is not. The
entire process of simulation was repeated 100 times.
Two measures are used to quantify the efficiency of a tree-making method inrecovering the topology of the model tree. One is the proportion (PC) of correct trees(topologies) obtained. The other is the average distortion index (Tateno et al. 1982)
based on Robinson and Foulds’ (198 1) metric on tree comparison. The distortion
index (&) is twice the number of branch interchanges required for a reconstructed
tree to be converted to the true tree. Here, we consider only unrooted trees.

C

a

C

0.5a-I

C

a

C
C
C

a
Q

i___

6a+b

i

L
0.502

C
P

b

I

I

2
3
4
5
6
7
8

(B)
FIG. 6.-Modeltrees (A) and (B) under the assumption of constant rate of nucleotide substitution

4 16

Saitou and Nei

Table 3P, and dT (in parentheses) for Six Tree-making Methods for
the Case of a = 0.01, b = 0.04, and c = 0.07
MODEL TREEBa

MODEL TREEAa
METHOD

300

600

900

300

600

900

UPGMA:
pb . . .
..
d” . . . . . . .

14 (3.18)
15 (3.18)

36 (1.72)
34 (1.74)

58 (0.98)
56 (1.04)

14 (4.54)
13 (4.56)

36 (2.74)
35 (2.70)

51 (1.68)
52 (1.60)

39 (1.76)
38 (1.92)

73 (0.58)
72 (0.62)

95 (0.10)
95 (0.10)

24 (2.86)
19 (2.94)

51 (1.30)
48 (1.42)

67 (0.76)
64 (0.86)

DW:
p .
... .
d .... ...

42 (1.70)
37 (1.74)

75 (0.54)
74 (0.58)

96 (0.08)
95 (0.10)

26 (2.36)
28 (2.36)

55 (1.12)
58 (1.06)

79 (0.48)
79 (0.46)

LI:
p ....
d ... ....

41 (1.58)
36 (1.84)

71 (0.70)
66 (0.82)

94 (0.12)
89 (0.24)

40 (2.04)
39 (2.10)70 (0.78)
70 (0.78)

90 (0.22)
90 (0.26)

ST:
p ... ....
d ........

48 (1.26)
44 (1.48)

75 (0.54)
70 (0.62)

97 (0.06)
96 (0.08)

45 (1.66)
43 (1.62)75 (0.62)
74 (0.64)

91 (0.22)
91 (0.22)

NJ:
p ..,..
d ..
..,48 (1.36)
41 (1.60)

76 (0.54)
70 (0.62)

97 (0.06)
96 (0.08)

46 (1.64)
45 (1.62)

76 (0.60)
75 (0.60)

91 (0.20)
91 (0.20)

MF:
p ........
d.....

‘As shown in fig. 6.
b Trees reconstructed from data on the proportion of different nucleotides between the sequences compared.
’Trees reconstructed from the Jukes-Cantor distance.

Table 3 shows the results for the case of a = 0.01, b = 0.04, and c = 0.07, where
the B for all OTUs is 0.16 for both model trees. It is clear that in all tree-making
methods PC increases as the number of nucleotides used (n) increases, whereas dT
decreases. This is of course due to the fact that the sampling error of the distance
between a pair of OTUs decreases as n increases. The PC and dT values obtained byusing p and d are nearly the same, though p tends to show a better performance inrecovering the correct topology, particularly for model tree (A).
In the case of model tree (A) UPGMA shows the poorest performance in termsof both criterion PC and criterion dr. Even when 900 nucleotides are used, the proportion of correct trees obtained is -57%. The other five tree-making methods show
a much better performance than UPGMA, and when 900 nucleotides are used, PCis
-95%. Interestingly, all of them show a similar performance for all n’s examined. In
the case of model tree (B), UPGMA again shows a poorer performance than any othermethod. In this case, however, all the five methods do not show the same performance.
Rather, the NJ and the ST methods are better than the LI method, which is in turnbetter than the DW and MF methods.
The results for the case of a = 0.02, b = 0.13, and c = 0.19 are presented in table
4. The D for this case is 0.42 for model tree (A) and 0.43 for model tree (B). For
model tree (A), UPGMA shows an improved performance compared with the case intable 3. However, all other methods show a small value of PCand a larger value of dT
than those in table 3. This is apparently due to the fact that there are more backward

Neighbor-joining

Method4 17

Table 4
P, and dT (in parentheses)

for Six Tree-making

Methods

forthe Case of a = 0.02, b = 0.13, and c = 0.19
MODEL TREE Aa

MODEL TREE Ba

METHOD

300

600

900

300

600

900

UPGMA:p...
d...

15 (3.24)
15 (3.28)

50 (1.32)
49 (1.34)

62 (0.82)
61 (0.84)

11 (4.62)
13 (4.50)

28 (2.94)
30 (2.90)

54 (1.48)
57 (1.44)

MF:
p . .....
d .. ....

34 (2.38)
30 (2.70)

65 (0.82)
62 (1.02)

79 (0.44)
76 (0.54)

10 (4.00)
9 (4.12)25 (2.22)
22 (2.28)

43 (1.48)
43 (1.48)

DW:
p . ....
d .......

27 (2.40)
27 (2.52)

66 (0.96)
62 (1.02)

77 (0.54)
70 (0.70)

17 (3.54)
18 (3.54)

39 (1.92)
36 (1.98)

54 (1.10)
53 (1.16)

23 (2.60)
20 (2.82)

44 (1.34)
33 (1.78)

67 (0.80)
55 (1.12)

25 (3.54)
20 (3.70)

50 (1.52)
49 (1.54)

81 (0.52)
81 (0.50)

35 (2.06)
26 (2.42)

67 (0.74)
61 (0.96)

82 (0.38)
78 (0.48)

34 (2.40)
31 (2.50)

60 (1.08)
58 (1.16)

82 (0.38)
83 (0.36)

36 (2.14)
26 (2.38)

64 (0.88)
58 (1.08)

83 (0.34)
78 (0.48)

34 (2.32)
33 (2.56)

63 (0.96)
61 (1.04)

82 (0.36)
83 (0.34)

LI:
P

p
d

..

..
....

NJ:
p .. ..
d . ...
NOTE.-Notations

are as in table 3.

a As shownin fig. 6.

and parallel substitutions involved in this case. Nevertheless, UPGMA still shows apoorer performance than all other methods except LI, which is less efficient than
UPGMA for the case of n = 600. The NJ, ST, DW, and MF methods give similar
results, though the first two methods give slightly better results than the others for n= 900. We also note that p gives a better result than d for all methods but UPGMA,
for which both p and d give essentially the same results. In the case of model tree (B),the PC values for UPGMA are not necessarily higher than those in table 3, but they
are higher than those for the MF method for the same case. The DW method also
shows a rather poor performance, though it is slightly better than the UPGMA andMF methods. The NJ and ST methods again show the best performance, but their PC
values are slightly lower than those for the case of table 3. The LI method is quite
good but not as good as the NJ and ST methods. Interestingly, p and d give similar
results for all methods, unlike the case of model tree (A).
Table 5 shows the results for the case of a = 0.03, b = 0.34, c = 0.42, and D= 0.92 for tree (A) and 0.91 for tree (B). Compared with the two previous cases, thefrequency of backward and parallel substitutions is expected to be much higher becauseof the larger D, values used. Therefore, we used n = 500, 1,000, and 2,000 for thiscase. Yet, the PCvalues are smaller than those for the two previous cases. The relative
merits of different tree-making methods for the case of model tree (A) are more orless the same as those for the case of table 4, except that the LI method tends to show
a poorer performance than UPGMA. When n = 500, the MF and DW methods show
a slightly higher value of PC than the ST and NJ methods, but for the other two n4 18

Saitou and Nei

Table 5
PC and dT (in parentheses) for Six Tree-making Methods for
the Case of Q = 0.03, b = 0.34, and c = 0.42
MODEL TREE Aa

MODEL TREE Ba

500

1,000

2,000

500

1,000

2,000

9 (3.78)
9 (3.78)

27 (2.10)
27 (2.10)

62 (0.86)
62 (0.88)

10 (5.20)
11 (5.30)

18 (3.76)
18 (3.74)

54 (1.32)
55 (1.26)

P

15 (4.02)
13 (4.42)

41 (1.82)
34 (2.14)

62 (0.92)
55 (1.14)

3 (5.68)
3 (5.72)

17 (3.64)
13 (3.80)

28 (2.40)
26 (2.48)

p......
d ... .. .

16 (3.78)
15 (4.22)

46 (1.54)
40 (1.96)

63 (0.82)
58 (0.98)

4 (5.42)5 (5.50)

18 (3.28)
18 (3.48)

41 (1.72)
35 (1.82)

p ......
d3 (4.26)
3 (4.84)

37 (2.00)
25 (2.60)

53 (1.18)
39 (1.66)

15 (4.48)
12 (4.72)

28 (2.98)
27 (3.06)

70 (0.90)
66 (1.02)

P

10 (3.56)
6 (4.06)44 (1.62)
40 (1.82)

68 (0.76)
56 (1.04)

13 (4.00)
10 (4.32)

36 (2.34)
34 (2.34)

74 (0.62)
71 (0.72)

p .
....
d ....

11 (3.70)

44 (1.68)
38 (2.00)

67 (0.80)
57 (1.06)

13 (4.46)
14 (4.44)

34 (2.38)
32 (2.42)75 (0.62)
73 (0.72)

METHODUPGMA:P ...
d
MF:“‘‘-”.

LI:

NOTE.--Notations

5 (4.24)
are as in table 3.

a As shown in fig. 6.

values they show more or less the same performance. Data on p again give a better
result for the five methods (except for UPGMA) than do those on d. In the case ofmodel tree (B), the MF method shows a poorer performance than UPGMA, which
now gives results similar to the DW method. However, the P, values for the LI, ST,
and NJ methods are substantially higher than those for UPGMA and the DW methods.
Although the above computer simulations were done for a limited number ofcases, the results obtained may be summarized as follows: (1) The efficiency of the NJ
method in recovering the true unrooted tree is virtually the same as that of the ST
method. (2) The NJ and ST methods perform well for both model tree (A) and model
tree (B), whereas the DW and MF methods are good only for tree (A) and the LI
method is good only for tree (B). For both model trees, UPGMA is rather poor inrecovering the true unrooted tree. (3) In the case of model tree (A), data on p tend to
give slightly better results than those on d, except for UPGMA. For model tree (B),however, both p and d give similar results.
Conclusion (3) above indicates that data on p are better than those on d forconstructing a topology, particularly when the OTUs used form a topology similar tomodel tree (A). However, since p is not a linear function of nucleotide substitutions,
it does not provide good estimates of branch lengths unless the p values are very small.
It is therefore advised that once a topology is obtained by using data on p, branch
lengths should be estimated by using data on d.
Tateno et al. ( 1982) and Sourdis and Krimbas ( 1987) conducted similar computersimulation studies, comparing the efficiency of the UPGMA and the DW and MF
methods as well as Fitch and Margoliash’s ( 1967) method for model tree (A). AlthoughNeighbor-joining Method

4 19

the parameter values used in their simulations are different from ours, their conclusions
with respect to unrooted trees are more or less the same as ours.
Varying Rate of Nucleotide Substitution
When the rate of nucleotide substitution varies from evolutionary lineage toevolutionary lineage, the probability of obtaining the correct tree is expected to belower than that for the case of rate constancy. To see the effect of this factor on PC,
we conducted another computer simulation.
In this simulation, we used the two model trees ([A] and [B]) given in figure 7.
The topologies of trees (A) and (B) in fig. 7 are identical, respectively, with those oftrees (A) and (B) in figure 6. The value given for each branch of these trees is the
expected branch length (the expected number of nucleotide substitutions per site). The
expected branch lengths for tree (A) in figure 7 were obtained under the assumptionthat b in figure 6(A) varies according to the gamma distribution with mean 0.04 andvariance 0.08 (see Tateno et al. 1982 for the justification of this procedure). Similarly,
the expected branch lengths for tree (B) in figure 7 were obtained under the assumptionthat c in figure 6(B) varies according to the gamma distribution with mean 0.07 andvariance 0.14. The value of a and the expectation of D over all branches were 0.01
and 0.016, respectively. Therefore, the simulations for model trees (A) and (B) correspond, respectively, to those for trees (A) and (B) in table 3. Once the expected
length of a particular branch was determined, the actual number of nucleotide substitutions for that branch was obtained by using the Poisson distribution. The eight
nucleotide sequences thus obtained were used for the construction of phylogenetic
trees. This process was repeated 100 times. In this simulation, only the case of 600
nucleotides was examined, and the trees were constructed by using the p values.
The results of this simulation are presented in table 6. One striking feature inthis simulation is that the performance of UPGMA was very poor and that in noneof the 100 replications was the correct tree obtained for both model tree (A) and model
tree (B). This is in sharp contrast to the case of rate constancy (table 3), in which the
PCfor UPGMA is 36% when n = 600. The effect of varying rate on the PCvalue is lessnoticeable for the other tree-making methods. The PC values for the LI method aresomewhat lower than those for the case of constant rate (see tables 3 and 6). In theremaining four methods, the PCvalues are virtually the same for both cases of constant

I.04

L
.ol

.ol

.ol

07

.05

.16

.09

FIG. 7.-ModelEU

trees (A) and (B) under the assumption of varying rate of nucleotide substitution

420

Saitou and Nei

Table 6
P, and dT (in parentheses) for Six Tree-making Methods for
the Case of Varying Rate of Nucleotide Substitution
MethodUPGMA: p
MF:p . . . .
DW: p . . .
LI: p . . . .
ST:p
....
NJ:p
....

Model Tree Aa

Model Tree B”

0 (8.06)
77 (0.50)
69 (0.72)
46 (1.30)
77 (0.50)
75 (0.56)

0 (9.74)
57 (1.46)
59 (1.26)
45 (1.68)
69 (0.82)
72 (0.78)

Nom.--Notations
are as in table 3.
‘As shown in fig. I.

and varying rates of nucleotide substitution. Therefore, the conclusions obtained forthe case of constant rate also apply to the case of varying rate as far as the NJ, ST,
MF, and DW methods are concerned.
Discussion

Y

Unlike the standard algorithm for minimum-evolution
trees, the NJ methodminimizes the sum of branch lengths at each stage of clustering of OTUs starting witha starlike tree. Therefore, the final tree produced may not be the minimum-evolution
tree among all possible trees. However, it should be noted that the real minimumevolution tree is not necessarily the true tree. Saitou and Nei (1986) have shown thatthe minimum-evolution or maximum-parsimony tree often has an erroneous topology
and that the maximum-parsimony method of tree making is not always the best in
recovering the true topology. It seems to us that the relative efficiencies of different
tree-making methods should eventually be evaluated by computer simulation. Our
computer simulation has shown that the NJ method is quite efficient compared withother tree-making methods that produce a single parsimonious tree.
We have shown that the estimates of branch lengths of the tree obtained by the
NJ method are least-squares estimates determined at each stage of clustering of OTUs.
This does not mean that these estimates are identical with those that are obtainable
by the least-squares method for all branches of the final tree topology. Nevertheless,this property gives some assurance about the reliability of the estimates of branch
lengths. Particularly when the number of OTUs is four or less, the branch lengths areexactly least-squares estimates, as is clear from equation (A4) below.
Our procedure of estimating branch lengths is essentially the same as that of
Fitch and Margoliash ( 1967). Some estimates of branch lengths may therefore become
negative. If one is reluctant to accept negative estimates, there are two ways to eliminate
them. One is to impose the condition that all branches be positive and then to reestimate
the branch lengths. The other is to assume that negative estimates are due to sampling
error and that the real values are zero rather than negative. Under this assumption,
one may simply convert all negative estimates to zero. The second method is justified
if we note that the absolute values of negative estimates are usually very small.
A computer program for constructing a tree by using the NJ method is available
from the authors on request.

Neighbor-joining

Method42 1

Acknowledgments

We thank Clay Stephens for his comments and John Sourdis for his help incomputer simulation. This study was supported by research grants from the National
Institutes of Health and the National Science Foundation.
APPENDIXA

Least-Squares Estimation of the Branch Lengths

Let us consider the tree of figure 2(b). If we use matrix notation, the problem isto obtain the least-squares solution of the linear equation Ax = d, where x is a column
vector of N + 1 branch lengths (x’ = [LIX, LzX, LsY, L4r, . . . , LNY, Lxr]), d is acolumn vector of N(N - 1)/2 pairwise distances (d’ = [I&, I&, Q4, D15, . . . DIN,
IJN]), and A is an [N(N - 1)/2] X (N + 1) matrix. The. . . D,,_
023,
024,
- - - D2iv,
element of the ith row and the jth column of matrix A is given by1
au=
0

if the ith distance includes the jth branch.

[ otherwise

An example of A for N = 5 is shown below:-110000
101001100101100011011001010101
010011001100001010
-0
0 0

A=

1

10

,1

The least-squares solution of the equation Ax = d is given by solving the equation
A’Ax = A’d. It becomes x L = B-‘A’d, where B = A’A. The general expressions ofsymmetric matrices B and B-’ are
N-l
1

1
N-l

1
1

!

...

1
1

N-2
N-2

....

1.

.

2.

N-l
2

2
2(N - 2)

. . .

.
i

1.

.
i

N-l.
i

.:.

N-2

N-2

2

..e

a
b
0
B-l=

.

b
a
0

0
0
c

.

0
0
d

0
0
d

-0
*.
**
l

l

. . . . .
. . . . .
;, ;, >;;I.:.;je e f f f l_

.

0
0
d

e
e
f

(Al)

WV

f g

422

Saitou and Nei

where a = N/4(N - 2), b = (N - 4)/4(N - 2), c = (2N2 - 11N + 16)/2(N - 2)2(N- 3), d = -(N - 4)/2(N - 2)2(N - 3), e = -V&f= -1/2(N - 2)(N - 3), and g = (N- 1)/4(N - 3). Therefore, xL becomes1

Llx=-D,2+21

L2x=-D,2+
2

Li,=-

N!ZU’

-‘(P+
’ (N-2)2

j&P-Q>,

6434

g&Q-PI,

WV

N-4

Q)-

(N-

2)2(N-(3 5 is N)

3) ’

(A3d)

i(P+Q)-fD”-(N_2;(N_3)K
Lxy= 2(N-2)

where P = 25 3 D,j, Q = Zc 3 D,, Ui = 2;;(A3c)D, (i 2 3), and V = C3 sj < k Djks Notethat equations (A3a) and (A3b) are equivalent to equations (6a) and (6b), respectively.
Thus, the sum of branch lengths (Si2) for the topology in which OTUs 1 and 2 areclustered becomes
S12 = LlX+ L2X+ 5 Liy+ Lxy=&(p+~)++h+&v.

bw

i=3

Equation (A4) is equivalent to equation (4).
APPENDIX B

Branch Lengths for a Purely Additive Tree

Let us consider the tree given in figure 1. If the tree is purely additive, 012 = LIA
and Dlj - D, = LIA - La (3 I j I N). Substituting these equations into
equation (6a), we have+ Lu

L~x=;(LIA+L4)+
$-+

W-

WL4

-

= LA.

Ldl

(A%

The estimated branch length (Llx) is identical with the true one (LIA). The same thingcan be proven for L 2x. Therefore, the node X is identical with the node A in the tree
in figure 1.
If OTUs 1 and 2 are neighbors, they are combined into a single OTU, ( l-2).
Suppose that OTUs (l-2) and 3 are a new pair of neighbors. The estimates of branch
lengths for AB and 3B can then be obtained correctly, as shown below. Since the tree
is purely additive,
D(,-2)3=(03+023)/z

=

[(&A

+

~543)

+

&-I

+

L/13)1/2= 012/z

+

P13

-

LIA)

(A64

and
D (1-W -D3j=D,2/2+(DU-LIA)-(L3B+LBj)=D12/2+LAB-L3B(j24).

(A6b)

Neighbor-joining Method

423

Substituting these into equation (6a), we have1[(N-3)(D,2/2+L~~-L3~)1+ 2(N- 3)=;D~z+h.

Since LAX = Lt1_2jx- D&2, L AX= LAB. On the other hand, as before, it easily can beshown that L 3x = LjB. Therefore, X = B.
The above argument can be applied to any situation if the additivity of branch
lengths is maintained.
APPENDIX C

The Smallest Sii Gives the True Neighbors

In the following, we show that for a purely additive tree OTUs 1 and 2 are trueneighbors when Si2 is smallest among all &‘s. We first show this for the case of four
OTUs and then use the principle of induction to prove that it is generally true.
Using the results presented in the Criterion for the Minimum-Evolution
Treesection, we can state that the condition for S12 to be smallest among the six Sij’s forfour OTUs isand

Our task is to show that if Si2 is smallest, OTUs 1 and 2 are true neighbors. In thecase of four OTUs, OTUs 3 and 4 are also neighbors if OTUs 1 and 2 are neighbors(see fig. 4). We prove our assertion by showing that when S12 is smallest, only OTUs
1 and 2 (and OTUs 3 and 4) are neighbors. To prove this, we first assume that OTUs
1 and 3 (and 2 and 4) are neighbors. We then should have
03

+

024

=

@‘I

+

b3)

+

@‘2

+

b4),

from formula (7), in which bi is the branch length between the ith OTU and its nearest
interior node and a is the length between two interior nodes. Since a > 0, (D13 + 024)
should be smaller than (Di2 + D34). However, this contradicts formula (A8). Therefore,
OTUs 1 and 3 cannot be neighbors. Similarly, it can be shown that OTUs 1 and 4
are not neighbors. Therefore, only OTUs 1 and 2 (and OTUs 3 and 4) are the neighbors.
For the cases of more than four OTUs, we use the induction principle. Assumingthat OTUs 1 and 2 are true neighbors when S12 is smallest among all &‘s for the caseof N - 1 OTUs, we prove that the same rule applies in the case of N OTUs.
Suppose that S12 is smallest among all &‘s when there are N OTUs. If we ignore
the Nth OTU, OTUs 1 and 2 are, by assumption, neighbors for the remaining N - 1
OTUs. Therefore, there are three possible pairs of neighbors when the Nth OTU isadded: OTUs 1 and 2, OTUs 1 and N, and OTUs 2 and N. From equation (9), we
have
N-l
&N

-

&2

=

2
k=3

[(DIN

+

D2d

-

CD12

+

DN~I/[~(N-

91.

424

Saitou and Nei

C___-FIG. A 1.-Akpossible relationship for four OTUs (1, 2, N, and k). a, b, , bZ, and c are branch lengths

If OTUs 1 and N are neighbors, D 1N= bl -t bN, DZk = b2 + c, D12 = bl + bZ + a, and
DNk = bN + a + c (see fig. Al). Thus, (DIN + Dxk) - (D12 + DNk) = -2a irrespectiveof k, and & - S12 should be negative. This is contradictory to our assumption that
S12is smallest. Therefore, OTUs I, and N are not neighbors. Similarly, it can be shownthat OTUs 2 and N are not neighbors-andthus that OTUs 1 and 2 should be theneighbors. Since we know that our assertion is true for N = 4, it is true for any
N (24).
LITERATURECITEDBUNEMAN,P. 197 1. The recovery of trees from measurements of dissimilarity. Pp. 387-395 in
F. R. HODSON, D. G. KENDALL,and P. TAUTU, eds. Mathematics in the archeological andhistorical sciences. Edinburgh University Press, Edinburgh.
CASE, S. M. 1978. Biochemical systematics of members of the genus Rana native to western
North America. Syst. Zool. 27:299-3 11.
FAITH, D. P. 1985. Distance methods and the approximation of most-parsimonious trees. Syst.
Zool. 34:3 12-325.
FARRIS, J. S. 1972. Estimating phylogenetic trees from distance matrices. Am. Nat. 106:645668.
. 1977. On the phenetic approach to vertebrate classification. Pp. 823-850 in M. K.
HECHT, P. C. GOODY, and B. H. HECHT, eds. Major patterns in vertebrate evolution. Plenum,
New York.
FITCH, W. M. 198 1. A non-sequential method for constructing trees and hierarchical classifications. J. Mol. Evol. l&30-37.
FITCH, W. M., and E. MARGOLIASH. 1967. Construction of phylogenetic trees. Science 155:279-284.
HARTIGAN, J. A. 1973. Minimum mutation fits to a given tree. Biometrics 29:53-65.
JUKES, T. H., and C. R. CANTOR. 1969. Evolution of protein molecules. Pp. 2 l-l 32 in H. N.
MUNRO, ed. Mammalian protein metabolism. Vo13. Academic Press, New York.
KLOTZ, L. C., and R. L. BLANKEN. 198 1. A practical method for calculating evolutionary trees
from sequence data. J. Theor. Biol. 91:261-272.
LI, W.-H. 198 1. Simple method for constructing phylogenetic trees from distance matrices.
Proc. Natl. Acad. Sci. USA 78:1085-1089.
NEI, M. 1987. Molecular evolutionary genetics. Columbia University Press, New York.
ROBINSON,D. F., and L. R. FOULDS. 198 1. Comparison of phylogenetic trees. Math. Biosci.
53:131-147.
SAITOU,N. 1986. Theoretical studies on the methods of reconstructing phylogenetic trees fromDNA sequence data. Ph.D. diss. The University of Texas Health Science Center, Houston.
SAITOU,N., and M. NEI. 1986. The number of nucleotides required to determine the branching

Neighbor-joiningMethod 425
order of three species with special reference to the human-chimpanzee-gorilladivergence. J.
Mol. Evol. 24: 189-204.
SATTATH,S., and A. TVERSKY. 1977. Additive similarity trees. Psychometrika 42:3 19-345.
SOKAL,R. R., and P. H. A. SNEATH. 1963. Principles of numerical taxonomy. W. H. Freeman,
San Francisco.
SOURDIS, J., and C. KRIMBAS. 1987. Accuracy of phylogenetic trees estimated from DNA sequence data. Mol. Biol. Evol. 4: 159-168.
TATENO, Y., M. NEI, and F. TAJIMA. 1982. Accuracy of estimated phylogenetic trees frommolecular data. I. Distantly related species. J. Mol. Evol. l&387-404.
WALTER M. FITCH, reviewing
Received

Augusteditor5, 1986; revision

received

February

18, 1987.