A chiral uranyl-Kemp's tricarboxylate cubic framework: structuredirecting effect of counterions with three-fold rotational symmetry † Pierre Thuéry* a and Jack Harrowfield* b

In the presence of PPh3Me + cations, Kemp's tricarboxylate (kta 3-) complexes the uranyl cation to give [PPh3Me][UO2(kta)] (1), a triperiodic framework with cubic symmetry and srs topology. The cation is held by weak interactions into cavities with matching threefold rotational symmetry. Comparison with the diperiodic hemihydrate polymorph previously reported points to the disrupting role of OHO hydrogen bonds in the latter.

Among the extremely varied carboxylate-containing ligands which have been used to generate uranyl-based coordination polymers and frameworks, 1 the tricarboxylate derived from Kemp's triacid, cis,cis-1,3,5-trimethylcyclohexane-1,3,5-tricarboxylic acid (H3kta), is particularly appealing, having given several metallacyclic rings, cages and nanotubular complexes, besides more usual mono-and diperiodic polymers. 2 The most conspicuous absence in the series of complexes so far obtained is that of triperiodic frameworks. Adoption of a chair conformation of its six-membered ring by kta 3-provides a ligand of three-fold rotational symmetry which, in the conformation where the carboxylate groups are equatorially oriented, appears ideally suited to the formation of crystals with a diperiodic uranyl ion complex array of honeycomb (hcb) topology where the uranyl ion centres have a hexagonal-bipyramidal UO2( 2 O,O'-carboxylate)3 form. In contrast, the convergent, chair triaxial form of the ligand favours the formation of oligomeric closed species. Several honeycomb networks have also been found with the ligand in the boat diequatorial, monoaxial conformation, 2c-e and the hcb topology is thus the most common in uranyl complexes with this ligand. In the 1:1 UO2 2+ /kta 3- stoichiometry found in these networks, the complex is anionic and therefore must be accompanied by a countercation which is of course another influence within the structure. For structure-directing species incapable of direct coordination to carboxylate groups, their hydrogen bonding capacity appears to have a strong influence upon the form of their associated [UO2(kta)] -species, so that the present isolation of anhydrous [PPh3Me][UO2(kta)] (1) has, by comparison to its previously characterized hemi-hydrate (2), 2e provided an opportunity to assess the influence of strong (OH) and weak (CH) hydrogen bond donors upon the structure of the associated anionic metallapolymer. It is of course expected that the inherent difference between OHO and CHO interactions should be large but numerical factors do also need to be taken into account. In the event, the differences between the two structures are striking. Further, the differences between the two unsolvated compounds 1 and [PPh4][UO2(kta)] (3) 2e are equally striking and provide some indication of the finer detail involved in aliphatic and aromatic CHO interactions.

The complex [PPh3Me][UO2(kta)] (1) was synthesized ‡ under solvo-hydrothermal conditions similar to those giving complex 2, but for the addition of neodymium nitrate hexahydrate, intended to provide a heterometallic complex but eventually not included in the final product. Although difficult to rationalize, since nothing is known of the solution equilibria which may influence the solubility of any particular species, the effect of this modification in the synthetic conditions on the crystal structure is considerable, notwithstanding the fact that the formulas of 1 and 2 differ by only half a solvent water molecule.

Complex 1 crystallizes in the Sohncke 3 cubic space group P213, § with the uranium atom, located on a three-fold rotation axis (Wyckoff position 4a), in a hexagonal-bipyramidal environment (Figure 1) and it provides a not uncommon example of a uranyl-carboxylate polymer forming a chiral crystal. The kta 3-ligand, in the chair, triequatorial conformation, has also three-fold rotational symmetry and its three carboxylate groups are  2 O,O'-chelating. It is notable that the carboxylate groups are much closer to perpendicular to the mean plane of the ring than in 2 and 3 [dihedral angles 76.03 (12)° in 1 versus 24.4(3)-40.7(2)° in 2, and 21.5(8)-64.4(3)° in 3; Figure S1], the shape of the ligand being thus farther from discoidal. Despite the fact that both metal and ligand are three-coordinated (3-c) nodes as in the hcb networks previously reported, notably in 2, 2e the uninodal coordination polymer formed is here triperiodic and it has the point symbol {10 3 } and the srs (SrSi2) topological type 4 (Figures 2a andS2). This common topology has previously been found in uranyl complexes with glutarate 5 or cis-1,4-cyclohexanedicarboxylate, 6 with uranium atoms only as 3-c nodes. The srs net is one of the basic 3-c triperiodic nets (while the hcb net is the simplest 3-c diperiodic net) and it is noted to be the only one with vertices having three-fold rotational symmetry. 4 Although this topology frequently leads to interpenetration of enantiomorphic nets, this is prevented here by the bulky PPh3Me + cations which occupy all available spaces, the Kitaigorodski packing index given by PLATON 7 being 0.71 (0.35 with cations excluded). The cations are nestled within cavities centered on the three-fold axes (Figures 2b andS3), which entails CHO hydrogen bonds involving the cation methyl group and three carboxylate oxygen atoms [CO, 3.318(2) Å; C-HO, 158°] or aromatic protons and oxo groups, as well as CH(kta) interactions. The interactions between the methyl group of the cation and three carboxylate groups, in particular, are related to the tilting of the latter with respect to the mean cyclohexyl plane. There is no direct interaction between the PPh3Me + cations, all PP distances being larger than 9.8 Å. Although diperiodic uranyl polymer sheets orthogonal to [111] involving fused, 72-membered clover-shaped nona-uranacycles (Figure S4) can be discerned, only three of the nine equatorial UO6 entities in a given ring lie parallel to this plane, the others being involved in coordination providing crosslinks to sheets above and below to give the full triperiodic array.

While in solution both kta 3-and PPh3Me + appear, on the basis of their 1 H NMR spectra, 8 to have C3v symmetry, motional restrictions in the solid state can reduce this symmetry to C3 and this is the case in the crystal of 1. For kta 3-, this is due to the fact that the OCO planes of the carboxylate groups are tilted in the same sense with respect to the mean plane of the cyclohexane ring (Figure S5), the OO segments defining a left-handed helical array about the C3 axis in the case of the particular crystal studied. In this same crystal, the cation also shows C3 symmetry but here the CC vectors shown in Figure S5 define a right-handed helical array about the axis. While PPh3Me + in 2 is chiral and has a conformation close to one of C3 symmetry, the kta 3-ligand has a symmetry close to C3v and the complete crystal is achiral (space group Pbca). In this case, the three-fold symmetry of the cation is not reflected in its interactions beyond dispersion, a difference evident on the Hirshfeld surfaces generated with CrystalExplorer 9 , shown in Figure 3. These interactions of an unsymmetrical nature involve here also CHO hydrogen bonds between the cation methyl group and carboxylate oxygen atoms or between aromatic protons and oxo or carboxylate groups, as well as CH(kta) interactions. The result is that each cation is associated with three [UO2(kta)] -units, each one in a different way, rather than six, in three symmetric pairs, as in 1. This perturbation of the interactions with respect to those in the anhydrous crystal is seemingly associated with the presence in 2 of the water molecule which bridges, through hydrogen bond donation, two carboxylate oxygen atoms on every second uranyl centre. Interestingly, this presumably strong hydrogen bonding does not appear to inhibit further CH-donation to the same two carboxylate oxygen atoms. "Symmetry matching" of the components of a solid is clearly seen in the structure of 1 to be a significant influence on that structure, but its operation there passes through weak interactions which are not of a unique nature and can readily accommodate minor variations in the components of a crystal. The striking difference in structure between 1 and its hemi-hydrate 2 thus seems to have rather subtle origins, disruption of the symmetric interactions found in 1 by the additional water molecule in 2 being possibly a major influence.

The PPh4 + cation in 3 has an irregular chiral form and while it can in principle adopt a propeller form, this cannot have C3 symmetry (except if one aromatic ring is disordered), so that the ordered cation is symmetrically "incompatible" with any three-fold-symmetric form of the kta 3-ligand, which in fact in 3 has a conformation close to one with a plane of symmetry (Cs), meaning that here the chirality of the crystal (space group P21) is a reflection of that of the cation. Interactions beyond dispersion of the countercation with the hcb diperiodic polymer involve several CHO hydrogen bonds between aromatic protons and oxo or carboxylate groups making a rather different pattern to those seen for the other two complexes. It is noteworthy that formation of the hcb form of [UO2(kta)] -does not require the presence of counterions with three-fold rotational symmetry and the network can be considerably distorted with respect to this ideal symmetry, in particular when kta 3-adopts the boat conformation. 2e Being observed with counterions very different from one another, 2e the hcb topology may reflect a genuine preference of the uranyl ion for complexation in planar arrays. Complexes 1 and 3 provide examples of the class of chiral crystals generated not from chiral components but from configurationally labile species which can adopt chiral forms not separately observable in solution.

Although the triaxial form of the anions of Kemp's triacid in the chair conformation is particularly prone to give closed, polynuclear complexes with uranyl ions, 2e hcb diperiodic networks are the most common species of coordination polymer found with the triequatorial form. Complex 1 is the first example of a triperiodic framework obtained in this family, and it has the srs topology which retains the three-fold rotational symmetry of the ideal hcb net. The transition from the srs network in 1 to the hcb in the hemi-hydrate 2 is accompanied by a disruption of the three-fold symmetric weak interactions between polymer and counterions seen in 1 by the hydrogen bonded water molecule in 2. While multiple weak (CHO) interactions can provide an efficient means for the transmission of chirality from one component to another, the apparently subtle compositional variation associated with conversion from an anhydrous material to a hemi-solvate has dramatic consequences indicative of the dominance of classical over weak hydrogen bonding.

Conflicts of interest

There are no conflicts of interest to declare.

