12Phylogenetic relationship of a fossil macaque (Macaca cf. robusta) from the Korean

3

Peninsula to extant species of macaques based on zygomaxillary morphology

45
6

Tsuyoshi Ito1, Yung-jo Lee2, Takeshi D. Nishimura1, Mikiko Tanaka1, Jong-yoon Woo2,
Masanaru Takai1

7

8
9

10
11

1Department

of Evolution and Phylogeny, Primate Research Institute, Kyoto University,

Inuyama, Aichi 484-8506, Japan

2Instituteof Korean Prehistory, 2559, Yongam-dong, Sangdang-gu, Cheongju, Chungbuk

28763, Korea

1213* Corresponding author: Tsuyoshi Ito

14

Telephone: +81-568-63-0523, Fax: +81-568-61-5775
1

15

E-mail: ito.tsuyoshi.3a@kyoto-u.ac.jp

16

17
18

Keywords: Biogeography; Computed tomography; East Asia; Macaca; Phylogentic
morphometrics

19

20

Abstract21

Little is known about the biogeographical and evolutionary histories of macaques

22(Macaca spp.) in East Asia because the phylogenetic positions of fossil species remain

23

unclear. Here we examined the zygomaxillary remains of a fossil macaque (M. cf. robusta)

24

from the Durubong Cave Complex, South Korea, that dates back to the late Middle to Late

25

Pleistocene, to infer its phylogenetic relationship to extant species. We took 195 fixed- and26

semi-landmarks from the zygomaxillary regions of the fossil specimen and from 147

27

specimens belonging to 14 extant species. We then conducted a generalized Procrustes

28

analysis followed by a multivariate statistical analysis to evaluate the phenetic affinities of29

the fossil to the extant species and reconstructed the most parsimonious phylogenetic tree30

using a phylogenetic morphometric approach. We found that the fossil was most similar to2

31

M. fuscata (Japanese macaque) in the zygomaxillary morphospace although it was at the limit

32of the range of variation for this species. The second closest in the morphospace was the33

continental M. mulatta (rhesus macaque). Parsimonious reconstruction confirmed that the34

fossil was most closely related to M. fuscata, even after controlling for the effects of allometry.

35

These findings suggest that in the late Middle to Late Pleistocene, close relatives of M. fuscata

36that looked like the extant species were distributed on the Korean Peninsula, where no species

37

of macaques are found today. Thus, some morphological characteristics of M. fuscata may

38

have developed before its ancestor dispersed into the Japanese archipelago.

39

1. Introduction

40

The genus Macaca consists of approximately 20 extant species and occupies a wider

41

range of climates and habitats than any other genus of nonhuman primates (Fleagle, 2013).

42

Consequently, macaques have attracted much attention as an analogy to understand how

43humans left tropical regions and adapted to various other environments (Rae et al., 2003;

44

Márquez and Laitman, 2008; Hanya et al., 2011; Karen J, 2011; Tsuji et al., 2013; Ito et al.,452015). Furthermore, some East Asian species, particularly M. mulatta (rhesus macaque) and46

M. fuscata (Japanese macaque or snow monkey), are often used in biological laboratory3

47

studies (Sibal and Samson, 2001). To better understand these cold-adapted and well-studied

48

species, we need to understand their biogeographic and evolutionary histories, particularly at49

the northern limit of their distribution in East Asia.

50

The phylogeny of extant species of macaques has been well studied, with most molecular51

phylogenetic studies (Tosi et al., 2000; Li et al., 2009; Jiang et al., 2016), except

52

mitochondrial studies (Morales and Melnick, 1998; Liedigk et al., 2014), supporting the53

classification of Delson (1980). Delson’s classification subdivides macaques into four54

phylogenetic groups: fascicularis, sinica, silenus, and sylvanus groups. A molecular55

phylogenetic study by Jiang et al. (2016) suggested that the African sylvanus group first56

diverged ca. 5.5 million years ago (Mya), the silenus group diverged ca. 4.5 Mya, and the

57

fascicularis and sinica groups diverged ca. 3.5 Mya.

58

The most recently diverged groups (i.e., fascicularis and sinica groups) are distributed

59

in East Asia and partly overlap each other (Fooden, 1988, 2006). The fascicularis group60

consists of four extant species, including M. fuscata and M. mulatta, which are found in61

temperate regions: M. fuscata is distributed in the Japanese archipelago, which is at the62northernmost limit of extant nonhuman primates (Fooden and Aimi, 2005), while M. mulatta

63

is widely distributed from Afghanistan to southern China (Fooden, 2000). The sinica group4

64

consists of at least five (Sinha et al., 2005; possibly six or seven; Chakraborty et al., 2007; Li65

et al., 2015; Fan et al., 2016) extant species, with some of the large-bodied species, such as

66

M. thibetana, M. assamensis, and M. leucogenys, distributed in high-altitude areas in southern

67

China (Fooden, 1982, 1983; Li et al., 2015). Thus, both phylogenetic groups inhabit relatively

68

cold environments in East Asia. However, no species currently inhabit northern China and69

the Korean Peninsula, isolating M. fuscata from the other continental species.

70

However, some excavations of Pleistocene fossil macaques have been reported from71

northern China and the Korean Peninsula (Schlosser, 1924; Young, 1934; Zhang et al., 1986;72

Pan and Jablonski, 1987; Park and Lee, 1998; Takai et al., 2008). The first record of a fossil

73

macaque from East Asia was M. anderssoni, which was described from a nearly complete

74

face excavated from the Early Pleistocene sediment of Mianchi, Henan Province, China

75(Schlosser, 1924). A decade later, another species, M. robusta, was described from a partial

76

maxilla excavated from the Middle Pleistocene sediment of Choukoutien, Beijing, China

77(Young, 1934). Since then, nearly all of the fossil specimens that have been discovered in78

northern China or the Korean Peninsula from the Early to Middle Pleistocene have been79

referred to as one of these two fossil species. Although M. anderssoni has a larger dental size

80

than M. robusta, it is sometimes considered a junior synonym of the latter (Simons, 1970;5

81

Delson, 1980; but see Fooden, 1990). These fossils are believed to fill the biogeographical82

gap in the current distribution of macaques, which will potentially elucidate their evolutionary83

histories. However, the phylogenetic relationships of these fossils and extant species remain

84

controversial (Delson, 1977, 1980; Jablonski and Pan, 1988; Fooden, 1990; Pan and85

Yanzhang, 1995; Park and Lee, 1998; Ito et al., 2014b), particularly with regard to whether

86

they are phylogenetically related to members of the sinica group or the fascicularis group.

87

The fossil records from the Korean Peninsula is key to understanding the evolution of88

M. fuscata in particular. Some fossil macaques have been reported from the Middle and Late

89

Pleistocene sediments of several localities in the Korean Peninsula (Park and Lee, 1998;

90

Fooden and Aimi, 2005; Lee and Woo, 2005; Lee, 2006; Lee and Takai, 2012; Lee et al.,912013). Most of these fossils are isolated teeth or partial fragments of skeletons, but the fossil

92

specimen excavated by the Chungbuk National University Museum Team from the Durubong

93(=Turupong) Cave Complex, Cheongju City, Chungbuk Province, South Korea, is a well-94

preserved zygomaxillary region (Fig. 1). This was identified as M. cf. robusta by Park and95

Lee (1998), who reported that its morphology was intermediate between M. robusta and M.

96

fuscata but more similar to M. robusta with regard to tooth size. Lee and Takai (2012) stated97

that the upper molars of the Korean fossils were relatively large, making them comparable
6

98

with those of M. anderssoni, while the lower molars were relatively small compared with99

other fossil and extant macaques. Furthermore, they also stated that most of the Korean fossils

100

retain some accessory cusps (distoconules, interconulus, 6th cusp, and 7th cusp), which are

101frequently found in the molars of M. fuscata (Lee and Takai, 2012). Fooden and Aimi (2005)

102

demonstrated that measurements of the Korean fossil molars were within the range of

103

variation shown by M. fuscata. Together, these findings imply that the Korean fossils may be

104

related to the ancestors of M. fuscata or the other fossil species, namely M. anderssoni and/or

105

M. robusta; and that more than one lineage may have been distributed in the Korean Peninsula106

during the Middle to Late Pleistocene. However, most previous studies have based their107

conclusions on tooth morphology alone, and no phylogenetic assessment of the fossil species

108

in comparison with the various extant species of macaques has been conducted to date.

109

In this study, we compared the morphology of the zygomaxillary region in the Korean

110

fossil specimen that was discovered in the Durubong Cave Complex with that of extant

111

species of macaques. This is the most well-preserved fossil macaque specimen in Korea.

112

Zygomaxillary morphology is likely to be informative for the reconstruction of phylogenetic

113

relationships, because zygomaxillary shape can be partially differentiated among the four114

phylogenetic groups in macaques (Ito et al., 2014a) and has some phylogenetic signal (Ito et7

115al., 2014a). First, we evaluated the pattern of zygomaxillary variation in extant species and116

estimated the evolutionary and phylogenetic significance of this variation. We then assigned

117

the fossil to the morphospace of the extant species and performed a cladistic analysis to infer

118

the evolutionary and phylogenetic relationships of this fossil to the extant species. To do this,119

we used semi-landmark-based geometric morphometrics to capture the surface topography120

of the zygomaxillary region, because anatomically defined landmarks were limited on the121

fossil. Moreover, we examined the nasal cavity morphology using computed tomography

122(CT), because it has been reported that this reflects phylogenetic relationships (Ito et al.,123

2014b; Nishimura et al., 2014; Ito and Nishimura, 2016). Based on the findings of these124

analyses, we discuss the biogeographic and evolutionary histories of East Asian macaques.

125

2. Materials and methods

126

2.1. Molecular phylogeny

127

The phylogeny that was used for the phylogenetic comparative analyses was estimated128