Notes and References ‡ Synthesis of complex 1. H3kta (13 mg, 0.05 mmol), [UO2(NO3)2(H2O)2]•4H2O (25 mg, 0.05 mmol), [Nd(NO3)3]•6H2O (22 mg, 0.05 mmol), PPh3MeBr (18 mg, 0.05 mmol), N,Ndimethylformamide (0.2 mL), and demineralized water (0.6 mL) were placed in a 10 mL tightly closed glass vessel and heated at 140 °C under autogenous pressure in a sand bath, giving a few light yellow crystals of complex 1 within one week. § The data were collected at 100(2) K on a Bruker D8 Quest diffractometer using an Incoatec Microfocus Source (IS 3.0 Mo) and a PHOTON III area detector, and operated with APEX3. 10 The data were processed with SAINT, [START_REF]SAINT, ver. 8.40A; Bruker Nano[END_REF] and empirical absorption corrections were made with SADABS. [START_REF]SADABS, ver[END_REF] The structure was solved by intrinsic phasing with SHELXT, [START_REF] Sheldrick | [END_REF] and refined by full-matrix least-squares on F 2 with SHELXL, 14 using the ShelXle interface. 15 Drawings were made with ORTEP-3 16 and VESTA, 17 and the topological analysis with ToposPro. 18 Crystal data for 1: C31H33O8PU, M = 802.57, cubic, space group P213, a = 14.2093(8) Å, V = 2868.9(5) Å 3 , Z = 4. Refinement of 127 parameters on 2929 independent reflections out of 106653 measured reflections (Rint = 0.053) led to R1 = 0.010, wR2 = 0.025, min = -0.39, max = 0.50 e Å -3 , Flack parameter = 0.000(4).

Fig. 1

 1 Fig.1(a) View of complex 1. Displacement ellipsoids are drawn at the 50% probability level and hydrogen atoms are omitted. Symmetry codes: i = y + 1/2, 1/2 -z, 1 -x; j = 1 -z, x -1/2, 1/2 -y; k = y, z, x; l = z, x, y. Selected bond lengths (Å) and angles (°): U1-O1 1.779(3), U1-O2 1.781(3), U1-O3 2.4498(14), U1-O4 2.4757(15), O1-U1-O2 180, O3-U1-O4 52.69(5), O4-U1-O3 i 67.59(5). (b) and (c) Two views of the framework, down [100] or[111] (parallel to the C3 axis), with uranium coordination polyhedra colored yellow and counterions omitted.

Fig. 2

 2 Fig. 2 (a) Nodal representation of the framework with srs in 1 down a direction slightly rotated from [100] (yellow, uranium nodes; blue, ligand nodes). (b) Space filling representation of the inclusion of the PPh3Me + cation within the cavity centered on the three-fold rotation axis. Yellow, uranium; red, oxygen; blue, carbon in kta 3-; green, all atoms of the counterion. Hydrogen atoms are omitted.

Fig. 3

 3 Fig. 3 Hirshfeld surface of the PPh3Me + cation mapped with dnorm in complexes 1 (a) and 2 (b) viewed down the P-CH3 bond. The red spots near the methyl group correspond to CHO hydrogen bonds.


Nicolas Brucato
email: nicolasbrucato@gmail.com

Mathilde André

Roxanne Tsang

Lauri Saag

Jason Kariwiga

Kylie Sesuki

Teppsy Beni

William Pomat

John Muke

Vincent Meyer

Anne Boland

Jean- François Deleuze

Herawati Sudoyo

Mayukh Mondal

Luca Pagani

Irene Gallego Romero

Mait Metspalu

Murray P Cox

Matthew Leavesley

François-Xavier Ricaut
email: francois-xavier.ricaut@univ-tlse3.fr

Papua New Guinean genomes reveal the complex settlement of north Sahul

Keywords: Papuan, Human genome, Demographic history, Sahul, Oceania

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

New Guineans are descendants of one of the most exceptional and ancient migrations in human history, but the nature of their initial settlement of New Guinea is still veiled in mystery [START_REF] Allen | A different paradigm for the initial colonisation of Sahul[END_REF]. When modern human first arrived in Southeast Asia 50-70thousand years ago (kya) [START_REF] Westaway | An early modern human presence in Sumatra 73,000-63,000 years ago[END_REF], the lower sea level during a colder climate caused two gigantic ancient continents to emerge: Sunda joined many current southeast Asian islands with mainland Asia, while Australia and New Guinea formed a single landmass called Sahul. These vast territories were separated by a substantial water crossing, Wallacea, containing a scattering of small islands and creating the major ecological barrier known as the Wallace line (Wallace 1869).

The first anatomically modern humans to cross from Sunda to Sahul numbered probably no more than a few hundred, but migration through the Wallacean islands was nonetheless likely a deliberate effort (Bird, et al. 2019;Bradshaw, et al. 2019). Several initial settlement paths are possible based on paleogeological reconstructions: a northern route through what today are the Maluku islands [START_REF] Kealy | Least-cost pathway models indicate northern human dispersal from Sunda to Sahul[END_REF][START_REF] Norman | An early colonisation pathway into northwest Australia 70-60,000 years ago[END_REF], and a southern route through Flores and Timor [START_REF] Birdsell | The recalibration of a paradigm for the first peopling of Greater Australia. Sunda and Sahul: Prehistoric Studies in Southeast Asia[END_REF][START_REF] Norman | An early colonisation pathway into northwest Australia 70-60,000 years ago[END_REF]. Although not mutually exclusive, for decades the consensus has been that a single migration from Sunda to Sahul created the entire biological diversity of both New Guineans and Indigenous Australians today [START_REF] Malaspinas | A genomic history of Aboriginal Australia[END_REF][START_REF] O'connell | When did Homo sapiens first reach Southeast Asia and Sahul?[END_REF]). Both populations have unique genomic diversities marked by their high percentage of genetic inheritance coming from Denisovans (3-3.5%) (Jacobs, et al. 2019;[START_REF] Reich | Denisova Admixture and the First Modern Human Dispersals into Southeast Asia and Oceania[END_REF]Sankararaman, et al. 2016;[START_REF] Vernot | Excavating Neandertal and Denisovan DNA from the genomes of Melanesian individuals[END_REF]. However, recent studies now suggest that the arrival of just one human group is unlikely to be sufficient to explain the deep genetic diversity still observed in this region [START_REF] Allen | A different paradigm for the initial colonisation of Sahul[END_REF][START_REF] Pedro | Papuan mitochondrial genomes and the settlement of Sahul[END_REF]. In part, this argument draws from the uniparental genetic markers found in New Guineans and Indigenous Australians. While both populations undoubtedly share a common maternal genetic ancestry, marked by the presence of haplogroups M*,N* and R*, the diversity of their respective mitochondrial DNA (mtDNA) sub-haplogroups indicate stark genetic differences [START_REF] Redd | Peopling of Sahul: mtDNA variation in aboriginal Australian and Papua New Guinean populations[END_REF].

Most mtDNA sub-haplogroups are exclusive to either New Guinea (e.g., P1, P2, P4a, M27, M28, Q*) [START_REF] Pedro | Papuan mitochondrial genomes and the settlement of Sahul[END_REF] or Australia (e.g., P8, P4b, M42, S*) [START_REF] Tobler | Aboriginal mitogenomes reveal 50,000 years of regionalism in Australia[END_REF], in agreement with local evolution of some of these region-specific haplogroups. Coalescence dates between Australian and New Guinean maternal lineages suggest a split at the time of the initial settlement of Sahul [START_REF] Tobler | Aboriginal mitogenomes reveal 50,000 years of regionalism in Australia[END_REF], or potentially before, as shown by the coalescence dates of P2'P8'P10 (49kya), P4 (50kya) [START_REF] Pedro | Papuan mitochondrial genomes and the settlement of Sahul[END_REF]. Similar patterns are observed for paternal genetic lineages on the Y chromosome (Bergstrom, et al. 2016;Nagle, et al. 2016). Such different uniparental genetic patterns suggest that either New Guineans and Indigenous Australians separated very early after the initial migration to Sahul and that uniparental lineages segregated accordingly, or that both populations descend from two initially separate migrations [START_REF] Redd | Peopling of Sahul: mtDNA variation in aboriginal Australian and Papua New Guinean populations[END_REF], perhaps at different times as suggested by some archaeological records [START_REF] Clarkson | Human occupation of northern Australia by 65,000 years ago[END_REF].

Once in Sahul, modern human dispersals within what is now the island of New Guinea began by at least 49kya [START_REF] Summerhayes | Human adaptation and plant use in highland New Guinea 49,000 to 44,000 years ago[END_REF]. Given the paucity of archaeological records, how New Guinea was settled remains poorly resolved. Of particular note, the initial entry point into New Guinea has proven particularly difficult to locate, leading to alternative hypotheses of migration via the New Guinean north coast, through the mountainous highlands, along the southern foothills, or even further south via a vast plain that is now submerged and forms the Arafura Sea [START_REF] Summerhayes | The archaeology of forest exploitation and change in the tropics during the Pleistocene: The case of Northern Sahul (Pleistocene New Guinea)[END_REF]. Archaeological records perhaps favor the Arafura plain pathway, since the oldest known archaeological sites in New Guinea are located in its southeast region. The suggestion is that even earlier sites, marking those first movements, now lie flooded under the Arafura Sea [START_REF] Summerhayes | The archaeology of forest exploitation and change in the tropics during the Pleistocene: The case of Northern Sahul (Pleistocene New Guinea)[END_REF].

Nevertheless, the alternative hypotheses remain feasible. Although the highlands were not hospitable during the Upper Pleistocene, long-term human presence is attested there from at least 25kya [START_REF] Summerhayes | The archaeology of forest exploitation and change in the tropics during the Pleistocene: The case of Northern Sahul (Pleistocene New Guinea)[END_REF].Major rivers carved into the New Guinean geology may have formed early natural corridors during initial settlement, allowing human mobility as they still do along the Markham and Ramu Rivers in the northeast, the Fly River in the south and the Sepik River in the north [START_REF] Pawley | Papuan pasts: Cultural, linguistic and biological histories of Papuan-speaking peoples: Pacific Linguistics, Research School of Pacific and Asian Studies, The[END_REF]. The Sepik route has been supported by close genetic relationships between modern New Guinea highlanders and lowland groups from the Sepik region, but this association may date to a later time period [START_REF] Bergstrom | A Neolithic expansion, but strong genetic structure, in the independent history of New Guinea[END_REF].

After the initial settlement of New Guinea, independent development of agricultural practices in New Guinea between 10-7 kya led to a unique Neolithic revolution [START_REF] Denham | Origins of agriculture at Kuk Swamp in the highlands of New Guinea[END_REF][START_REF] Golson | Ten thousand years of cultivation at Kuk Swamp in the highlands of Papua New Guinea[END_REF], driving demographic expansion in the highlands that heavily impacted the genetic and linguistic diversity of these populations [START_REF] Bergstrom | A Neolithic expansion, but strong genetic structure, in the independent history of New Guinea[END_REF][START_REF] Greenhill | TransNewGuinea.org: An Online Database of New Guinea Languages[END_REF][START_REF] Reesink | Explaining the linguistic diversity of Sahul using population models[END_REF]. By the time Austronesian settlers arrived around 3kya [START_REF] Bellwood | Prehistory of the Indo-Malaysian Archipelago[END_REF][START_REF] Lipson | Reconstructing Austronesian population history in Island Southeast Asia[END_REF], New Guineans had already developed a complex network of diverse cultures [START_REF] Pawley | Papuan pasts: Cultural, linguistic and biological histories of Papuan-speaking peoples: Pacific Linguistics, Research School of Pacific and Asian Studies, The[END_REF].

Given its extraordinary diversity today, comparable to large continental areas elsewhere, the genetics of New Guinean populations are poorly explored, and scenarios of initial and postsettlement movements are lacking. We have therefore generated 58 high quality whole genome sequences (WGS) from Papua New Guinean individuals (Supplementary table 1).

Coupling these with a dataset of 905 genomes from diverse human groups, but particularly individuals from neighboring Indonesia and Australia (Supplementary table 2), we reconstruct the key genetic events during the Upper Pleistocene that led to the settlement of New Guinea. Our two primary questions are: was the migration to Sahul conducted by one or several human groups, and what migratory routes did they take? And how did New Guineans subsequently settle the diverse environments of the modern island, and what were the dynamics of this settlement?

Results

Samples and diversity

Fifty-eight Papua New Guinean participants were included in this study to increase the geographical coverage of New Guinea, especially in the southern regions underrepresented in published whole genome datasets (Supplementary figure 1 andsupplementary table 1).

Whole genome sequencing data were obtained at an average depth of 36x, characterized for 41,250,479 SNPs. Uniparental haplogroups were obtained from mitochondrial DNA and Y chromosome data (Supplementary table 1). New Guinean mitochondrial DNA diversity was previously described [START_REF] Pedro | Papuan mitochondrial genomes and the settlement of Sahul[END_REF], showing mostly haplogroups specific to the region (Q*:38%; P*:29%) as well as the presence of haplogroups frequent in Austronesian populations (B4a1a1:20%). New Guinean Y chromosome diversity is dominated by M1 (39%) and S* (18%) haplogroups frequent in the region, and haplogroups of Austronesian origin (O1:14%) as previously described [START_REF] Karafet | Major east-west division underlies Y chromosome stratification across Indonesia[END_REF][START_REF] Mona | Patterns of Y-chromosome diversity intersect with the Trans-New Guinea hypothesis[END_REF]. Autosomal data were combined with published genomes, notably from Indonesia including the western half of New Guinea Island) (Jacobs, et al. 2019), northeast Australia (Mallick, et al. 2016) and the Bismarck Archipelago [START_REF] Vernot | Excavating Neandertal and Denisovan DNA from the genomes of Melanesian individuals[END_REF], and with genotyping data from North Maluku [START_REF] Kusuma | The last sea nomads of the Indonesian archipelago: genomic origins and dispersal[END_REF] and imputed genotyping data from Papua New Guinea [START_REF] Bergstrom | A Neolithic expansion, but strong genetic structure, in the independent history of New Guinea[END_REF], resulting in a dataset of 966 individuals (including genomes from Denisovan, Neandertal and Chimpanzee, Supplementary table 2). By increasing the number of New Guinean whole genome sequencing data as well as the geographic resolution of New Guinean genetic diversity, we created a high quality dataset (e.g., with good phasing [START_REF] Delaneau | A linear complexity phasing method for thousands of genomes[END_REF]) and imputation [START_REF] Howie | A flexible and accurate genotype imputation method for the next generation of genome-wide association studies[END_REF])), thus allowing us to perform fine-scale demographic analyses.

Characterization of genetic structure

We analyzed genetic patterns in individuals from New Guinea, northeast Australia, the Bismarck Archipelago and Wallacea to identify key regional genetic groups. A coancestry matrix was built by scanning each pair of genomes to detect Identity-By-Descent (IBD) fragments using fineSTRUCTURE (Lawson, et al. 2012)(Supplementary Fig. 2). Most genomes clustered according to their respective geographical locations, defining 15 regional groups, including seven regional groups from the Wallacea/Sahul region: New Guinean southeast lowlands, New Guinean north and south lowlands, New Guinean east highlands, New Guinean west highlands, northeast Australia, the Bismarck Archipelago, and Wallacea. Of these seven regional groups of interest, we also characterized 58 subgroups based on the fineSTRUCTURE dendrogram (Fig. 1A and Supplementary table 2). Each of these 58 subgroups defines a 'population' based on genetic similarity (Lawson and Falush 2012), a definition that can differ from a linguistic or anthropological perspective but which is more accurate regarding the analyses performed in our work.

As previously described [START_REF] Bergstrom | A Neolithic expansion, but strong genetic structure, in the independent history of New Guinea[END_REF][START_REF] Hudjashov | Complex Patterns of Admixture across the Indonesian Archipelago[END_REF], recent gene flow from Asian groups, mostly linked to the Austronesian dispersal ~3kya, has substantially impacted regional genetic diversity. Since our study focused on Upper Pleistocene movements, we decided to mask this later Asian genetic ancestry using haplotypic information with PCAdmix [START_REF] Brisbin | Sabatier user on 16 August 2021 along each chromosome in individuals with admixed ancestry from two or more populations[END_REF]. Principal Component Analysis (PCA) shows that this procedure was efficient; prior to masking the Asian genetic ancestry, all genomes from New Guinea, Wallacea, northeast Australia and the Bismarck Archipelago were spread along the first principal component towards Asian genomes, while after masking, these genomes instead clustered together (Supplementary Fig. 3). The dataset was also analyzed using ADMIXTURE prior and after the masking of Asian tracks (Supplementary figure 4).

Asian genetic ancestry is present at an average of 22.7% in New Guinean lowlands, the Bismarck Archipelago, Wallacea and northeast Australia. It ranges between 0% and 62.1% in individuals from New Guinean lowlands, but absent in New Guinean highlanders, as previously reported [START_REF] Bergstrom | A Neolithic expansion, but strong genetic structure, in the independent history of New Guinea[END_REF], and between 0% and 49.2% in individuals from the Bismarck Archipelago, between 38.8% and 63.7% in individuals from Wallacea and between 8.1% and 12.1% in individuals northeast Australia, as previously reported [START_REF] Hudjashov | Complex Patterns of Admixture across the Indonesian Archipelago[END_REF]Mallick, et al. 2016) (Supplementary figure 4). After the masking, the average Asian genetic ancestry is greatly reduced (Supplementary figure 4), which allows us to infer demographic scenarios that are not strongly biased by recent Asian genetic admixture.

We first described the masked dataset using a PCA focused on individuals in the 7 regional groups of interest (Fig. 1B). All genomes from the same geographical regions are grouped together in distinctive clusters. Clusters of New Guinean lowlanders, Wallacean islanders and northeast Indigenous Australians fall relatively close together. New Guinean highlanders, especially individuals from the New Guinean east highlands, are separated on the second principal component. Individuals from the Bismarck Archipelago form two clusters separated by both main components. This apparent genetic structure coincides with F ST genetic distances (Supplementary Fig. 5). The genetic diversity of the Bismarck Archipelago is particularly differentiated from other New Guineans, Wallaceans and northeast Indigenous Australians (0.18<F ST <0.29; Supplementary Fig. 5). Within New Guinea, larger genetic distances are found between New Guinean lowlanders and highlanders (0.05<F ST <0.18; Supplementary Fig. 5) than between any other New Guinean sub-groups, as previously reported [START_REF] Bergstrom | A Neolithic expansion, but strong genetic structure, in the independent history of New Guinea[END_REF]. Wallaceans are genetically closer to New Guineans (0.04<F ST <0.26) than to Bismarck Islanders or northeast Indigenous Australians (0.06<F ST <0.318). We further confirm this geographical structure using EEMS [START_REF] Petkova | Visualizing spatial population structure with estimated effective migration surfaces[END_REF], which places the biggest genetic barriers between the Bismarck Archipelago and New Guinea, between the New Guinean lowlands and highlands, as previously noted [START_REF] Bergstrom | A Neolithic expansion, but strong genetic structure, in the independent history of New Guinea[END_REF], and between northeast Australia and New Guinea (Supplementary Fig. 6). We do not detect any strong genetic barriers in Wallacea.

All of these analyses clearly show that there is strong geographical structure of human genetic diversity in Sahul today. However, surprisingly, it also reveals a relative genetic continuity between New Guinea and Wallacea, suggesting that Wallacea might have represented less of a barrier to humans than other floral and faunal species.