using 11 mitochondrial and 53 nuclear DNA sequences (Supplementary Online Material

129[SOM] Table S1). The DNA sequences were obtained from the 10KTrees webserver (Arnold

130

et al., 2010; http://10ktrees.fas.harvard.edu/Primates/) and Perelman et al. (2011). The
8

131

sequences were aligned using MAFFT version 7 (Katoh and Standley, 2013) and then132

concatenated using SequenceMatrix software (http://gaurav.github.io/taxondna/). BEAST 2

133

software (Bouckaert et al., 2014) was used for phylogenetic inference using a site model that134

was estimated by the Bayesian approach with the bModelTest package (Bouckaert, 2015).

135

Markov Chain Monte Carlo simulations were performed for 10,000,000,000 generations with136

a sampling frequency of 100,000. Five species belonging to the subtribe Papionina

137(Cercocebus torquatus, Lophocebus aterrimus, Mandrillus sphinx, Papio hamadryas, and138

Theropithecus gelada) were used as an outgroup. The maximum clade credibility tree was139

chosen after removing a 10% burn-in using TreeAnnotator version 2.4.7 (SOM Fig. S1).

140

2.2. Cranial sample

141

For the fossil specimen, we used the zygomaxillary specimen of a fossil macaque stored142

at the Chungbuk Natural History Museum, Chungbuk, South Korea (2ㅜㄷ-12-2169; Fig. 1),

143

which has been partially broken and repaired. The fully erupted third molars and relatively

144

large canines indicate that this individual was an adult male. This specimen originated from145

the Durubong (=Turupong) Cave Complex, Cheongju City, Chungbuk Province, South Korea.

146

The faunal composition of Durubong Cave Complex (Ursus arctosi, Crocuta ultima,

147

Dicerorhinus cf. choukoutiensis, and M. cf. robusta) suggests that this specimen dates back9

148

to the late Middle Pleistocene to the Late Pleistocene period (Park and Lee, 1998), while the149

faunal complex together with the results of pollen analyses suggests that this was a warm

150

(interglacial or interstadial) period (Lee and Woo, 2005).

151

For the extant samples, we used the dry crania with almost-erupted or fully erupted third152

molars from 147 adult males belonging to 14 extant species (Table 1). These specimens are

153

stored at eight institutes in Japan and USA, which are listed in SOM Table S2. Sulawesi

154

macaques were not included in the sample because of their unusual cranial morphology

155(Albrecht, 1978). The samples consisted of wild (n = 80), captive (n = 35), and unknown-

156

origin specimens (n = 32), but any pathological specimens, which showed evidence of severe157

alveolar pyorrhea and/or unusual deflection of muzzles, were not included in the samples.

158

ANOVA showed that the differences in zygomaxillary shape among wild, captive, and159

unknown-origin specimens were significant in some cases, but were much smaller compared160

with interspecies differences (SOM Table S3).

161

2.3. CT examinations

162

CT images were acquired for the fossil specimen and 134 specimens of extant species.

163

The zygomaxillary fragment of the Korean fossil specimen (2ㅜㄷ-12-2169) was scanned

10

164

using a helical scanner (HiSpeed; GE Medical Systems Inc., Waukesha, WI, USA) at the165veterinary medical center of Chungbuk National University, Chungbuk, South Korea, with a166pixel size of 0.1 mm and a slice thickness of 1.0 mm. The original volume data were then167

reconstructed to give a slice thickness of 0.3 mm. CT scans of the extant species were168

obtained previously (Ito and Nishimura, 2016) and are available at the Digital Morphology

169

Museum, Kyoto University (http://dmm.pri.kyoto-u.ac.jp), with the exception of specimens

170

that were not permitted for release due to copyright restrictions (SOM Table S2), or from

171

MorphoSource (http://www.morphosource.org).

172

The internal structure of the maxilla was observed using Amira 5.5. In particular, we173

investigated the degree of lateral expansion of the nasal cavity and the maxillary thickness,

174

which are considered to reflect phylogenetic relationships in macaques (Ito et al., 2014b).

175

2.4. Geometric morphometrics

176

Three-dimensional (3D) surface models of the face were obtained from the Smithsonian177

3D collection of primates (http://humanorigins.si.edu/evidence/3d-collection/primate) or

178

constructed from the serial CT images using Amira 5.5 software (FEI Visualization Sciences

179

Group, Bordeaux, France). Details of the 3D surface modeling were described in a previous

11

180

paper (Ito and Nishimura, 2016). Small superficial holes, if any, were manually filled using

181

Amira 5.5 or Rapidform XOS3 64. Where the left side of the face was broken, the 3D surface

182

model was horizontally flipped using MeshLab software (http://www.meshlab.net). Because

183

the repair of the fossil specimen was slightly imperfect (i.e., the joins of fragments were not184

smooth; Fig. 1), we slightly re-repaired the fossil specimen virtually using Amira 5.5 software185

so that the joins of fragments were smooth and the repaired specimen had a more natural

186

appearance (Fig. 2).

187

In total, 195 landmarks were taken on the 3D surface models (left side of the face) using

188

Stratovan Checkpoint software (Stratovan Corporation, Sacramento, CA, USA; Fig. 3). This

189

included 22 anatomically defined landmarks, 42 semi-landmarks on curves, and 131 semi-

190

landmarks on the surface (SOM Table S4). Missing landmarks (16 landmarks in nine191

specimens) were estimated by deforming the configurations of the weighted average of the192

complete dataset for the same species using a thin-plate spline (TPS) interpolation calculated

193

from the available landmarks. The semi-landmarks were slid on curves or surfaces so that the194bending energy of TPS deformation was minimal (Gunz et al., 2005). These procedures were195

performed using the package “Morpho” version 2.5.1 (Schlager, 2017) in R software version

196

3.4.1 (R Core Team, 2017). Following sliding, the semi-landmarks were treated in the same
12

197

way as the homologous landmarks in subsequent analyses.

198

All of the final coordinates were subjected to generalized Procrustes analysis (GPA),

199

which registered the landmark configurations across all individuals to produce Procrustes

200

coordinates. The centroid size was simultaneously calculated as the square root of the sum of201

the square distances of all landmarks from their centroid. To evaluate non-allometric202

variations, we also produced size-adjusted shape components by calculating the residuals of203

the multiple regression of Procrustes coordinates on natural-logarithmically transformed

204

centroid size (lnCS) using MorphoJ version 1.06d (Klingenberg, 2011), wherein 16.2% of

205

the shape variance was explained by lnCS (p < 0.0001).

206

We then calculated species mean values. To do this, GPA was conducted for each207

species, and the mean shapes of the 15 operational taxonomic units (the 14 extant species and208

the fossil specimen) were then subjected to another GPA. A size-adjusted shape component

209

was then calculated in the same way as mentioned above but on a species mean basis, wherein

210

27.2% of the shape variance was explained by lnCS (p = 0.0012).

211

The individual values were used in the phenetic analyses, while the species mean values

212

were used in the phylogenetic comparative and cladistic analyses. Each of these analyses was

13

213

conducted twice using the raw and size-adjusted shape data.

214

2.5. Phenetic analyses

215

The phenetic affinity of the fossil specimen was evaluated to assign it to the216

morphospace of the extant species. To detect shape vectors that exhibited interspecies

217

differences, between-group principal component analyses (bgPCAs) were performed. The

218

bgPCA is the projection of data onto the principal components of the group means

219(Mitteroecker and Bookstein, 2011). In this study, the species means for extant species were220

included in the calculation of the covariance matrix, following which the individuals of extant

221

species and the fossil specimens were projected onto the bgPCA axes. This step was

222

performed using the “plotTangentSpace” function of the “geomorph” package with some223

modifications. Shape changes along the bgPC axes were visualized using R, Amira 5.5, and224

Rapidform XOS3 64.

225

Procrustes distances were calculated from the fossil specimen to all extant individuals

226

to quantify its phenetic affinity to extant species using the R package “shapes” version 1.2.0

227(Dryden, 2016).

14

228

2.6. Phylogenetic comparative analyses

229

Phylogenetic comparative analyses were performed to evaluate the phylogenetic

230

significance and evolutionary polarity of zygomaxillary shape. A phylomorphospace

231(Sidlauskas, 2008) was obtained to reconstruct the evolutionary changes in zygomaxillary232

shape along the tree using the “plotGMPhyloMorphoSpace” function of the “geomorph”

233

package with some modifications. A phylomorphospace is a PCA of the species mean values

234

and the inferred ancestral shapes, with the phylogeny mapped onto the plane of the PCs. This235

was produced using only the extant species (and the inferred ancestral nodes), following

236

which the fossil specimen was extrapolated into the PC axes. The phylogenetic signal in237

zygomaxillary shape was evaluated using the “physignal” function of the “geomorph”

238

package. The significance of the phylogenetic signal was then tested using permutation

239

procedures with 9,999 iterations.

240

2.7. Cladistic analyses

241

A phylogenetic morphometric approach (Catalano et al., 2010; Goloboff and Catalano,2422011) was used to estimate the phylogenetic relationship between the fossil specimen and the243

extant species examined. This approach parsimoniously detects the ancestral points that

15

244

minimize the distance between the ancestor and descendant points along the tree and has been

245

shown to have a better performance than other landmark-based cladistic approaches

246(Catalano and Torres, 2016). We performed the cladistic analyses using the species mean

247

values of the Procrustes coordinates with the Windows GUI version of TNT 1.5 software

248(Goloboff et al., 2008; Goloboff and Catalano, 2016), setting M. sylvanus as an outgroup and249

searching for the best tree with and without monophyly constraints. Considering previous

250

molecular phylogenetic studies of Alu elements (Li et al., 2009; Jiang et al., 2016), the251

constraints were forced as follows: [([arc, ass, thi], [rad, sin]), ([fus, mul, cyc], fas)], (leo, sil,252

sib, nem); the fossil was allowed to move freely (for abbreviations, see Table 1). Landmark

253

optimization was set as six grid cells and one nesting with iterations. A traditional search was254

conducted, which consisted of Wagner trees with 10 replications followed by tree bisection

255

reconnection. Resampling was performed to evaluate phylogenetic uncertainty, wherein

256

symmetric resampling with 33% change probability was conducted, and the frequency

257

difference (GC values; Goloboff, 2003) was calculated.

16

258

3. Results259

3.1. Internal structures

260

The internal structures of the maxilla in the extant species are shown in Figure 4. The

261

large-bodied species of the sinica group, namely M. arctoides, M. assamensis, and M.

262

thibetana, often exhibited a laterally expanded nasal cavity with thin maxillary walls, which

263

was well represented at the M1 level (this corresponds approximately to the level at which the264nasal cavity shows a maximum width). By contrast, other species, including members of the265

fascicularis, silenus, and sylvanus groups, usually had limited lateral expansion of the nasal

266

cavity and a relatively thick maxillary wall.

267

In the fossil specimen, the lateral and upper walls of the nasal cavity were broken (Fig.

268

5). Nevertheless, based on the surrounding structures, it appeared that the nasal cavity was269

not as greatly expanded laterally as those of the large-bodied species of the sinica group, and270

the maxillary walls were relatively thick.

271

3.2. Phenetic affinities of zygomaxillary shape272

The bgPCA for the raw shape data showed allometric and non-allometric variations in273

the zygomaxillary shape (Fig. 6). bgPC1 and bgPC2 accounted for 47.2% and 18.0% of the
17

274

total interspecies variance, respectively; bgPC1 was positively correlated with lnCS (r = 0.70,

275

p < 0.001), whereas bgPC2 was not (r = 0.10, p = 0.238). Individuals with more positive276

bgPC1 scores (and thus larger lnCS) tended to have muzzles that were longer, superior-

277

inferiorly shorter, and narrower, had a more prominent canine alveolus, and were more

278

ventrally oriented at the anterior portion (Fig. 7). The intercept of this allometric trajectory279

was significantly different among the species groups (ANCOVA, p = 0.010; the sylvanus

280

group was excluded from this analysis due to its small sample size), but their distributions

281

highly overlapped each other. Individuals with more positive bgPC2 scores tended to have282

less prognathic and higher muzzles, with a wider posterior portion compared with the anterior283

portion. The fascicularis and sylvanus groups had high bgPC2 scores compared with the other284

two groups, although their distributions also considerably overlapped.

285

The bgPCA for the size-adjusted shape data exhibited non-allometric variations (Fig. 6),286

with bgPC1 and bgPC2 accounting for 27.2% and 16.4% of this interspecies variance,

287

respectively. In contrast to the bgPCA for the raw shape data, bgPCA for the size-adjusted

288

shape data did not show a variation in facial elongation (Fig. 7). A more positive bgPC1 score289

was associated with a more inferior deflection of the anterior portion of the muzzle against290

the posterior portion. The shape changes, as well as interspecies group differences along the
18

291

bgPC2 axis, were almost identical to those shown for bgPC2 of the raw shape data, but in the292

opposite direction.

293

For the fossil specimen, lnCS and bgPC1 were approximately intermediate among the294

macaques (Fig. 6). The bgPC2 showed that the fossil specimen had a less prognathic and295

higher muzzle, and the posterior portion of the muzzle was wider than the anterior portion,

296

making it similar to members of the fascicularis and sylvanus groups. In particular, the fossil

297

specimen was closest to or within the ranges of variation for M. fuscata and M. mulatta, as

298

seen in the scatterplots of lnCS and bgPC scores. Figure 8 shows the Procrustes distances299

from the fossil specimen to each of the examined individuals of extant species. For both the300

raw and size-adjusted shape data, the fossil specimen was the most similar to M. fuscata and301

the least similar to species of the silenus group. The other species were positioned

302

intermediate between these two extremes, but M. thibetana was similar to the fossil when303

size was adjusted.

304

3.3. Phylogenetic comparisons305

A moderate phylogenetic signal was detected in the raw shape data (K = 0.42, p =306

0.0163) but was marginally significant in the size-adjusted shape data (K = 0.36, p = 0.0627).

19

307

The phylomorphospace demonstrated evolutionary changes in the zygomaxillary shape308(Fig. 9). For the raw shape data, PC1 and PC2 accounted for 50.4% and 17.5% of the variance,309

respectively, and PC1 was positively correlated with lnCS (r = 0.79, p = 0.0009), whereas

310

PC2 was not (r = 0.06, p = 0.831). Shape changes along these PC axes showed similar patterns

311

to the bgPCA of the raw shape data (SOM Fig. S2). For the size-adjusted shape data, PC1

312

and PC2 accounted for 27.7% and 26.4% of the variance, respectively. More positive PC1

313

scores were associated with a more prominent canine alveolus and with the anterior portion314

of the muzzles being deflected more inferiorly against the posterior portion. Shape changes

315

along PC2 were almost identical to those shown in the raw shape data, but in the opposite

316

direction.

317

The fossil specimen was positioned near the branch of M. fuscata for both the raw and318

size-adjusted shape data (Fig. 9). The clade that included M. fuscata, M. mulatta, and M.

319

cyclopis, as well as M. sylvanus, deviated from the ancestral nodes in phylomorphospace.

320

Although the degree of deviation was limited compared with the two lineages, M. thibetana

321

also approached the derived fields of the phylomorphospace, particularly when size was322

adjusted. Such deviations indicated that these lineages independently acquired similar derived

323

features (i.e., a less prognathic and higher muzzle, whose posterior portion was narrower than
20

324

the anterior portion).

325

3.4. Cladistic analyses

326

Cladistic analyses were performed based on the 3D landmark data from the327

zygomaxillary region (Fig. 10). The unconstrained models reconstructed trees whose

328

topologies were quite different from the molecular phylogenetic trees (Li et al., 2009; SOM

329

Fig. S1), indicating that there is some homoplasy in the zygomaxillary morphology. The best330

tree scores, which indicate the sum of the differences between the shapes of the ancestor and331

descendants for all tree branches (Goloboff and Catalano, 2016), were 3.21664 and 3.41040

332

for the raw and size-adjusted shape data, respectively. When the analyses were repeated with333

monophyly constraints, which searched for the most parsimonious scenario within the334

constraints, the best tree scores were 3.70956 and 3.85731 for the raw and size-adjusted shape335

data, respectively. In all four cases, the phylogenetic proximity between the fossil specimen

336

and M. fuscata was supported.

337

4. Discussion338

4.1. Patterns of zygomaxillary variation in extant species

339

We observed allometric variations in zygomaxillary shape among the extant species
21

340

examined. As commonly found in papionin crania (Singleton, 2002; Frost et al., 2003; Gilbert

341

and Grine, 2010; Ito et al., 2014a), we confirmed that larger individuals tended to have more342

elongated and inferiorly deflected muzzles. Such allometric shape components should be

343

interpreted with care in phylogenetic studies because they can obscure phylogenetic signal.

344

For example, although bgPC1 showed that an elongated and inferiorly deflected muzzle is345characteristic of members of the silenus group, this feature is actually just a scaled-up version

346

of that found in small-bodied taxa. Similarly, the shape components that were represented by347

PC1 of the phylomorphospace were also allometric and therefore problematic. Although348

allometric shape components are potentially phylogenetically informative (Gilbert and Rossie,

349

2007; Gilbert et al., 2009), we were unable to detect clear differences in allometric vectors

350

among the four phylogenetic groups. Furthermore, although size itself may provide some351

phylogenetic information, it is inappropriate to rely on a single trait (i.e., size) when inferring352

phylogenetic relationships (see Pearson et al., 2015).

353

Non-allometric shape components were observed in two ways. First, bgPC2 and PC2 of354

phylomorphospace based on the raw shape data were independent of size, which can be

355interpreted as non-allometric. Second, the size-adjusted shape components, which were the356

residuals of a multiple regression of raw Procrustes shape components on lnCS (Monteiro,
22

357

1999; Drake and Klingenberg, 2008; Klingenberg, 2011), were free from allometric358

associations, although they had also partially lost some phylogenetic information in the359

process (see Gilbert, 2011). These non-allometric shape components did not include facial

360

elongation but could be characterized by the degrees of height and prognathism of the muzzle,361

as well as some localized shape changes in the zygomaxillary region.

362

The phylogenetic significance and evolutionary polarity of zygomaxillary shape were363

evaluated. The phylomorphospace suggested that a higher and less prognathic muzzle was a364

derived condition in macaques as shown in PC2, and this was observed in the clade that365

includes M. fuscata, M. mulatta, and M. cyclopis of the fascicularis group, as well as M.

366

thibetana and M. sylvanus. These lineages all inhabit relatively cold environments, and367

therefore this shape change may reflect parallel adaptive evolution to a cold environment,

368

although further testing is required to verify this hypothesis. The degree of change along this369

PC2 axis, however, varies among lineages; M. fuscata is located at the extreme end of this

370

axis (and thus has the highest and least prognathic muzzle), followed by M. sylvanus, close

371

relatives of M. fuscata (M. mulatta and M. cyclopis), and M. thibetana. The difference in the372

degree of the shape change among phylogenetic groups could be explained by phylogenetic

373

constraints, as suggested by the finding that some phylogenetic signal is detected in the
23

374

zygomaxillary shape, although it is marginally significant for size-adjusted shape data.

375

In addition, previous studies have indicated that variations in the nasal cavity can be

376

explained by phylogeny but not by allometry (Ito et al., 2014b; Ito and Nishimura, 2016).

377

Here, the most notable characteristics are a laterally expanded nasal cavity and thin lateral378

maxillary walls at the anterior portion of the muzzle, which are often seen in the three large-379

bodied species of the sinica group (M. arctoides, M. assamensis, and M. thibetana). These

380

internal features are considered to be derived traits in macaques (Ito et al., 2009; Ito et al.,381

2014b). Certainly, nasal cavity morphology also varies in relation to ecological adaptation382

and developmental plasticity. For example, in modern humans, individuals in colder

383

environments tend to have a narrower nasal cavity (Noback et al., 2011), while in Japanese384

macaques, individuals in colder environments tend to have a larger cross-sectional area of the385

cavity (Rae et al., 2003); such clinal variations are generally considered to be a consequence386

of adaptation to climate. Additionally, cold environment-reared rats have smaller nasal

387

cavities than a corresponding control group, which reflects developmental plasticity (Rae et388al., 2006). Nevertheless, it is difficult to explain, solely by such environmental factors, why389

an expanded nasal cavity is often found in the three-large-bodied species of the sinica group390

but is barely observed in other cold environment-inhabiting species such as M. sylvanus and
24

391

M. fuscata. Therefore, the derived traits of internal structures that are found in members of392

the sinica group as well as those of zygomaxillary shape that are typically found in members393

of the fascicularis and sylvanus groups are considered to be informative for estimating the394

phylogenetic relationships of a fossil specimen.

395

4.2. Characteristics of the Korean fossil specimen and their phylogenetic significance396

The Korean fossil specimen was most similar to M. fuscata in several characteristics.

397

However, it must be remembered that a phenetic affinity does not necessarily reflect398

phylogenetic relatedness. The nasal cavity of the Korean fossil specimen exhibited moderate

399

lateral expansion, making it similar to members of the fascicularis group and some of the400

other groups. However, this is not considered to be a derived character (Ito et al., 2014b) and401

so cannot be used to imply their phylogenetic relatedness. The Korean fossil specimen also402

had a somewhat similar level of facial elongation to some other species, such as M. fuscata.

403

However, since this feature is largely associated with size, its reliability for phylogenetic

404

interpretation is also questionable. The long-standing controversies around the phylogenetic

405

relationships of fossil macaques (Delson, 1977, 1980; Jablonski and Pan, 1988; Fooden,

406

1990; Pan and Yanzhang, 1995; Park and Lee, 1998; Ito et al., 2014b) could be partly ascribed407

to the use of such confusing characteristics.
25

408

The Korean fossil specimen has a vertically high and less prognathic muzzle, which is409

the same as the clade that includes M. mulatta, M. cyclopis, and M. fuscata, as well as M.

410

sylvanus. This shape component is independent of allometry and is probably a derived

411

condition in macaques. It is not likely that the Korean fossil is closely related to the African

412

species M. sylvanus. Among the three species of the clade, the Korean fossil was most similar413

to M. fuscata, followed by M. mulatta, with features that lay within the range of intraspecific414

variation for each. Furthermore, cladistic analyses also confirmed phylogenetic relatedness

415

between the Korean fossil specimen and M. fuscata. Although M. thibetana is similar to the416

Korean fossil, particularly for the size-adjusted shape data, the degree of similarity is limited

417

compared with M. fuscata. Therefore, it is reasonable to consider that the Korean fossil

418

specimen is phylogenetically related to the clade that includes M. mulatta, M. cyclopis, and419

M. fuscata, particularly to M. fuscata.

420

It should be noted here that this finding does not necessarily suggest a phylogenetic

421

closeness between members of the fascicularis group and the Chinese fossil species, namely

422

M. robusta and M. anderssoni. Although the fossil macaques discovered from the Korean

423

Peninsula are sometimes referred to as M. robusta (Park and Lee, 1998; Lee and Woo, 2005),

424

some or all of them may be actually unrelated to this fossil species. Further studies are
26

425

expected to elucidate the phylogenetic relationships among the fossil specimens in East Asia.

426

4.3. Biogeographical and evolutionary implications

427

The present study suggests that the Korean fossil specimen is phylogenetically related428

to the fascicularis, rather than the sinica, group of macaques. Delson (1980) and Ito et al.

429(2014b) suggested that the proto-members of the sinica group, for example, M. anderssoni,

430

were largely distributed in northern China in the Early Pleistocene and that members of this431

group may have retreated southward and been replaced by latecomers, that is, members of432

the fascicularis group, likely by the Middle Pleistocene (Delson, 1980; Ito et al., 2014b). This

433

scenario is congruent with the present inference that the macaques living on the Korean

434

Peninsula in the late Middle to Late Pleistocene were related to the fascicularis group.

435

The late Middle to Late Pleistocene period is likely later than the divergence time436

between M. fuscata and M. mulatta, which is considered to be ca. 0.4–1.4 Mya based on437

molecular studies and fossil evidence from Japan (Fooden and Aimi, 2005; Fooden, 2006;

438

Chu et al., 2007; Jiang et al., 2016). Therefore, the Korean macaque is not the direct ancestor

439

of M. fuscata but rather belongs to a collateral line that diverged from the main line of M.

440

fuscata, M. mulatta, or the common ancestral lineage of the two species. They are probably

27

441

remnants that had been retained in the Korean Peninsula. Alternatively, the Korean macaque442

may be a subpopulation of Japanese macaques that had returned from the Japanese443

archipelago. In light of the parsimony principle, the former scenario is more probable. In

444

either case, the Korean lineage could have become extinct due to an inability to escape the445

Korean Peninsula during the Late Pleistocene glacial period.

446

The zygomaxillary shape of the Korean fossil specimen is intermediate between M.

447

fuscata and M. mulatta but more similar to the former. It has previously been suggested that448

the distinct morphological characteristics of M. fuscata developed after its dispersal into the449

Japanese archipelago (Fooden, 2006). However, assuming that the Korean macaque is not the450

lineage that has returned from the Japanese archipelago, it follows that the morphological451

evolution of these features started before this dispersal. These findings fill in some of the452biogeographical and evolutionary gaps between M. fuscata and continental lineages,453

potentially shedding light on cold adaptation in primates.

454

4.4. Limitations of this study

455This study has several important limitations, particularly related to sample composition.

456

First, the sample size is small for both extant (n = 2–53/species) and fossil species (n = 1).

28

457

This limitation risks an underestimation of intraspecific variation and reduces the reliability458

of species mean value estimation (Cardini and Elton, 2007; Schillaci and Schillaci, 2009).

459

Second, the sample includes not only wild-caught specimens but also captive and unknown-

460

origin specimens. We demonstrated that the morphological differences among specimen

461

sources were small compared with interspecies differences, but were significant in some cases.

462

Therefore, our interpretation of the results, which mainly suggested a phylogenetic closeness463

between the Korean lineage and M. fuscata, is not necessarily conclusive. In particular, M.

464

thibetana should remain a candidate for being the sister taxon of the Korean lineage, because

465

M. thibetana is somewhat similar to the Korean fossil in terms of zygomaxillary shape,

466

particularly when size is adjusted. Further studies with a larger sample size, which ideally467

include only wild-caught specimens, will clarify our findings.

468

5. Conclusion

469

This study investigated the zygomaxillary morphology of the fossil macaque from the470

Middle to Late Pleistocene of the Korean Peninsula compared with living macaque species.

471

We found that the Korean fossil is the most similar in zygomaxillary morphology to M.

472

fuscata, followed by M. mulatta, and is phylogenetically related to the clade that includes M.

29

473

mulatta, M. cyclopis, and M. fuscata (particularly to M. fuscata). This finding, unlike

474

traditional scenarios, suggests that the morphological characteristics of M. fuscata could have475

been formed before this group dispersed to the Japanese archipelago, potentially shedding476

light on how primates adapted to a cold environment.

477

Acknowledgments

478

We would like to convey our gratitude to JY Seong, SY Lee, BI Yun, KW Lee, and the479

staff of the Chungbuk Natural History Museum and the Institute of Korean Prehistory for480

their help with accessing the fossil specimens and for providing information about them, and481

C. Dongwoo and Y. Lee for their assistance with CT scanning of the specimens. We also

482

thank Tet. Hayashi, Ter. Hayashi, S. Yachimori, E. Ashida, D. Shimizu, T. Takano, Y.

483

Shintaku, H. Taru, H. Hirotani, H. Takahashi, N. Shigehara, and many field workers for484

collecting/preserving the remains of macaques and/or permitting us to access the collections

485

that they manage. We thank N. Ogihara, S. Kondo, M. Nakatsukasa, Thaung Htike, Zin-

486

Maung-Maung-Thein, and W. Yano for their help with CT scanning of extant species and N.

487

Egi for her help with the 3D surface modeling. We thank three anonymous reviewers and488

editors for their constructive comments on our previous version of this manuscript, which

30

489

considerably improved the clarity of this paper. We thank the Smithsonian’s Division of490

Mammals (Dr. Kristofer Helgen) and Human Origins Program (Dr. Matt Tocheri) for the491scans of USNM specimens used in this research (http://humanorigins.si.edu/evidence/3d-492

collection/primate). These scans were acquired through the generous support of the493

Smithsonian 2.0 Fund and the Smithsonian's Collections Care and Preservation Fund. Lynn

494

Copes, Lynn Lucas, and the MCZ provided access to these scans, originally appearing in

495

Copes and Kimbel (2016) and Copes et al. (2016), funding for the collection of which was

496

provided by NSF DDIG #0925793, and a Wenner-Gren Foundation Dissertation Grant #8102

497(both to Lynn Copes). These scans were downloaded from MorphoSource.org, a web-498

accessible archive for 3D digital data housed by Duke University. This study was funded by

499

the Fujiwara Natural History Foundation (to TI), the Keihanshin Consortium for Fostering500

the Next Generation of Global Leaders in Research (K-CONNEX) (to TI), and JSPS grants-

501

in-aid for Scientific Research (Grant 26650171 to TDN).

502

References

503

Albrecht, G.H., 1978. The Craniofacial Morphology of the Sulawesi Macaques: Multivariate

504

Approaches to Biological Problems. Contributions to Primatology 13,. Karger, Basel.

31

505
506

Arnold, C., Matthews, L.J., Nunn, C.L., 2010. The 10 k Trees Website: A new online resourcefor primate phylogeny. Evolutionary Anthropology 19, 114-118.

507

Bouckaert, R., Heled, J., Kuhnert, D., Vaughan, T., Wu, C.H., Xie, D., Suchard, M.A., Rambaut,

508

A., Drummond, A.J., 2014. BEAST 2: a software platform for Bayesian evolutionary509

analysis. PLoS Computational Biology 10, e1003537.

510
511
512
513
514
515
516
517

Bouckaert, R., 2015. bModelTest: Bayesian site model selection for nucleotide data. bioRxiv,
020792.
Cardini, A., Elton, S., 2007. Sample size and sampling error in geometric morphometric studies
of size and shape. Zoomorphology 126, 121-134.
Catalano, S.A., Goloboff, P.A., Giannini, N.P., 2010. Phylogenetic morphometrics (I): the useof landmark data in a phylogenetic framework. Cladistics 26, 539-549.
Catalano, S.A., Torres, A., 2016. Phylogenetic inference based on landmark data in 41
empirical data sets. Zoological Scripta 46, 1-11.

518

Chakraborty, D., Ramakrishnan, U., Panor, J., Mishra, C., Sinha, A., 2007. Phylogenetic

519

relationships and morphometric affinities of the Arunachal macaque Macaca munzala, a520

newly described primate from Arunachal Pradesh, northeastern India. Molecular521

Phylogenetics and Evolution 44, 838-849.
32

522

Chu, J.H., Lin, Y.S., Wu, H.Y., 2007. Evolution and dispersal of three closely related macaque523

species, Macaca mulatta, M. cyclopis, and M. fuscata, in the eastern Asia. Molecular524

Phylogenetics and Evolution 43, 418-429.

525
526

Copes, L.E., Kimbel, W.H., 2016. Cranial vault thickness in primates: Homo erectus does nothave uniquely thick vault bones. Journal of Human Evolution 90, 120-134.

527

Copes, L.E., Lucas, L.M., Thostenson, J.O., Hoekstra, H.E., Boyer, D.M., 2016. A collection528

of non-human primate computed tomography scans housed in MorphoSource, a repository529

for 3D data. Scientific Data 3, 160001.

530

Delson, E., 1977. Vertebrate paleontology, especially of nonhuman primates, in China. In:

531

Howells, W.W., Tsuchitani, P.J. (Eds.), Paleoanthropology in the People's Republic of

532

China. National Academy of Science, Washington, D.C., pp. 40-65.

533

Delson, E., 1980. Fossil macaques, phyletic relationships and a scenario of deployment. In:

534

Lindburg, D.G. (Ed.), The Macaques: Studies in Ecology, Behavior and Evolution.

535

Nostrand-Reinhold, New York, pp. 10-30.

536

Drake, A.G., Klingenberg, C.P., 2008. The pace of morphological change: historical

537

transformation of skull shape in St Bernard dogs. Proceedings of the Royal Society of

538

London B: Biological Sciences 275, 71-76.
33

539
540

Dryden, I.L., 2016. shapes: Statistical Shape Analysis. R package version 1.1-13. https://cran.rproject.org/package=shapes.

541

Fan, P., Liu, Y., Zhang, Z., Zhao, C., Li, C., Liu, W., Liu, Z., Li, M., 2016. Phylogenetic position542

of the white-cheeked macaque (Macaca leucogenys), a newly described primate from

543

Southeastern Tibet. Molecular Phylogenetics and Evolution 107, 80-89.

544

Fleagle, J.G., 2013. Primate Adaptation and Evolution. 3rd edn. Academic Press, San Diego.

545

Fooden, J., 1982. Taxonomy and evolution of the sinica group of macaques: 3. Species and

546
547
548
549
550
551
552
553
554
555

subspecies accounts of Macaca assamensis. Fieldiana Zoology 10, 1-52.
Fooden, J., 1983. Taxonomy and evolution of the sinica group of macaques: 4. Species accountof Macaca thibetana. Fieldiana Zoology 17, 1-20.
Fooden, J., 1988. Taxonomy and evolution of the sinica group of macaques: 6. Interspecificcomparisons and synthesis. Fieldiana Zoology 45, 1-44.
Fooden, J., 1990. The bear macaques, Macaca arctoides: a systematic review. Journal of
Human Evolution 19, 607-686.
Fooden, J., 2000. Systematic review of the rhesus macaques, Macaca mulatta (Zimmermann,
1780). Fieldiana Zoology 96, 1-180.
Fooden, J., Aimi, M., 2005. Systematic review of Japanese macaques, Macaca fuscata (Gray,
34

556
557
558

1870). Fieldiana Zoology 104, 1-198.
Fooden, J., 2006. Comparative review of fascicularis-group species of macaques (Primates:
Macaca). Fieldiana Zoology 107, 1-44.

559

Frost, S.R., Marcus, L.F., Bookstein, F.L., Reddy, D.P., Delson, E., 2003. Cranial allometry,560

phylogeography, and systematics of large-bodied papionins (primates: Cercopithecinae)

561

inferred from geometric morphometric analysis of landmark data. Anatatomical Record

562

275A, 1048-1072.

563

Gilbert, C.C., 2011. Phylogenetic analysis of the African papionin basicranium using 3-D

564

geometric morphometrics: the need for improved methods to account for allometric effects.

565

American Journal of Physical Anthropology 144, 60-71.

566

Gilbert, C.C., Grine, F.E., 2010. Morphometric variation in the papionin muzzle and the567

biochronology of the South African Plio-Pleistocene karst cave deposits. American

568

Journal of Physical Anthropology 141, 418-429.

569

Gilbert, C.C., Rossie, J.B., 2007. Congruence of molecules and morphology using a narrow

570

allometric approach. Proceedings of the National Academy of Sciences of the United

571

States of America 104, 11910-11914.

572

Gilbert, C.C., Frost, S.R., Strait, D.S., 2009. Allometry, sexual dimorphism, and phylogeny: a
35

573

cladistic analysis of extant African papionins using craniodental data. Journal of Human

574

Evolution 57, 298-320.

575
576
577
578
579
580
581
582

Goloboff, P., 2003. Improvements to resampling measures of group support. Cladistics 19, 324332.
Goloboff, P.A., Catalano, S.A., 2011. Phylogenetic morphometrics (II): algorithms forlandmark optimization. Cladistics 27, 42-51.
Goloboff, P.A., Catalano, S.A., 2016. TNT version 1.5, including a full implementation ofphylogenetic morphometrics. Cladistics 32, 221-238.
Goloboff, P.A., Farris, J.S., Nixon, K.C., 2008. TNT, a free program for phylogenetic analysis.
Cladistics 24, 774-786.

583

Gunz, P., Mitteroecker, P., Bookstein, F.L., 2005. Semilandmarks in three dimensions. In: Slice,584

D.E. (Ed.), Modern Morphometrics in Physical Anthropology. Kluwer Academic

585

Publishers, New York, pp. 73-98.

586

Hanya, G., Menard, N., Qarro, M., Ibn Tattou, M., Fuse, M., Vallet, D., Yamada, A., Go, M.,

587

Takafumi, H., Tsujino, R., Agetsuma, N., Wada, K., 2011. Dietary adaptations of588

temperate primates: comparisons of Japanese and Barbary macaques. Primates 52, 187-

589

198.
36

590

Ito, T., Nishimura, T.D., 2016. Enigmatic diversity of the maxillary sinus in macaques and its

591

possible role as a spatial compromise in craniofacial modifications. Evolutionary Biology

592

43, 414-426.

593

Ito, T., Nishimura, T.D., Senut, B., Koppe, T., Treil, J., Takai, M., 2009. Reappraisal of Macaca

594

speciosa subfossilis from the late Pleistocene of northern Vietnam based on the analysis595

of cranial anatomy. International Journal of Primatology 30, 646-661.

596

Ito, T., Nishimura, T., Takai, M., 2014a. Ecogeographical and phylogenetic effects on

597

craniofacial variation in macaques. American Journal of Physical Anthropology 154, 27-

598

41.

599

Ito, T., Nishimura, T.D., Ebbestad, J.O.R., Takai, M., 2014b. Computed tomography

600

examination of the face of Macaca anderssoni (Early Pleistocene, Henan, northern China):

601

implications for the biogeographic history of Asian macaques. Journal of Human

602

Evolution 72, 64-80.

603

Ito, T., Nishimura, T.D., Hamada, Y., Takai, M., 2015. Contribution of the maxillary sinus to604

the modularity and variability of nasal cavity shape in Japanese macaques. Primates 56,605

11-19.

606

Jablonski, N.G., Pan, Y., 1988. The evolution and palaeobiogeography of monkeys in China.
37

607

In: Whyte, P., Aigner, J., Jablonski, N.G., Taylor, G., Walker, D., Wang, P.X. (Eds.), The

608

Palaeoenvironment of East Asia from the Mid-Tertiary. Centre of Asian Studies, Hong

609

Kong, pp. 849-867.

610

Jiang, J., Yu, J., Li, J., Li, P., Fan, Z., Niu, L., Deng, J., Yue, B., 2016. Mitochondrial genome611

and nuclear markers provide new insight into the evolutionary history of macaques. PLoS

612

ONE 11, e0154665.

613
614

Karen J, W., 2011. Climatic and altitudinal influences on variation in Macaca limb morphology.
Anatomy Research International 2011, 1-18.

615

Katoh, K., Standley, D.M., 2013. MAFFT multiple sequence alignment software version 7:

616

improvements in performance and usability. Molecular Biolology and Evolution 30, 772-

617

780.

618
619

Klingenberg, C.R., 2011. MorphoJ: an integrated software package for geometricmorphometrics. Molecular Ecology Resources 11, 353-357.

620

Lee, Y.-J. (Ed.), 2006. The Paleolithic Culture of Jungwon Region, Korea. Institute for621

Jungwon Culture, Chungbuk National University, Cheongju [Korean with English

622

abstracts].

623

Lee, Y.-J., Takai, M., 2012. The Middle to Late Pleistocene macaque fossils from central Korea.
38

624

In: Drozdov, N.I., Lee, Y.-J., Woo, J.-Y. (Eds.), The 17th International Symposium:

625

Suyanggae and Her Neighbors in Kurtak. Institute of Korean Prehistory and Krasnoyarsk

626

State Pedagogical University, Cheongju, pp. 116–119.

627
628

Lee, Y.-J., Woo, J.-Y. (Eds.), 2005. Paleolithic Sites of the Jungwon Region, Korea. Chungbuk
National History Museum, Chungbuk [Korean with English abstracts].

629

Lee, Y.-J., Woo, J.-Y., Lee, S.-W., Byeong-Il, Y., Park, J.-Y., 2013. Progress Report of Gunang

630

Cave Site, Danyang (IV) – 2011 Year Excavation (5th). Institute of Korean Prehistory

631

Danyang County-Research Report 38 [Korean with English abstract].

632

Li, C., Zhao, C., Fan, P.F., 2015. White‐ cheeked macaque (Macaca leucogenys): A new

633

macaque species from Modog, southeastern Tibet. American Journal of Primatology 77,634

753-766.

635

Li, J., Han, K., Xing, J., Kim, H.S., Rogers, J., Ryder, O.A., Disotell, T., Yue, B., Batzer, M.A.,

636

2009. Phylogeny of the macaques (Cercopithecidae: Macaca) based on Alu elements.

637

Gene 448, 242-249.

638
639
640

Liedigk, R., Roos, C., Brameier, M., Zinner, D., 2014. Mitogenomics of the Old World monkey
tribe Papionini. BMC Evolutionary Biology 14, 176.
Márquez, S., Laitman, J.T., 2008. Climatic effects on the nasal complex: a CT imaging,
39

641

comparative anatomical, and morphometric investigation of Macaca mulatta and Macaca

642

fascicularis. Anatomical Record 291, 1420-1445.

643

Mitteroecker, P., Bookstein, F., 2011. Linear discrimination, ordination, and the visualization

644

of selection gradients in modern morphometrics. Evolutionary Biology 38, 100-114.

645

Monteiro, L.R., 1999. Multivariate regression models and geometric morphometrics: the search

646

for causal factors in the analysis of shape. Systematic Biology 48, 192-199.

647

Morales, J.C., Melnick, D.J., 1998. Phylogenetic relationships of the macaques

648(Cercopithecidae: Macaca), as revealed by high resolution restriction site mapping of649

mitochondrial ribosomal genes. Journal of Human Evolution 34, 1-23.

650

Nishimura, T.D., Ito, T., Yano, W., Ebbestad, J.O.R., Takai, M., 2014. Nasal architecture in

651

Procynocephalus wimani (Early Pleistocene, China) and implications for its phyletic

652

relationship with Paradolichopithecus. Anthropological Science 122, 101-113.

653
654
655
656
657

Noback, M.L., Harvati, K., Spoor, F., 2011. Climate-related variation of the human nasal cavity.
American Journal of Physical Anthropology 145, 599-614.
Pan, R., Yanzhang, P., 1995. Body size of Macaca anderssoni and its offspring. Human
Evolution 10, 283-287.
Pan, Y.R., Jablonski, N.G., 1987. The age and geographical distribution of fossil cercopithecids
40

658

in China. Human Evolution 2, 59-69.

659

Park, S., Lee, Y., 1998. Pleistocene faunal remains from Saekul/Chonyokul at Turupong cave

660

complex with special emphasis on the large mammalian fossils. In: Xu, Q., Lee, Y. (Eds.),

661

International Symposium for the Celebration of Chinese Academician Jia Lanpo's 90th

662

Birthday: Suyanggae and Her Neighbours. Science Press, Beijing, pp. 55-70.

663

Pearson, A., Groves, C., Cardini, A., 2015. The ‘temporal effect’in hominids: Reinvestigating664

the nature of support for a chimp-human clade in bone morphology. Journal of Human

665

Evolution 88, 146-159.

666

Perelman, P., Johnson, W.E., Roos, C., Seuánez, H.N., Horvath, J.E., Moreira, M.A., Kessing,

667

B., Pontius, J., Roelke, M., Rumpler, Y., 2011. A molecular phylogeny of living primates.

668

PLoS Genetetics 7, e1001342.

669
670

R Core Team, 2017. R: A language and environment for statistical computing, R foundationfor Statistical Computing, Vienna.

671

Rae, T.C., Hill, R.A., Hamada, Y., Koppe, T., 2003. Clinal variation of sinus volume in

672

Japanese macaques (Macaca fuscata). American Journal of Primatology 59, 153-158.

673

Rae, T.C., Vidarsdottir, U.S., Jeffery, N., Steegmann, A.T., 2006. Developmental response to674

cold stress in cranial morphology of Rattus: implications for the interpretation of climatic
41

675

adaptation in fossil hominins. Proceedings of the Royal Society of London B: Biological676

Sciences 273, 2605-2610.

677

Schillaci, M.A., Schillaci, M.E., 2009. Estimating the probability that the sample mean is within

678

a desired fraction of the standard deviation of the true mean. Journal of Human Evolution

679

56, 134-138.

680

Schlager, S., 2017. Morpho and Rvcg – Shape Analysis in R. In: Zheng, G., Li, S., Szekely, G.

681(Eds.), Statistical Shape and Deformation Analysis. Academic Press, New York, pp. 217-

682

256.

683

Schlosser, M., 1924. Fossil primates from China. Palaeontologia Sinica Series C 1, 1-16.

684

Sibal, L.R., Samson, K.J., 2001. Nonhuman primates: a critical role in current disease research.

685
686
687

ILAR Journal 42, 74-84.
Sidlauskas, B., 2008. Continuous and arrested morphological diversification in sister clades ofcharaciform fishes: a phylomorphospace approach. Evolution 62, 3135-3156.

688

Simons, E.L., 1970. The deployment and history of Old World monkeys (Cercopithecidae,

689

Primates). In: Napier, J.R., Napier, P.H. (Eds.), Old World Monkeys: Evolution,

690

Systematics, and Behavior. Academic Press, New York, pp. 97-137.

691

Singleton, M., 2002. Patterns of cranial shape variation in the Papionini (Primates:
42

692

Cercopithecinae). Journal of Human Evolution 42, 547-578.

693

Sinha, A., Datta, A., Madhusudan, M.D., Mishra, C., 2005. Macaca munzala: A new species

694

from Western Arunachal Pradesh, Northeastern India. International Journal of695

Primatology 26, 977-989.

696

Takai, M., Maschenko, E.N., Nishimura, T.D., Anezaki, T., Suzuki, T., 2008. Phylogenetic697

relationships and biogeographic history of Paradolichopithecus suskini Trofimov 1977, a

698

large-bodied cercopithecus monkey from the Pliocene of Eurasia. Quaternery International

699

179, 108-119.

700

Tosi, A.J., Morales, J.C., Melnick, D.J., 2000. Comparison of Y chromosome and mtDNA

701

phylogenies leads to unique inferences of macaque evolutionary history. Molecular702

Phylogenetetics and Evolution 17, 133-144.

703
704

Tsuji, Y., Hanya, G., Grueter, C.C., 2013. Feeding strategies of primates in temperate and alpine
forests: comparison of Asian macaques and colobines. Primates 54, 201-215.

705

Young, C.C., 1934. On the Insectivora, Chiroptera, Rodentia and Primates other than

706

Sinanthropus from locality 1 at Choukoutien. Palaeontologia Sinica Series C 8, 1-139.

707

Zhang, Z., Wei, H., Xu, Z., 1986. Mammalian fossils. In: City, M.o.L.P.a.M.o.B. (Ed.),

708

Miaohoushan: A Site of Early Paleolithic in Benxi Country, Liaoning. Wenwu Press,
43

709

Beijing, pp. 35–66 (in Chinese).

710

711

712

Figure legends

713

Figure 1. The zygomaxillary specimen (2ㅜㄷ-12-2169) from the Durubong Cave Complex

714(Site 2), South Korea. (a) Left lateral view, (b) right lateral view, (c) occlusal view, and715(d) dorsal view. The scale is in centimeters.

716

Figure 2. (a) Original surface model and (b) corrected surface model of the Durubong specimen.

717

The arrows denote fixing of the orientation of fragments (green and blue regions) and718

filling of gaps in the main body (red region).

719

Figure 3. Landmarks used in this study. (a) Frontal view, (b) dorsal view, (c) left lateral view,720

and (d) occlusal view. Anatomically defined landmarks: red; curve semi-landmarks:721

blue; and surface semi-landmarks: green (see also SOM Table S4).

722

Figure 4. Coronal computed tomography (CT) images at the level of M1 in extant macaques.

723(a–d) Species of the sinica group. Arrows denote the degree of lateral
44

724

expansion/constriction of the nasal cavity and the thickness of the lateral bony wall of725

the maxilla. The scale is in centimeters. An asterisk denotes the maxillary sinus.

726

Figure 5. Coronal computed tomography (CT) image (bottom left) and corresponding line-

727

drawn diagram at the level of M1 (bottom right) in the Durubong fossil macaque728

specimen (2ㅜㄷ-12-2169). Upper image illustrates where the cross-section is located.

729

The areas filled in gray and the dots on the line diagram indicate the bony structure and730

plaster, respectively. Arrows denote the limited expansion of the nasal cavity and the731

thick lateral bony wall of the maxilla. The scale is in centimeters.

732

Figure 6. Relationships between bgPC scores and the lnCS of the zygomaxillary region in the733

fossil macaque specimen and extant macaque species. (a) bgPC1 based on the raw shape734

data vs. lnCS, (b) bgPC2 based on the raw shape data vs. lnCS, (c) bgPC2 vs. bgPC1

735

based on the raw shape data, and (d) bgPC2 vs. bgPC1 based on the size-adjusted shape736

data. Red (dashed-dotted lines), sinica group; green (dashed lines), fascicularis group;

737

blue (dotted lines), silenus group; yellow (solid lines), sylvanus group; black, the fossil

738

specimen (2ㅜㄷ-12-2169).

739

Figure 7. Zygomaxillary shape changes along the bgPC axes. (a) Image of the fossil specimen

45

740

showing the locations of the polygons illustrated in this figure. (b) bgPC1 and (c) bgPC2

741

based on the raw shape data. (d) bgPC1 and (e) bgPC2 based on the size-adjusted shape742

data. The upper row represents the lateral views of the zygomaxillary region, while the743lower column represents the dorsal views. The negative (left, blue) and positive (right,

744red) extremes, and their superimposition (middle) are shown.

745

Figure 8. Procrustes distances from the fossil macaque specimen (Macaca cf. robusta) to all

746

examined individuals of extant species of macaques based on (a) the raw shape data and747(b) the size-adjusted shape data. Color indicates the species group (see the legend of748

Fig. 6).

749

Figure 9. Phylomorphospace of (a) the raw shape data and (b) the size-adjusted shape data.

750

Color indicates the species group (see the legend of Fig. 6). Gray-filled circles indicate

751

the most recent common ancestor of all macaques.

752

Figure 10. Parsimonious phylogenetic trees inferred from the three-dimensional landmark data753

from the zygomaxillary region of the fossil macaque specimen (2ㅜㄷ-12-2169) and754

extant species of macaques based on (a) the raw shape data with a monophyly constraint,

755(b) the size-adjusted shape data with a monophyly constraint, (c) the raw shape data

46

756

without a monophyly constraint, and (d) the size-adjusted shape data without a757

monophyly constraint. Numbers near the nodes indicate GC support value when the758

landmarks were resampled. Asterisk indicates the fossil specimen.

759

47

a

b

c

d

a

b

a

b

pr

wpa

c

ns

d

bca
bia
lia

lp4
lca lp3

aif
bp4
bp3

dm3

bm2
bm1

bm3

lif pif

lm2

lm3

lm1
gpf

M. arctoides(KUPRI-4227)

M. assamensis(KUPRI-2303)

M. thibetana(KPN-NF1001819)

M. radiata(KUPRI-6079)

M. fuscata(KUPRI-8886)

M. mulatta(KUPRI-4408)

M. cyclopis(DKY-1688)

M. fascicularis(KUPRI-4703)

*

*

●

●

●
●

●●

●●
●● ●
●
●
●●
●●
●●●
●
● ● ●●●●●●●●● ● ●●●●
●
● ●
●
●
●
●
●
●
●
●
●
●

●

●●

●
●
● ●
●● ●● ●●●●● ● ●●●
●
●
●
● ●
● ●●
●
●
●
●
● ●●
●● ●
●
●● ●
●
● ●
●
●
●
●
●
●
●
●●●
●
●
●
●
●
●

bgPC2
0.00

0.05

●●●

● ●
●
● ●
●

−0.05

0.05
bgPC1
0.00

●

5.2

5.3

5.4

5.5
5.6
lnCS

5.7

5.8

●

●

5.2

5.3

5.4

5.5
5.6
lnCS

5.7

5.8

d
0.10

c

●

●

●
●

−0.10

−0.10 −0.05

arc
ass
thi
rad
sin
● fus
mul
cyc
fas
● leo
sil
sib
nem
syl
rob

0.10

b
0.10

a

●

●

0.05

●

●

● ●
●
●

●

●

●

●

bgPC2
0.00

bgPC2
0.00

0.05

●
●
●
●
●
● ● ●● ●●● ●● ●● ●●● ●●●●● ●●●●●
●
●
● ●●●
●
●
●
●●
● ●
●
●
●
●● ●

●

●

●
●

●

●

●

●
●

−0.10

●

●●

●

●

●
● ●
●

●
●

●
●●● ●
●●●
●
● ● ●
●● ●
●
●● ●
● ●●
●●
●●
●
● ● ●
●
●●●
●
● ●
●
●

−0.05

−0.05

●
●

●

●

●

−0.10

−0.05

●

0.00
bgPC1

0.05

0.10

−0.10

−0.05

0.00
bgPC1

0.05

a

b

c

-

bgPC1

+

-

bgPC2

+

-

bgPC2

+

e

d

-

bgPC1

+

a
sylvanus gr.
silenus gr.

fascicularisgr.

sinica gr.

b

syl
nem
sib
sil
leo
fas
cyc
mul
fus
sin
rad
thi
ass
arc
0.06

●
●

●
●

0.10
0.14
Procrustes distance

0.18

syl
nem
sib
sil
leo
fas
cyc
mul
fus
sin
rad
thi
ass
arc
0.06

●
●
●

●

●

●
●

0.10
0.14
Procrustes distance

0.18

b
0.02 0.04

0.10

a

syl

mul
cyc

thi
ass

sin fas

leo
arc

0.05

0.00
PC1

cyc
mul

nem

0.05

fas

sil
ass

PC2
0.02

sib
sil

rad

rad

sin
leo
sib

thi

syl
0.06

0.05

PC2
0.00

0.05

fus

nem
arc

0.10

fus
0.04

0.00 0.02 0.04 0.06
PC1

a

b
100

92

100

＊

94

8

100

100
100
100

100

c

d
17

100

17

100

29

100

61

100

100
93

100

100
78
50

100
100

100
100

100

16

100

100
100
61
100

100

＊

＊

＊

Table 1 The sample used in this study.
Species group

SpeciesAbbreviation

sinica group

M. arctoidesarc

9

M. assamensis

ass

10

M. thibetanathi

5

M. radiatarad8

M. sinicasin

2

M. fuscata

fus

53

M. mulatta

mul

11

M. cyclopis

cyc

10

M. fascicularis

fas

11

M. leonina

leo

2

M. silenussil

4

M. siberu

sib

9

M. nemestrinanem

8

M. sylvanussyl

5fascicularis groupsilenus groupsylvanus groupnunknown

M. cf. robustarob

1Table S2 Specimens used in this study.

Species

Institutea Catalog number ID for analysis

Data sourceb

M. arctoides

DKY

1473

DKY_1473

PRI CT

unknown

JMC

3111

JMC_3111

PRI CT

captive

NMNH

256825

NMNH_256825

Smithsonian 3D collection

wild

PRI

1871

PRI_1871

KAS CT

captive

PRI

2793

PRI_2793

PRI CT

captive

PRI

291

PRI_291

PRI CT

captive

PRI

4227

PRI_4227

KAS CT

unknown

PRI

5846

PRI_5846

PRI CT

captive

PRI

6434

PRI_6434

PRI CT

unknown

PRI

1893

PRI_1893

PRI CT

wild

PRI

2178

PRI_2178

KAS CT

unknown

PRI

2208

PRI_2208

PRI CT

unknown

PRI

2302

PRI_2302

PRI CT

unknown

PRI

2303

PRI_2303

PRI CT

unknown

PRI

3057

PRI_3057

PRI CT

captive

PRI

4498

PRI_4498

PRI CT

unknown

PRI

5015

PRI_5015

PRI CT

unknown

PRI

5230

PRI_5230

KAS CT

unknown

PRI

5342

PRI_5342

PRI CT

unknown

M. assamensis

Doi

OriginM. cyclopis

M. fascicularis

DKY

1688

DKY_1688

PRI CT

unknown

NMNH

296795

NMNH_296795

Smithsonian 3D collection

wild

PRI

1358

PRI_1358

PRI CT

unknown

PRI

4201

PRI_4201

PRI CT

captive

PRI

4233

PRI_4233

PRI CT

unknown

PRI

4486

PRI_4486

PRI CT

wild

PRI

528

PRI_528

PRI CT

captive

PRI

5847

PRI_5847

PRI CT

unknown

PRI

6431

PRI_6431

PRI CT

captive

PRI

9520

PRI_9520

PRI CT

unknown

DKY

995

DKY_995

PRI CT

unknown

NMNH

121511

NMNH_121511

Smithsonian 3D collection

wild

NMNH

317191

NMNH_317191

Smithsonian 3D collection

wild

NMNH

573504

NMNH_573504

Smithsonian 3D collection

wild

PRI

3046

PRI_3046

PRI CT

unknown

PRI

4434

PRI_4434

KAS CT

unknown

PRI

4477

PRI_4477

PRI CT

unknown

PRI

4478

PRI_4478

PRI CT

unknown

PRI

4703

PRI_4703

PRI CT

unknown

PRI

5225

PRI_5225

PRI CT

captive

PRI

6680

PRI_6680

PRI CT

captive

M. fuscata

HMCZ

37709

HMCZ_37709

MorphoSource

doi:10.17602/M2/M3043 wild

HNCC

03-10

HNCC_03-10

PRI CT

wild

HNCC

04-14

HNCC_04-14

PRI CT

wild

HNCC

04-60

HNCC_04-60

PRI CT

wild

HNCC

05-16

HNCC_05-16

PRI CT

wild

HNCC

05-17

HNCC_05-17

PRI CT

wild

HNCC

05-43

HNCC_05-43

PRI CT

wild

HNCC

06-23

HNCC_06-23

PRI CT

wild

HNCC

06-40

HNCC_06-40

PRI CT

wild

HNCC

06-59

HNCC_06-59

PRI CT

wild

HNCC

06-60

HNCC_06-60

PRI CT

wild

PRI

1571

PRI_1571

PRI CT

wild

PRI

2168

PRI_2168

PRI CT

wild

PRI

2586

PRI_2586

PRI CT

wild

PRI

3452

PRI_3452

PRI CT

wild

PRI

4076

PRI_4076

PRI CT

wild

PRI

4113

PRI_4113

PRI CT

wild

PRI

4335

PRI_4335

PRI CT

wild

PRI

4338

PRI_4338

PRI CT

wild

PRI

4341

PRI_4341

PRI CT

wild

PRI

5868

PRI_5868

PRI CT

wild

PRI

6162

PRI_6162

PRI CT

wild

PRI

6470

PRI_6470

PRI CT

wild

PRI

6474

PRI_6474

PRI CT

wild

PRI

6498

PRI_6498

PRI CT

wild

PRI

6504

PRI_6504

PRI CT

wild

PRI

6833

PRI_6833

PRI CT

wild

PRI

7317

PRI_7317

PRI CT

wild

PRI

7381

PRI_7381

PRI CT

wild

PRI

7382

PRI_7382

PRI CT

wild

PRI

7385

PRI_7385

PRI CT

wild

PRI

8644

PRI_8644

PRI CT

wild

PRI

8658

PRI_8658

PRI CT

wild

PRI

8873

PRI_8873

PRI CT

wild

PRI

8885

PRI_8885

PRI CT

wild

PRI

8886

PRI_8886

PRI CT

wild

PRI

8889

PRI_8889

PRI CT

wild

PRI

8892

PRI_8892

PRI CT

wild

PRI

8893

PRI_8893

PRI CT

wild

PRI

8914

PRI_8914

PRI CT

wild

PRI

8974

PRI_8974

PRI CT

wild

PRI

9328

PRI_9328

PRI CTwild

M. leonina
M. mulatta

PRI

9332

PRI_9332

PRI CT

wild

PRI

9340

PRI_9340

PRI CT

wild

PRI

9361

PRI_9361

PRI CT

wild

SNC

H22-170

SNC_H22-170

PRI CT

wild

SNC

H23-077

SNC_H23-077

PRI CT

wild

SNC

H23-285

SNC_H23-285

PRI CT

wild

TPM

M1656

TPM_M1656

PRI CT

wild

TPM

M1889

TPM_M1889

PRI CT

wild

TPM

M3958

TPM_M3958

PRI CT

wild

TPM

M529

TPM_M529

PRI CT

wild

TPM

M659

TPM_M659

PRI CT

wild

AMNH

MO-11090

AMNH_MO-11090 MorphoSource

wild

NMNH

241022

NMNH_241022

Smithsonian 3D collection

wild

DKY

1682

DKY_1682

PRI CT

unknown

HMCZ

26475

HMCZ_26475

MorphoSource

KAS

41

KAS_0041

KAS CT

captive

PRI

218

PRI_218

PRI CT

captive

PRI

2200

PRI_2200

PRI CT

unknown

PRI

223

PRI_223

PRI CT

captive

PRI

224

PRI_224

PRI CT

captive

PRI

242

PRI_242

PRI CT

captivedoi:10.17602/M2/M3052 wild

M. nemestrina

M. radiata

M. siberu

PRI

3523

PRI_3523

PRI CT

unknown

PRI

4408

PRI_4408

PRI CT

captive

PRI

580

PRI_580

PRI CT

captive

DKY

2110

DKY_2110

PRI CT

unknown

NMNH

123144

NMNH_123144

Smithsonian 3D collection

wild

NMNH

154367

NMNH_154367

Smithsonian 3D collection

wild

PRI

1849

PRI_1849

PRI CT

unknown

PRI

2299

PRI_2299

PRI CT

captive

PRI

2454

PRI_2454

KAS CT

unknown

PRI

3055

PRI_3055

KAS CT

unknown

PRI

4225

PRI_4225

PRI CT

captive

JMC

1483

JMC_1483

KAS CT

captive

PRI

3052

PRI_3052

PRI CT

unknown

PRI

6079

PRI_6079

PRI CT

captive

PRI

6684

PRI_6684

PRI CT

captive

PRI

6685

PRI_6685

PRI CT

captive

PRI

6686

PRI_6686

PRI CT

captive

PRI

7140

PRI_7140

PRI CT

captive

PRI

9532

PRI_9532

PRI CT

captive

NMNH

546835

NMNH_546835

Smithsonian 3D collection

wild

PRI

1324

PRI_1324

PRI CTwild

M. silenus

M. sinica
M. sylvanus

M. thibetana

PRI

1338

PRI_1338

PRI CT

wild

PRI

1344

PRI_1344

PRI CT

wild

PRI

1347

PRI_1347

PRI CT

wild

PRI

5072

PRI_5072

PRI CT

wild

PRI

5079

PRI_5079

PRI CT

wild

PRI

5080

PRI_5080

PRI CT

wild

PRI

5095

PRI_5095

PRI CT

wild

JMC

2488

JMC_2488

PRI CT

captive

JMC

4864

JMC_4864

PRI CT

captive

JMC

5740

JMC_5740

PRI CT

captive

NMNH

574135

NMNH_574135

Smithsonian 3D collection

wild

NMNH

15259

NMNH_15259

Smithsonian 3D collection

wild

PRI

1886

PRI_1886

PRI CT

wild

JMC

1392

JMC_1392

PRI CT

captive

JMC

4798

JMC_4798

KAS CT

captive

JMC

5644

JMC_5644

KAS CT

captive

JMC

6330

JMC_6330

PRI CT

captive

NMNH

255979NMNH_255979

Smithsonian 3D collection

wild

AMNH

M-84472

AMNH_M-84472

MorphoSource

wild

JMC

5722

JMC_5722

PRI CT

captive

KPM

NF1001819

KPM_NF1001819

PRI CT

unknown

Unknown (M.

NMNH

241163

NMNH_241163

Smithsonian 3D collection

wild

PRI

4230

PRI_4230

KAS CT

unknown

CNUM

2ㅜㄷ-12-216

CNUM_2-12-2169

CNU CTwild

cf. robusta)
a CNUM, the Chungbuk National University Museum, Chungcheongbuk-do, Korea; PRI, the Primate Research Institute, Kyoto

University, Inuyama, Japan; JMC, the Japan Monkey Centre, Inuyama, Japan; KPM, the Kanagawa Prefectural Museum of Natural
History, Odawara, Japan; DKY, the Department of Anatomy (Macro), Dokkyo Medical University, Mibu, Japan; KAS, the Laboratoryof Physical Anthropology, Graduate School of Science, Kyoto University, Kyoto, Japan; AMNH, the American Museum of Natural
History, New York, USA; NMNH, the National Museum of Natural History, Washington DC, USA; HMCZ, the Harvard Museum of
Comparative Zoology, Cambridge, USA.
b PRI CT, the data scanned at the Primate Research Institute, Kyoto University, Inuyama, Japan; KAS CT, the data scanned at the theLaboratory of Physical Anthropology, Graduate School of Science, Kyoto University, Kyoto, Japan; CNU CT, the data scanned at theChungbuk National University, Chungcheongbuk-do, Korea; the data downloaded from Smithsonian 3D collection(http://humanorigins.si.edu/evidence/3d-collection/primate); the data downloaded from MorphoSource (http://morphosource.org).

Table S3 Test for the differences in zygomaxillary shape among the different origins of specimens.
Response variable

Explanatory variable

R2

SS

Df

F

P

All (wild, captive, and unknown-origin specimens)
bgPC1 of raw shape data

Origin

0.009

0.001

2

0.59

0.556

Species

0.488

0.138

14

8.92

< 0.001

0.145

131

Residuals
bgPC2 of raw shape data

Origin

0.007

0.001

2

0.47

0.628

Species

0.626

0.136

14

15.69

< 0.001

0.081

131

Residuals
bgPC1 of size-adjusted shape data

Origin

0.044

0.005

2

3.02

0.052

Species

0.303

0.045

14

4.07

< 0.001

0.103

131

Residuals
bgPC2 of size-adjusted shape data

Origin

0.039

0.002

2

2.63

0.076

Species

0.670

0.124

14

18.96

< 0.001

0.061

131

Residuals
Procrustes coordinates of raw shape data

Origin

0.015

0.019

2

1.95

< 0.001

Species

0.335

0.422

14

6.13

< 0.001

0.644

131

Residuals
Procrustes coordinates of size-adjusted shape data

Origin

0.021

0.021

2

2.44

< 0.001

Species

0.304

0.314

14

5.15

< 0.001

0.571

131

Residuals

Subset (wild and captive specimens)
bgPC1 of raw shape data

Origin

0.001

0.000

1

0.11

0.739

Species

0.579

0.125

14

9.81

< 0.001

0.091

100

Residuals
bgPC2 of raw shape data

Origin

0.001

0.000

1

0.08

0.778

Species

0.712

0.125

14

17.68

< 0.001

0.051

100

Residuals
bgPC1 of size-adjusted shape data

Origin

0.046

0.003

1

4.83

0.030

Species

0.377

0.038

14

4.31

< 0.001

0.063

100

Residuals
bgPC2 of size-adjusted shape data

Origin

0.025

0.001

1

2.59

0.111

Species

0.732

0.116

14

19.49

< 0.001

0.043

100

Residuals
Procrustes coordinates of raw shape data

Origin

0.013

0.012

1

2.763

< 0.001

Species

0.398

0.382

14

6.184

< 0.001

0.441

100

Residuals
Procrustes coordinates of size-adjusted shape data

Origin

0.017

0.013

1

3.344

< 0.001

Species

0.360

0.280

14

5.037

< 0.001

0.397

100

Residuals

ANOVA (Type II) was performed for the bgPC scores and Procrustes coordinates, wherein origin of specimen (captive, wild, andunknown) and species were used as explanatory variables. This was done for all sample and for the subset of sample (wild andcaptive). For these analyses, "procD.lm" function of "geomorph" package and "etasq" function of "heplots" package in R were used.

Table S4 Landmarks used in this study.
Abbreviation

Definition

Anatomically-defined landmarks
wpa

Point corresponding to largest width of piriform aperture.

ns

Nasospinale.

aif

Antero-sperior point of incisive foramen.

pif

Posterior-most point of incisive foramen.

lif

Lateral-most point of incisive foramen.

pr

Prosthion.

bia

Postero-labial point of lateral incisor alveolus.

lia

Postero-lingual point of lateral incisor alveolus.

bca

Antero-labial point of canine alveolus.

lca

Antero-lingual point of canine alveolus.

bp3

Labial P3: most mesial point on P3 alveolus, projected labially onto alveolar margin.

lp3

Lingual P3: most mesial point on P3 alveolus, projected lingually onto alveolar margin.

bp4

Labial P4: most mesial point on P4 alveolus, projected labially onto alveolar margin.

lp4

Lingual P4: most mesial point on P4 alveolus, projected lingually onto alveolar margin.

bm1

Labial M1: contact point between P4 and M1, projected labially onto alveolar margin.

lm1

Lingual M1: contact point between P4 and M1, projected lingually onto alveolar margin.

bm2

Labial M2: contact point between M1 and M2, projected labially onto alveolar margin.

lm2

Lingual M2: contact point between M1 and M2, projected lingually onto alveolar margin.

bm3

Labial M3: contact point between M2 and M3, projected labially onto alveolar margin.

lm3

Lingual M3: contact point between M2 and M3, projected lingually onto alveolar margin.

dm3

Distal M3: most posterior point on M3 alveolus, projected labially onto alveolar margin.

gpf

Most posterior point on the margin of greater palatine foramen.

Semi-landmarks on curves
Eight points on the superior midline of nasal bone.
Six points on the superior midline of premaxilla, which is between ns and pr.
Seven points on the inferior midline of maxilla, which is from pif to approximately the meeting point of maxilla andpalatine.
Five points on the lateral margin of nasal aperture, which is between wpa and ns.
Ten points on the inferior margin of orbital opening.
Seven points on the inferior margin of zygomatic process.

Semi-landmarks on surface
One-hundred and thirty-two points on the lateral side of zygomaxilla.
Some of these definitions are in accordance with Cardini et al. (2007) and Ito and Nishimura (2016).

1
2
3
4

Figure S1.

Most likely phylogenetic tree from the BEAST2 analysis for macaques and outgrouptaxa using mtDNA and nuclear sequences. Blue horizontal bars represent the posterior 95% CI for the node
ages.

5

1

a

b

-

PC1

+

-

PC2

+

PC2

+

d

c

-

PC1

+

-

6
7
8
9

Figure S2.

Shape changes along PC axes in the phylomorphospace. (a) PC1 and (b) PC2 based on

the raw shape data. (c) PC1 and (d) PC2 based on the size-adjusted shape data. For more details, see thelegend of Figure 6 in the main text.

2