Identification of gene flow

We analyzed the genetic diversity of our dataset to identify possible connections between each individual in the seven regional groups from New Guinea, northeast Australia, the Bismarck Archipelago and Wallacea using ADMIXTURE [START_REF] Alexander | Fast model-based estimation of ancestry in unrelated individuals[END_REF]). The genetic diversity in our dataset is best deconstructed into eight genetic components (Supplementary Fig. 789), four of which are dominated by New Guinean, northeast Indigenous Australian, Bismarck Archipelago and Wallacean groups (Fig. 1C and Supplementary Fig. 9). This strong genetic structure is confirmed by the genetic ancestries defined in New Guinean east and west highlanders as well as in the Bismarck Archipelago (red, orange and purple components respectively). The fourth genetic component (yellow) is more widespread, as it is present in individuals from Wallacea, the New Guinean lowlands, northeast Australia and the Bismarck Archipelago. It is the main component in Wallacea and New Guinea, and reaches its highest value in the southeast region (31-63%; Fig. 2A, Supplementary Fig. 7 and9).

To determine the history of regional mixing, population splits and admixture events were estimated without a priori assumptions with TREEMIX (Pickrell and Pritchard 2012)(Fig. 1D and Supplementary Fig. 10). An optimal model is obtained including five gene flows, one of which is in the Wallacea/Sahul region: from ancestors of northeast Indigenous Australians to a group ancestral to most Bismarck Archipelago islanders (Fig. 1D and Supplementary Fig. 10). This may partly explain the shared genetic ancestry observed in the ADMIXTURE plot (purple, Fig. 1C).

The limited number of admixture events is confirmed by f3-statistics performed between all pairs of the 58 sub-groups [START_REF] Patterson | Ancient admixture in human history[END_REF]. Only seven genetic admixture events are detected (Supplementary table 3), four of which show multiples waves of admixture. For all of these gene flows, we estimated the dates of the admixture events using MALDER [START_REF] Loh | Inferring admixture histories of human populations using linkage disequilibrium[END_REF]. Gene flows are detected between sub-groups of New Guinean west highlands, mostly between 9.9-11.9kya (Supplementary table 3). Multiple gene flow events between New Guinean east highlanders and New Guinean sub-groups from the south lowlands occurred between 25.4-30.2kya and 2.6-4.1kya. During the same two time frames, we also observe gene flows between a New Guinean sub-group from the south lowlands and a Wallacean sub-group from the Lesser Sunda Islands (Supplementary table 3). Although limited in number, these gene flows occurred at different time frames spanning several millennia, indicative of a population dynamics model involving multiple periods of isolation and contact.

Modeling the genetic scenario of the migration in Wallacea

Our analyses show that New Guinean groups appear to be more closely related to populations in Wallacea than to groups in northeast Australia and the Bismarck Archipelago (Fig. 1D and Supplementary Fig. 5). Since human populations settled Sahul from Sunda through Wallacea (Allen and O'Connell 2020), we expected that genetic diversity from New Guinea, northeast Australia and the Bismarck Archipelago would branch from the diversity of Wallacean sub-groups. This highest genetic affinity between New Guineans and Wallaceans could be due to ancient gene flow between these populations (Supplementary table 3), or to a genetic divergence between northeast Indigenous Australians and Wallaceans prior to the genetic divergence between New Guineans and Wallaceans (Fig. 1D). We therefore explored different scenarios that could have led to the observed genetic patterns.

We first sought to identify the most likely route taken by the first settlers of Sahul. We used an f3-outgroup approach to compare the genetic diversity present in New Guinean, northeast Indigenous Australian and Bismarck Archipelago islander sub-groups to each Wallacean sub-group, using African Yoruba as an outgroup (Supplementary Fig. 11). All subgroups from New Guinea, northeast Australia and the Bismarck Archipelago have f3statistics values equally distant to all Wallacean sub-groups located on the possible southern route to Sahul (Lesser Sunda Islands and Flores Island)(Supplementary Fig. 11A-C). However, when we included a sub-group located on the northern route (North Maluku islands, Kei Island), we detected a stronger genetic link with all studied sub-groups (Supplementary Fig. 11D-F). Most populations from New Guinea and the Bismarck Archipelago are significantly genetically closer to the population from North Maluku using f4-statistics (f4(Yoruba,X;Lesser Sunda Islands, North Maluku)>0.001, Z-score>3, Supplementary Fig. 12).

We found higher f3-statistics and more significant f4-statistics values for New Guinean subgroups with North Maluku than for sub-groups from the Bismarck Archipelago and northeast Australia (Fig. 2A and Supplementary Fig. 12). This suggests that the ancestors of New Guineans, Bismarck Archipelago islanders and northeast Indigenous Australians are all more closely related to North Maluku islanders than to groups from Lesser Sunda islands or Flores Island, and that this link is probably more important in New Guineans. This analysis supports two non-exclusive hypotheses: (i) the northern route as the most likely path for the initial settlement of Sahul, (ii) secondary gene flows between groups of North Maluku and New Guinea, the Bismarck Archipelago and northeast Australia.

We modeled these hypotheses using qpGraph [START_REF] Patterson | Ancient admixture in human history[END_REF]. For completeness, we tested several models that consider migrations from Sunda to Sahul along either the southern route (represented by the Lesser Sunda Islands) or the northern route (represented by the North Maluku Islands) (Supplementary Fig. 13), as well as different possibilities of admixture between Wallacean and New Guinean groups as suggested by our other analyses (Fig. 1C, Supplementary Fig. 6, Supplementary table 3). Since our analyses suggested that Wallaceans share a closer genetic link with New Guineans than northeast Indigenous Australians (Fig. 1D and 2A, Supplementary Fig. 6), we decided to model demographic scenarios with one or two migrations from Sunda to Sahul. Considering one migration from Sunda to Sahul, meaning that both New Guineans and Indigenous Australians share a common ancestor more recently than with Wallaceans, we obtained one model fitting the data (Z=-2.558, Figure 2B and Supplementary Fig. 13). It describes a common ancestry to both New Guineans and Indigenous Australians, after the split from Wallaceans, followed by major secondary gene flows from New Guinea to Wallacea (between 61-79%, Figure 2B and Supplementary Fig. 13). In this model no preferable route in Wallacea can be distinguished. We then tested models considering two migrations from Sunda to Sahul, in which northeast Indigenous Australians and New Guineans descended separately from groups moving along the northern or southern routes (Supplementary Fig. 14). We obtained one significant model: two migrations to Sahul, including one via the northern route, an admixture event between the ancestors of New Guineans and northeast Indigenous Australians, and a subsequent gene flow from ancestral groups from New Guinea to North Maluku (52%, Z = 2.337, Supplementary Fig. 14). In this model the ancestral group to New Guineans is related to the northern route, via North Maluku, but no route can be proposed for the ancestral Indigenous Australian group. Our qpGraph analyses thus support both models of migrations from Sunda to Sahul, with one or two migrations, followed by high secondary gene flows from New Guinea to Wallacea. We refined our models to include a group from the Bismarck Archipelago (represented by New Britain), considering it either to be a sister group to New Guineans or the result of a separate admixture event with ancestral northeast Indigenous Australians (Supplementary Fig. 15), as suggested by other analyses (Fig. 1D). For both models, with one or two migrations to Sahul, we obtained models fitting the data when the Bismarck Archipelago group results from a separate admixture event between ancestral New Guineans and northeast Indigenous Australians (One migration: Z = -2.570, Two migrations: Z = 2.301, Fig. 2B and Supplementary Fig. 15). This corresponds to the gene flow observed between an ancestral Indigenous Australian group and ancestral Bismarck Archipelago group detected by TREEMIX (Fig. 1D).

Modeling the genetic scenario of the settlement of Sahul

We added further detail to this scenario by reconstructing the effective population sizes of these groups, and estimating divergence dates between all studied groups using MSMC2 [START_REF] Schiffels | Inferring human population size and separation history from multiple genome sequences[END_REF]. Since archaic admixture can be a confounding factor in these analyses (Mallick, et al. 2016), we masked both Neandertal and Denisovan genomic tracks using ChromoPainter (Lawson, et al. 2012). The archaic genetic ancestry in the newly sequenced genomes ranges between 1.8% and 2.8% of Denisovan ancestry and ranges between 2.3% and 2.5% of Neandertal ancestry, corresponding to previous estimations (Jacobs, et al. 2019;[START_REF] Vernot | Excavating Neandertal and Denisovan DNA from the genomes of Melanesian individuals[END_REF]). All estimated effective population size curves show a bottleneck around 60-70 kya, as found for European and Asian populations (Supplementary Fig. 16). Around 50kya, estimated effective population sizes for populations from Wallacea and Sahul are between 1671-2100 individuals, and between 2540-2999 individuals in Eurasia and 9506 individuals in Africa. All divergence dates for Wallacea/Sahul groups from Africa, represented by Yoruba genomes, show similar results (74 ± 4.2kya, Supplementary Fig. 17). We note here that this date is older than that obtained between Eurasian and African genomes (66.5±3.7kya), a result previously reported and interpreted as a potential methodological bias or the signal of an even earlier human migration from Africa [START_REF] Bergstrom | A Neolithic expansion, but strong genetic structure, in the independent history of New Guinea[END_REF][START_REF] Malaspinas | A genomic history of Aboriginal Australia[END_REF]Mallick, et al. 2016;[START_REF] Pagani | analyses inform on migration events during the peopling of Eurasia[END_REF]). We cannot confirm any of these interpretations but the similar dates of cross-coalescence with African genomes and the similar profiles of effective population size curves clearly show that all genomes from Wallacea, New Guinea, northeast Australia and the Bismarck Archipelago descend from a common ancestor, probably located in Sunda, as represented by our TREEMIX analysis (Fig. 1D). We then estimated divergence dates between each of the studied groups (Fig. 2C). Genomes from New Guinean and Wallacean groups diverged between 23.7-26.7kya. Groups from New Guinea and the Bismarck Archipelago show genetic divergence between 17.5-19.8kya.

Genomes of New Guineans diverged from the genomes of northeast Indigenous Australians between 13.5-15.3kya. We found that genomes from northeast Indigenous Australians diverged from the genomes of Wallaceans 13.3-15.2kya, and from genomes of Bismarck Archipelago islanders between 15.2-17.3kya. We caution that these analyses are sensitive to recent admixture events, and the masking of Asian genomic tracks reduces the amount of data to accurately estimate divergence dates. It likely affected results on Wallacean genomes for which the masked Asian genetic data were large (>50%). Moreover, the limited number of Indigenous Australian genomes likely resulted in an insufficient phasing quality for the two available genomes which could strongly affect MSMC2 analyses, as previously noted [START_REF] Malaspinas | A genomic history of Aboriginal Australia[END_REF]. A previous study on Indigenous Australians genomes [START_REF] Malaspinas | A genomic history of Aboriginal Australia[END_REF](not available for this study) found that all Indigenous Australians groups diverged from each other around 25kya, according to MSMC2 analyses, and from New Guinean highlanders around 30-40kya. Our estimation of the divergence between Indigenous Australians and New Guinean is even more recent, potentially caused by the identified gene flow between groups of northeast Australia and the Bismarck Archipelago (Fig. 1D). Although these dates cannot be interpreted as being part of the initial migration to Sahul, they strongly indicate long-term contacts between all of these groups, as seen in our models. The population dynamics between ancestral groups was probably important, over long distances and periods of time, leading to the emergence of modern populations like New Guineans.

Modeling the genetic scenario of the settlement of New Guinea

Our TREEMIX analysis show that New Guinean sub-groups branch together. The root of the tree lies in the genetic diversity present today in New Guinean southeast lowlanders (Fig. 1D). This result reflects the ADMIXTURE analysis, which shows that the main genetic component in Wallacea is also the main component in southeast New Guinea, reinforcing the hypothesis of an ancient link (yellow component, Fig. 1C and Supplementary Fig. 9). Subsequently, sub-groups from the New Guinean north and south lowlands form two different genetic clusters that are genetically drifted from the New Guinean southeast Downloaded from https://academic.oup.com/mbe/advance-article/doi/10.1093/molbev/msab238/6349182 by Univeristé Paul Sabatier user on 16 August 2021 lowlanders (Fig. 1D). This suggests that both clusters emerged from two separate migrations, south and north, initiated from the southeast New Guinean region. We modeled the settlement of New Guinea based on our results showing that most of New Guinean genetic diversity could be described by genetic drift with limited occurrences of gene flows (Fig. 1D and Supplementary table 3). We used qpGraph to determine the best scenario fitting data from the five New Guinean groups (southeast lowlanders, south lowlanders, north lowlanders, east highlanders and west highlanders) and one African outgroup (Yoruba). First we modeled the most probable entry point into New Guinea, considering possible migrations from the north lowlands, the south lowlands and the southeast lowlands, with or without admixture events (Supplementary Fig. 18). One model is significant (Z-score = -0.06); this places southeast lowlands as the most likely entry point but only if there is a secondary admixture event with ancestral New Guinean north lowlanders (Supplementary Fig. 18). This model corresponds with results from other analyses, which positions the New Guinean southeast group at the root of the New Guinean diversity (Fig. 1D).

Highlanders are the last to branch off in our TREEMIX analysis (Fig. 1D). Our analyses showed that New Guinean west and east highlanders form two clearly separated clusters (Fig. 1A and 1B, Supplementary Fig. 2 and5), and that both New Guinean highlander groups drifted from New Guinean lowlanders located in the southern region (Fig. 1D). We first added a New Guinean east highlander group to model dispersal in the highlands from the north, south or southeast lowlands; using qpGraph (Supplementary fig. 19). The model best fitting the data connects New Guinean east highlanders to a branch leading to the New Guinean north lowlanders (Z-score=-0.068, Supplementary Fig. 19). Finally we added a New Guinean west highlander group, defined as a sister group to New Guinean east highlanders or coming from a second migration from the lowlands (Supplementary Fig. 20). The best model places New Guinean west highlanders with New Guinean south lowlanders (Z-score=2.759, Supplementary Fig. 20). This branching describes two separate migrations leading to the genetic diversity observed in New Guinean east and west highlanders, which would explain the stark genetic difference observed in other analyses (Fig. 1B and1C). However we note here that it could appear to be in contrast to the TREEMIX analysis which showed the New Guinean east highlanders branched with some PNG South lowlands groups (Fig. 1D). In particular these South lowlands groups (PNG South lowlands subgroups 39 and 40, Fig. 1D) are geographically located close to North lowlands groups. We computed an f3-outgroup analysis to further explore the connection between PNG highlanders and lowlanders. It confirmed the closer link between PNG West highlanders and PNG South lowlanders but also showed that PNG East highlanders are almost genetically equidistant to both PNG North and South lowlanders, although still more closely related to PNG South highlanders (Supplementary figure 21). Based on our results, we postulate that PNG east highlanders came from a place intermediate to both South and North lowlanders, potentially in the current Morobe province. The qpGraph model also supports a secondary admixture event between New Guinean west highlanders and north lowlanders (Fig. 3A and Supplementary figure 20), which could explain the presence of the main New Guinean west highlanders genetic component in many lowlanders (Fig. 1C) as well as previous reports on the link between PNG north lowlanders and PNG highlanders [START_REF] Bergstrom | A Neolithic expansion, but strong genetic structure, in the independent history of New Guinea[END_REF]). This final model identifies the southeast lowlands as the entry point into New Guinea, followed by two migrations to the south and north lowlands, which in turn led to the independent settlements of the east and west highlands.

Given this scenario, we estimated divergence dates between each New Guinean group of interest using MSMC2 (Fig. 3B). Divergence dates between New Guinean lowlanders and New Guinean highlanders ranged between 8.8-12.8kya, which broadly correspond to previous estimations (12-16kya) [START_REF] Bergstrom | A Neolithic expansion, but strong genetic structure, in the independent history of New Guinea[END_REF]. West and east New Guinean highlander groups diverged between 7.6-8.8kya, similarly to previously reported (6-12kya) [START_REF] Bergstrom | A Neolithic expansion, but strong genetic structure, in the independent history of New Guinea[END_REF]). These rather recent dates of genetic divergence probably reflect long-term contacts between New Guinean groups following the initial settlement of the island, initiated as early as 49kya based on archaeological records [START_REF] Summerhayes | Human adaptation and plant use in highland New Guinea 49,000 to 44,000 years ago[END_REF]. They thus probably illustrate the gene flows between highlander groups, and between highlanders and lowlanders, which are estimated to have occurred at the same time (9.9-11.9kya, Supplementary table 3). Effective population sizes estimated with MSMC2 for each New Guinean group do not show any marked differences (Supplementary fig. 16). It was previously interpreted as a bias due to phasing quality [START_REF] Bergstrom | A Neolithic expansion, but strong genetic structure, in the independent history of New Guinea[END_REF]), but it could also mean that contacts between New Guinean highlanders and lowlanders were frequent enough so that populations grew similarly.

Discussion

Downloaded from https://academic.oup.com/mbe/advance-article/doi/10.1093/molbev/msab238/6349182 by Univeristé Paul Sabatier user on 16 August 2021

Our study provides the first models of the settlement of north Sahul and New Guinea based on whole genome data (Fig. 4). Our models support an early split of the ancestral groups of New Guineans from northeast Indigenous Australians, potentially in Wallacea. Among the possible migratory routes, our model brings new evidence for the use of the northern route, via the North Maluku islands. Our genetic data describe a settlement of New Guinea initiated from what is now its southeast lowland via the Arafura Plain, and long-term population contact between ancestral populations of New Guinea, northeast Australia, Wallacea and the Bismarck Archipelago.

Multiple migrations between Wallacea and Sahul

To reach Sahul, the first human settlers had to cross Wallacea, a region always composed of small islands, large water crossings and strong ecological barriers (Wallace 1869). These restrictions on movement are reflected in the substantial genetic structuring of current populations (Fig. 1 and Supplementary Fig. 5). But our results also show relative genetic continuity between groups from Wallacea to New Guinea, with stronger genetic barriers between New Guinea and Australia, as well as around the Bismarck Archipelago (Supplementary Fig. 6), as described in a previous report [START_REF] Bergstrom | A Neolithic expansion, but strong genetic structure, in the independent history of New Guinea[END_REF]. Wallacea was probably less a strong geographical frontier than a place of possible migrations to Sahul [START_REF] Allen | A different paradigm for the initial colonisation of Sahul[END_REF][START_REF] Birdsell | The recalibration of a paradigm for the first peopling of Greater Australia. Sunda and Sahul: Prehistoric Studies in Southeast Asia[END_REF][START_REF] Kealy | Reconstructing Palaeogeography and Inter-island Visibility in the Wallacean Archipelago During the Likely Period of Sahul Colonization, 65-45000Years Ago[END_REF]. The ancestral population to New Guineans, Bismarck Archipelago Islanders and northeast Indigenous Australians might have been composed of around 1300 individuals based on ecological information (Bradshaw, et al. 2019) and 1671-2100 individuals based on genetic data (Supplementary Fig. 16).

Our models, like those of others [START_REF] Malaspinas | A genomic history of Aboriginal Australia[END_REF][START_REF] Teixeira | Using hominin introgression to trace modern human dispersals[END_REF], show that this ancestral population split into different groups early during this migration, upon the arrival in Sahul 50kya (model A, fig. 3B and4), or even before in Wallacea, from which multiple migrations could have been conducted within a small time frame (i.e. in less than 900 years (Bradshaw, et al. 2019)) (model B, Fig. 3B and4). A model with a single migration is still valid given our genetic data but there is growing non-genetic evidence that the initial settlement of Sahul might have been made by multiple migrations (Bird, et al. 2019;Bradshaw, et al. 2021;Bradshaw, et al. 2019;[START_REF] Clarkson | Human occupation of northern Australia by 65,000 years ago[END_REF][START_REF] Crabtree | Landscape rules predict optimal superhighways for the first peopling of Sahul[END_REF]. This hypothesis has never been clearly explored using genomic data. Our study offers the first opportunity to test this model, and we found that our genomic data support a model with two migrations to Sahul to explain the current genetic diversity in New Guineans, Bismarck Archipelago Islanders and Indigenous Australians (model B, Fig. 3B and4). Among the possible routes from Wallacea to Sahul, this model describes the northern route through North Maluku as a preferable path taken by the ancestors of New Guineans (Fig. 2A). It is consistent with studies on paleogeological reconstructions of coastlines (Fig. 4) (Bird, et al. 2019;[START_REF] Kealy | Reconstructing Palaeogeography and Inter-island Visibility in the Wallacean Archipelago During the Likely Period of Sahul Colonization, 65-45000Years Ago[END_REF] as well as studies using ecological models (Bradshaw, et al. 2021).

We cannot propose a location for the second route given the limited number of available Indigenous Australian genomes, but it may have been different, and potentially used at different times during the Upper Pleistocene, as was recently described (i.e. possible entry though the northwest Sahul Shelf around 75kya (Bradshaw, et al. 2021)). Our genomic model B also agrees with previous studies on uniparental markers that propose stark differences in genetic diversity between New Guineans and northeast Indigenous Australians as a sign of multiple migrations [START_REF] Pedro | Papuan mitochondrial genomes and the settlement of Sahul[END_REF][START_REF] Redd | Peopling of Sahul: mtDNA variation in aboriginal Australian and Papua New Guinean populations[END_REF]. The genomic model B thus brings together a consensus view that New Guineans and northeast Indigenous Australians share a common ancestor, in Sunda, [START_REF] Malaspinas | A genomic history of Aboriginal Australia[END_REF][START_REF] Teixeira | Using hominin introgression to trace modern human dispersals[END_REF], therefore integrating recent perspectives on multiple migrations to Sahul [START_REF] Allen | A different paradigm for the initial colonisation of Sahul[END_REF][START_REF] Redd | Peopling of Sahul: mtDNA variation in aboriginal Australian and Papua New Guinean populations[END_REF]. Future work on modern and ancient DNA in Wallacea, New Guinea and Australia will surely help to disentangle the two models presented here.

Wallacea was clearly not a major barrier to human migrations (Supplementary Fig. 6) but it was also not necessarily a simple highway to Sahul (Fig. 4). We hypothesize that the multitude of Wallacean islands might have driven genetic differentiation between related groups. These human groups would have progressively migrated east, in a planned effort spanning over a millennium or more (Bradshaw, et al. 2019). This time was probably necessary for the populations to grow large enough before migrating further east (Bradshaw, et al. 2021), leading to the emergence of genetic differentiation. Several human populations could have crossed Wallacea to land on Sahul, using different routes (Bradshaw, et al. 2021).

During this very early period of the settlement of Sahul, it is likely that the demographic successes of these groups were different, as some archaeological sites might indicate [START_REF] Clarkson | Human occupation of northern Australia by 65,000 years ago[END_REF]. Admixture between ancestral populations, which countered the potential effects of inbreeding and drift in isolated groups, might have been an important Downloaded from https://academic.oup.com/mbe/advance-article/doi/10.1093/molbev/msab238/6349182 by Univeristé Paul Sabatier user on 16 August 2021 factor for the dispersal success of certain human groups in Sahul. A previous study found relatively late divergence dates between New Guinean highlanders and Indigenous Australians (30-40kya), as well as between Indigenous Australians groups (25kya) [START_REF] Malaspinas | A genomic history of Aboriginal Australia[END_REF]. It also found signals of post-divergence contacts between ancestral groups of New Guineans and Indigenous Australians in the north of Australia [START_REF] Malaspinas | A genomic history of Aboriginal Australia[END_REF], which agrees with the late genomic divergence dates that we obtain between southeast New Guinean lowlanders and northeast Indigenous Australians (13.5-15.3kya, Fig. 2C). These dates need to be interpreted with caution, because of technical issues raised above and by others [START_REF] Malaspinas | A genomic history of Aboriginal Australia[END_REF]), but they do show that Sahul, and especially the north region, was probably a place of interactions between human groups with related, but already differentiated, genomes. These interactions probably involved few individuals given the number of region-specific mitochondrial DNA [START_REF] Pedro | Papuan mitochondrial genomes and the settlement of Sahul[END_REF]) and Y chromosome lineages (Bergstrom, et al. 2016) shared between New Guineans and Indigenous Australians. This led to the emergence of the unique genetic diversities passed on to the modern populations of northeast Australia, the Bismarck Archipelago and New Guinea.

The complex settlement of New Guinea

New Guinea was settled rapidly after the initial migrations to Sahul as revealed by archaeological records [START_REF] Summerhayes | Human adaptation and plant use in highland New Guinea 49,000 to 44,000 years ago[END_REF]. This fast dispersal has made it challenging to determine the route taken to reach this territory. Our study strongly supports an entry point located in what is now the New Guinean southeast region (Fig. 4). The genetic distances between current New Guinean southeast lowlanders and Wallaceans are smaller than with any other New Guineans (Supplementary Fig. 5). The main genetic component in Wallaceans is also mostly present in groups from the New Guinean southeast lowlands (Fig. 1C and Supplementary Fig. 9). Among the different models tested, the only statistically significant scenario places New Guinean southeast lowlanders at the root of New Guinean genetic diversity, even considering secondary gene flows (Fig. 3A). This remarkably corresponds to the location of the oldest known archaeological sites, in the Ivané valley, dated at 45-49kya [START_REF] Summerhayes | Human adaptation and plant use in highland New Guinea 49,000 to 44,000 years ago[END_REF]). Our genetic model located the, apparently ancient, genetic background to the far east of New Guinea bringing new support for an initial migration via the ancient plain currently under the Arafura Sea [START_REF] Allen | A different paradigm for the initial colonisation of Sahul[END_REF]. Recent studies using ecological modeling identified the most likely route to New Guinea in the ancient plain currently under the Arafura Sea along the southeast coast of New Guinea (Bradshaw, et al. 2021;[START_REF] Crabtree | Landscape rules predict optimal superhighways for the first peopling of Sahul[END_REF], converging with our genetic model. Other hypothetical routes through the New Guinean north coast, the highlands or the south foothills could have been taken by other groups of settlers [START_REF] Summerhayes | The archaeology of forest exploitation and change in the tropics during the Pleistocene: The case of Northern Sahul (Pleistocene New Guinea)[END_REF]) but their demographic success, based on genetic data, must have been limited. During the Upper Pleistocene, the Arafura Plain south of modern New Guinea was dominated by savannah and lakes [START_REF] Summerhayes | The archaeology of forest exploitation and change in the tropics during the Pleistocene: The case of Northern Sahul (Pleistocene New Guinea)[END_REF]). It would have been a favorable ecosystem for human dispersal from the west coast to the east coast of Sahul, before reaching what is now the southeast coast of the island of New Guinea (Fig. 4).

Settlement within New Guinea was initiated from the lowlands. We found that the genetic diversity of groups from the New Guinean north and south lowlands derived separately from groups in the southeast region (Fig. 1D and3A). During the first stages of the settlement of New Guinea, the lowlands were probably the most hospitable places due to their access to littoral resources [START_REF] Summerhayes | The archaeology of forest exploitation and change in the tropics during the Pleistocene: The case of Northern Sahul (Pleistocene New Guinea)[END_REF]). However, the dry conditions and steep geography might have created locations unfavorable to permanent human settlements, potentially numerous in the highlands and the north lowlands [START_REF] Terrell | The'sleeping giant'hypothesis and New Guinea's place in the prehistory of Greater Near Oceania[END_REF]. Around 35kya, these adverse conditions appeared to have significantly reduced human mobility in New Guinea [START_REF] Allen | A different paradigm for the initial colonisation of Sahul[END_REF]. This would explain the relative slow population growth we observe at this period for all New Guinean groups, in comparison to European or African groups (Supplementary Fig. 16). A similar explanation was put forward to interpret the late coalescence dates of many New Guinean maternal lineages, postdating 30kya [START_REF] Pedro | Papuan mitochondrial genomes and the settlement of Sahul[END_REF]. This lack of resources could have induced New Guinean populations to explore other locations, such as the highlands [START_REF] Summerhayes | The archaeology of forest exploitation and change in the tropics during the Pleistocene: The case of Northern Sahul (Pleistocene New Guinea)[END_REF].

Our analyses show that the New Guinean highlands were settled by two separate migratory waves: one migration from the east foothills to the east highlands, and another migration from the south lowlands to the west highlands (Fig. 1D, 3A and 4). These two migrations could have used two major geological corridors defined by the Markham and Ramu Rivers in the east, and the Fly River in the south [START_REF] Summerhayes | The archaeology of forest exploitation and change in the tropics during the Pleistocene: The case of Northern Sahul (Pleistocene New Guinea)[END_REF]). Although we do not find support for the Sepik River as a major migratory route for the initial settlement of the highlands [START_REF] Bergstrom | A Neolithic expansion, but strong genetic structure, in the independent history of New Guinea[END_REF], all these natural corridors were likely a major path for secondary gene flows between New Guinean highlanders and New Guinean lowlanders, as indicated by our analyses (Fig. 3A). Many gene flows appear to have occurred during the early period of the settlement of the highlands, around 25.4-30.2kya, suggesting continuous contacts between highlanders and lowlanders (Fig. 3A and Supplementary table 3). This is consistent with archaeological data arguing for long-term, but not permanent, settlements in the highlands around 25kya [START_REF] Summerhayes | The archaeology of forest exploitation and change in the tropics during the Pleistocene: The case of Northern Sahul (Pleistocene New Guinea)[END_REF]. New Guinean populations probably lived in semi-permanent settlements, gathering resources in the highlands and the foothills, but keeping contacts with lowlander groups. These populations would have progressively fixed their settlement in the highlands, reducing contacts with the lowlanders, since no gene flow was detected between 20kya and 3kya, and increasing the genetic differentiation between highlanders and lowlanders (Fig. 1 and 3, Supplementary Fig. 5), as previously reported [START_REF] Bergstrom | A Neolithic expansion, but strong genetic structure, in the independent history of New Guinea[END_REF].

However an intense period of population dynamics occurred locally in the New Guinean highlands around 10kya. We detected several admixture events between west highlander groups at this time (Supplementary table 3). This period was marked by the development of agricultural practices and the emergence of a Neolithic culture, especially in the Waghi valley of the west highlands [START_REF] Denham | Origins of agriculture at Kuk Swamp in the highlands of New Guinea[END_REF][START_REF] Shaw | Emergence of a Neolithic in highland New Guinea by 5000 to 4000 years ago[END_REF]. As previously reported [START_REF] Bergstrom | A Neolithic expansion, but strong genetic structure, in the independent history of New Guinea[END_REF]), this was a major cultural change that likely influenced gene flows in the highlands and the dispersal of cultural traits, such as the Trans-New Guinean linguistic family [START_REF] Reesink | Explaining the linguistic diversity of Sahul using population models[END_REF]). Lowlanders might not have been significantly influenced by the development of a Neolithic culture in the highlands [START_REF] Pawley | Papuan pasts: Cultural, linguistic and biological histories of Papuan-speaking peoples: Pacific Linguistics, Research School of Pacific and Asian Studies, The[END_REF], but their population size was not necessarily much lower than those found in the highlands (Supplementary fig. 16) thanks to exchanges of goods and genes between the two regions (Supplementary table 3).

A second wave of gene flows between highlanders and lowlanders is only detected around 2.6-4.1kya (Supplementary table 3). At the same period, gene flow is also detected between New Guinean lowlanders and Wallaceans (Supplementary table 3). These events likely correspond to the influence of the Austronesian culture 3kya [START_REF] Lipson | Reconstructing Austronesian population history in Island Southeast Asia[END_REF], but rather than initiating new contacts between New Guinean groups, the arrival of Austronesians seems to have only reignited New Guinean population dynamics that already existed around 25kya (Supplementary table 3). Genes, material culture and ideas were already moving within New Guinea [START_REF] Pawley | Papuan pasts: Cultural, linguistic and biological histories of Papuan-speaking peoples: Pacific Linguistics, Research School of Pacific and Asian Studies, The[END_REF][START_REF] Torrence | Social networks and the spread of Lapita[END_REF]. The last stage of the settlement of New Guinea was probably another pulse in a network of population interactions built up over millennia.

  Downloaded from https://academic.oup.com/mbe/advance-article/doi/10.1093/molbev/msab238/6349182 by Univeristé Paul Sabatier user on 16 August 2021

  Downloaded from https://academic.oup.com/mbe/advance-article/doi/10.1093/molbev/msab238/6349182 by Univeristé Paul Sabatier user on 16 August 2021

Figures captions Figure 1 .

 captions1 Figures captions Figure 1. Genetic diversity of individuals from New Guinea, Wallacea, the Bismarck Archipelago and northeast Australia. (A) Geographic locations of all 58 genetic sub-groups defined by the fineSTRUCTURE analysis. (B) Principal Component Analysis representing 26.7% of the total variability of all individual genomes from the 58 sub-groups, masked for

Figure

  Figure 3

  Figure 3

Figure 4

 4 Figure 4

  Downloaded from https://academic.oup.com/mbe/advance-article/doi/10.1093/molbev/msab238/6349182 by Univeristé Paul Sabatier user on 16 August 2021

Downloaded from https://academic.oup.com/mbe/advance-article/doi/10.1093/molbev/msab238/6349182 by Univeristé Paul Sabatier user on 16 August 2021

Acknowledgments

We thank Georgia Kaipu (National Research Institute, Papua New Guinea), Christopher Kinipi (University of Papua New Guinea Clinic), Alois Kuaso and Kenneth Miamba (National Museum and Art Gallery, Papua New Guinea) for their help during the sampling campaigns.

We thank Chris Clarkson for his comments on the manuscript. We thank J. Stephen Lansing (Santa Fe Institute) for financially supporting some of the sequencing. We thank Monika Karmin for her help on the analysis of Y chromosome data. We thank Opale Coutant (University of Toulouse) and Adeline Morez (Liverpool University) for their help in data analysis. We acknowledge support from the GenoToul bioinformatics facility of Genopole Toulouse Midi-Pyrénées, France; the LabEx TULIP, France; and the Centre National de Recherche en Génomique Humaine (CNRGH), Evry, France. We especially thank all of our study participants. We acknowledge support from the National Geographic Society (Grant HJ-156R-17) and the Leakey Foundation. This research was supported by the European Union through Horizon 2020 research and innovation programme under grant no 810645 and through the European Regional Development Fund project no. MOBEC008. This work was supported by the French Ministry of Research grant Agence Nationale de la Recherche (ANR PAPUAEVOL), the French Ministry of Foreign and European Affairs (French Prehistoric Mission in Papua New Guinea), and the French Embassy in Papua New Guinea, as well as a fellowship from the Alexander von Humboldt Foundation.

Data and Code Availability Statement

Downloaded from https://academic.oup.com/mbe/advance-article/doi/10.1093/molbev/msab238/6349182 by Univeristé Paul Sabatier user on 16 August 2021

New Guinean genomes reveal possible scenarios describing the migrations from Sunda to Sahul and a settlement of New Guinea initiated in the southeast lowlands, followed by two separate dispersals to the highlands (Fig. 4). Our models will necessarily be confronted by future genetic and archaeological data, particularly in the western half of New Guinea Island and Australia. Paleogenomics in Oceania is an emerging field and it has already proven to be highly informative about the history of populations with New Guinean genetic ancestry [START_REF] Choin | Genomic insights into population history and biological adaptation in Oceania[END_REF]Skoglund, et al. 2016). Future research will certainly bring key information on the complex population dynamics that existed during the past 50,000 years in north Sahul.

Materials and Methods

Sampling

Fifty eight individuals from Papua New Guinea participated in the study (Supplementary All samples were collected from healthy unrelated adult donors, all of whom provided written informed consent (Supplementary text). DNA was extracted from saliva samples with the Oragene sampling kit according to the manufacturer's instructions. After a full presentation of the project to a wide audience, a discussion with each individual willing to participate ensured that the project was fully understood. Once participants had signed informed-consent forms, we interviewed them to obtain information on their date and place of birth, their spoken language(s), and similar data related to their genealogy (up to 2nd or 3rd generation) to establish local ancestry. Based on the participants' genealogy, we chose the samples to be sequenced so as to cover most of Papua New Guinean genomic diversity. (Mallick, et al. 2016); 60 Papuan genomes (EGA: EGAS00001001247 dbGAP: phs001085.v1.p1) [START_REF] Malaspinas | A genomic history of Aboriginal Australia[END_REF][START_REF] Vernot | Excavating Neandertal and Denisovan DNA from the genomes of Melanesian individuals[END_REF]) and 161 genomes from Indonesia (EGA: EGAS00001003054) (Jacobs, et al. 2019). SNP calling was performed on the combined dataset. Reads were aligned to the 'decoy' version of the GRCh37 human reference sequence (hs37d5) using BWA MEM [START_REF] Li | Toward better understanding of artifacts in variant calling from high-coverage samples[END_REF]. Duplicated reads were removed with picard-tools v.2.12.0 (http://broadinstitute.github.io/picard) and local realignment around indels were performed with GATK v. 3.5 [START_REF] Poplin | Scaling accurate genetic variant discovery to tens of thousands of samples[END_REF]). After alignment, the average sequencing depth across the samples was 36x, with a minimum of 31x coverage. Base-calling and per-sample gVCF files were generated using GATK HaplotypeCaller (reads mapping quality ≥ 20) [START_REF] Poplin | Scaling accurate genetic variant discovery to tens of thousands of samples[END_REF]. Multisample VCF files were obtained with CombineGVCFs and GenotypeGVCFs.

Datasets

Genotype calls were subsequently filtered for base depth (≥ 8x and ≤ 400x) and genotyping quality (≥ 30) using BCFtools v. 1.4 [START_REF] Li | Toward better understanding of artifacts in variant calling from high-coverage samples[END_REF]. We removed sites within segmental 2). All whole genomes were phased using SHAPEIT v. 2.r79033 [START_REF] Delaneau | A linear complexity phasing method for thousands of genomes[END_REF]) with default settings, except for iterations (50) and states ( 200) and used the HapMap phase 2 genetic map (International HapMap, et al. 2007) without a reference panel.

In addition to the whole genome dataset, we compiled another dataset including genotyping data from Papua New Guinea [START_REF] Bergstrom | A Neolithic expansion, but strong genetic structure, in the independent history of New Guinea[END_REF]. This genotyping data was phased and imputed using all phased New Guinean whole genomes as a reference panel (N=166) with IMPUTE2 v. 2.3.2 [START_REF] Howie | A flexible and accurate genotype imputation method for the next generation of genome-wide association studies[END_REF]. Ten Papuan WGS randomly selected were downgraded to the same set of SNPs as the genotyping data to serve as an internal control. All SNPs with an imputation score above 0.4 were kept, leaving 12,168,294 SNPs. The 10 downgraded genomes are over 98% consistent with actual genotypes from whole genomes for the same individuals. We also added genotyping data from North Maluku [START_REF] Kusuma | The last sea nomads of the Indonesian archipelago: genomic origins and dispersal[END_REF]), but we did not impute missing SNPs since no whole genomes are available for this group. Finally we included data from Neandertal [START_REF] Prufer | The complete genome sequence of a Neanderthal from the Altai Mountains[END_REF], Denisovan [START_REF] Reich | Genetic history of an archaic hominin group from Denisova Cave in Siberia[END_REF]) and Chimpanzee [START_REF] Mikkelsen | Initial sequence of the chimpanzee genome and comparison with the human genome[END_REF]. The final dataset of 966 genomes was phased with SHAPEIT v. 2.r79033 [START_REF] Delaneau | A linear complexity phasing method for thousands of genomes[END_REF]) without a reference panel with the same settings as above.

To focus on Upper Pleistocene genetic history, recent Asian ancestry was masked in admixed individuals from New Guinea, Wallacea, northeast Australia and the Bismarck Archipelago using PCAdmix v.1.0 [START_REF] Brisbin | Sabatier user on 16 August 2021 along each chromosome in individuals with admixed ancestry from two or more populations[END_REF]. Since both available Indigenous Australian genomes showed signs of recent European admixture (Mallick, et al. 2016), we also masked European genetic ancestry for all individuals. Two parental metapopulations were defined by randomly choosing 100 Asian-European individuals and 100 New Guinean individuals. For the latter group, only unadmixed New Guinean whole genome sequences were selected, as defined by ADMIXTURE v. 1.3 [START_REF] Alexander | Fast model-based estimation of ancestry in unrelated individuals[END_REF]). Local ancestry was estimated for all individuals from New Guinea, Wallacea, northeast Australia and the Bismarck Archipelago showing non-Eurasian genetic ancestry above 20%. Posterior probabilities for New Guinean ancestry above 0.9 and Asian ancestry below 0.1 were criteria to define SNPs with New Guinean ancestry. We applied quality controls using Plink v. 1.932 [START_REF] Chang | Second-generation PLINK: rising to the challenge of larger and richer datasets[END_REF] to filter for and exclude (1) close relatives by using an identity-by-descent (IBD) estimation with an upper threshold of 0.25 (second-degree relatives); (2) SNPs that failed the Hardy-Weinberg exact test (p < 10 -6 ) within each group; (3) samples with a call rate < 0.99 and displaying missing rates > 0.05 across all samples in each population. For analyses requiring linkage disequilibrium pruning, we used Plink to filter for variants in high linkage disequilibrium (r2 > 0.2), leaving 74,067 SNPs. Some Wallacean groups were pooled together because of the high rate of masked SNPs, according to the fineSTRUCTURE analysis and geography, in order to perform population-based analyses such as TREEMIX, f3 (ex:

Wallacea_3-4_6-7 refers to the pool of data from the subgroups Wallacea_1, Wallacea_2, Wallacea_6, Wallacea_7). Principal Components Analyses were run with smartpca from the EIGENSOFT v.7.2.1 package [START_REF] Patterson | Population structure and eigenanalysis[END_REF]. Populations were genetically defined using fineSTRUCTURE v. 2.07 (Lawson, et al. 2012). This package performs model-based Bayesian clustering of genotypes based on the shared IBD fragments between each pair of individuals, without self-copying, calculated with ChromoPainter v. 2.0 (Lawson, et al. 2012). Five million MCMC runs were performed. A coancestry heatmap and a dendrogram were inferred to visualize the number of statistically defined clusters that describe the data (Supplementary Fig. 2). 15 genetic groups were defined, and subdivided into 106 sub-groups, including 58 sub-groups from

Wallacea and Sahul (Supplementary Fig. 2 and Supplementary table 2). Genotyping data from North Maluku were defined as an additional subgroup (Wallacea_11) since the SNP density is much lower than whole genome sequencing data and no imputation could be performed.

Genetic barriers and corridors were defined using EEMS v. 1 [START_REF] Petkova | Visualizing spatial population structure with estimated effective migration surfaces[END_REF]). We used approximate geographic coordinates, a genetic dissimilarity matrix between populations and a map of New Guinea and the surrounding areas defining a grid of 700 modeled demes.

Depending on their location, several sub-groups can be included in one deme, whose size increased accordingly. We ran 5 million MCMC iterations before checking for convergence of the MCMC chain (Supplementary Fig. 6). Plots were generated in R v. 3.4.3 according to instructions in the EEMS manual. F ST genetic distances were obtained with the EIGENSOFT v.7.2.1 package [START_REF] Patterson | Population structure and eigenanalysis[END_REF].

Genetic ancestries present in the dataset were estimated by ADMIXTURE v. 1.3 [START_REF] Alexander | Fast model-based estimation of ancestry in unrelated individuals[END_REF]), with default settings, for components K = 2 to K = 15. Ten iterations with randomized seeds were run and compiled with CLUMPAK v. 1 [START_REF] Kopelman | Clumpak: a program for identifying clustering modes and packaging population structure inferences across K[END_REF]. We used the minimum average cross-validation value to define the most informative K components (here, K = 8), and the major modes defined by CLUMPAK v. 1 were reported.

Geographical distributions of the main genetic components found in Wallacea and Sahul were obtained using inverse distance weighting interpolation between data points with QGIS v. 3.16.1 (GDAL-SOFTWARE-SUITE 2013).

f3-statistics and f3-outgroup were computed with qp3pop from the ADMIXTOOLS v.6.0 package [START_REF] Patterson | Ancient admixture in human history[END_REF]) between all trios of populations. f4-statistics was performed with qp4pop from the ADMIXTOOLS v.6.0 package [START_REF] Patterson | Ancient admixture in human history[END_REF]. Admixture dates were estimated with MALDER v. 1.0 for all negatively significant scenarios from the f3statistics analysis [START_REF] Loh | Inferring admixture histories of human populations using linkage disequilibrium[END_REF]. All dates mentioned in the manuscript are converted into years based on 30 years per generation. TREEMIX v. 1.12 analysis [START_REF] Pickrell | Inference of population splits and mixtures from genome-wide allele frequency data[END_REF] was performed by setting blocks to 500 SNPs, accounting for linkage disequilibrium, and migration edges were added sequentially until the model explained 99% of the variance, rooting the tree with Chimpanzee (panTro5).

Scenarios of the settlement of Sahul were tested for the full dataset using qpgraph from the ADMIXTOOLS v.6.0 package [START_REF] Patterson | Ancient admixture in human history[END_REF]. One subgroup was chosen as representative population for a given regional cluster (as defined in Supplementary table 2), based on the number of individuals, the number of SNPs after masking for Asian SNPs and their location. Populations were added one-by-one, testing different scenarios with or without admixture events, following recommendations in [START_REF] Lipson | A Working Model of the Deep Relationships of Diverse Modern Human Genetic Lineages Outside of Africa[END_REF]. Models with a |Z-score|<3 were defined as fitting the data. Over one hundred different models were tested. Only models fitting the data and scenarios based on major hypotheses were represented. Results were visualized with GraphViz v. 2.46.0 (GraphViz). Population effective sizes and population divergence dates were estimated with MSMC2 v. 2.1.2 [START_REF] Schiffels | Inferring human population size and separation history from multiple genome sequences[END_REF]. Two individuals from each population were randomly chosen. Sample-specific masks were generated from the corresponding BAM files using bamCaller.py from MSMC-Tools packages [START_REF] Schiffels | MSMC and MSMC2: the multiple sequentially markovian coalescent[END_REF]. Individual masks for Asian admixture were generated from previously mentioned analyses with PCAdmix v.1.0 [START_REF] Brisbin | Sabatier user on 16 August 2021 along each chromosome in individuals with admixed ancestry from two or more populations[END_REF]. For all analyzed individual, the mask of Asian genetic tracks ranged between 0% and 49% of the total genome. The individual archaic masks were generated with ChromoPainter v.2.0 (Lawson, et al. 2012). Two metapopulations were defined: one metapopulation from 100 randomly selected Sub-Saharan Africans; and one metapopulation including Neandertal, Denisovan and Chimpanzee genomes. Only stretches of more than 4 SNPs 'painted' by the Sub-Saharan metapopulations were kept. A fourth individual mask was added based on human genome mappability hs37d5 [START_REF] Schiffels | Inferring human population size and separation history from multiple genome sequences[END_REF]. MSMC2 was run with a pattern of fixed time segments (1*5+40*1+5*3+3*5) and 25 iterations. Ambiguously phased SNPs were skipped using the -s flag. We used 1.25e-8 as the mutation rate.

Supplementary Materials

Supplementary data are available at Molecular Biology and Evolution online. All authors participated in the interpretation of the results and approved the submission.

Competing interest Declaration

The authors declare no competing interests.

All raw data are available on the European Genome-Phenome data repository:

EGAS00001005393.

Asian genetic ancestry. The two first components are represented (PC1 and PC2). (C) ADMIXTURE plot for eight genetic components summarized by genetic sub-groups (based on K=8 from Supplementary Fig. 7). The grey component combines four components mainly present in African and Eurasian groups. Numbers above bar plots identify genetic sub-groups referred in Supplementary table 2. (D) TREEMIX plot for all 58 sub-groups (based on Supplementary Fig. 7). All four panels used the same color scheme to identify the 58 genetic sub-groups. [START_REF] Summerhayes | The archaeology of forest exploitation and change in the tropics during the Pleistocene: The case of Northern Sahul (Pleistocene New Guinea)[END_REF]). The dashed black lines represent ecological lines in Wallacea.

The dark grey areas represent current land area, with light grey areas representing estimated land area around 50,000 years ago based on paleogeological reconstructions [START_REF] Coller | SahulTime: rethinking archaeological representation in the digital age[END_REF]).