New Guinea highland wild dogs are the original New Guinea singing dogs.

New Guinea singing dogs (NGSD) are identifiable by their namesake vocalizations, which are unlike any other canid population. Their novel behaviors and potential singular origin during dog domestication make them an attractive, but elusive, subject for evolutionary and conservation study. Although once plentiful on the island of New Guinea (NG), they were presumed to currently exist only in captivity. This conclusion was based on the lack of sightings in the lowlands of the island and the concurrent expansion of European- and Asian-derived dogs. We have analyzed the first nuclear genomes from a canid population discovered during a recent expedition to the highlands of NG. The extreme altitude (>4,000 m) of the highland wild dogs' (HWD) observed range and confirmed vocalizations indicate their potential to be a wild NGSD population. Comparison of single-nucleotide polymorphism genotypes shows strong similarity between HWD and the homogeneous captive NGSD, with the HWD showing significantly higher genetic diversity. Admixture analyses and estimation of shared haplotypes with phylogenetically diverse populations also indicates the HWD is a novel population within the distinct evolutionary lineage of Oceanic canids. Taken together, these data indicate the HWD possesses a distinct potential to aid in the conservation of NGSD both in the wild and under human care.

The New Guinea singing dog (NGSD) is a rare canid living in the New Guinea highlands that, in the wild, is the largest land predator on the island of New Guinea. The dogs produce a characteristic harmonic vocalization (1), described as a “wolf howl with overtones of whale song” (2). For years conservation biologists have concluded that the NGSD may be extinct or nearing extinction in the wild due to loss of habitat and the encroachment of mainland breed dogs and village dogs (1, 3, 4). In this context, nuclear DNA studies of highland wild dogs (HWD) from New Guinea, so named based on their habitat combined with the initial observation of a “wild dog” on Mount Scratchley in 1897 (5), which share striking morphologic similarity to both the NGSD and dingo, are of interest. However, because of their secretive nature and propensity to live at high altitudes distant from villages, HWD are rarely observed. Indeed, prior to 2016 they were photographed only twice (1989 and 2012) (2, 6). Although there was a report of a population of HWD in the HeLa Province of Papua New Guinea (PNG) in 2009, it was based on several assumptions and no animals were ever sighted (7).

The dogs of Oceania, unique populations found in Australia, New Zealand, and the islands nations of Melanesia, Micronesia, and Polynesia, originated from East Asian dog populations, with archaeological evidence supporting their arrival at least 3.5 kya (8). However, the dispersal timing of NGSD to New Guinea remains uncertain due to lack of archaeological evidence found on the island. NGSD were first described following collection of a specimen at an altitude of about 2,100 m in Central Province, PNG, in 1897 (9, 10). Originally classified as a distinct species, Canis hallstromi, their taxonomy remains controversial in part due to the availability of only captive specimens for genetic analysis and debate regarding their origin (3, 10, 11). Though genetically similar to the dingo, the NGSD represents a distinct population, as evidenced by both morphology and behavior (12). Resulting from a very small founder population, no more than 200 to 300 captive NGSD remain alive today, largely bred for conservation purposes. Therefore, the population of free-roaming HWD may not only represent a significant evolutionary unit important for conservation and management but possibly an important link to understanding dog domestication (3).

The need to study this presumably archaic lineage of dog is threefold. First, it is critical to determine if the captive NGSD are from the same population as the HWD or, alternatively, if the NGSD is truly extinct in the wild. Second, the primary aim of the NGSD Conservation Society (https://ngsdconservation.org/home/) is to maintain as much genetic variation within the population as possible. This is extremely challenging given the limited gene pool of the captive NGSD which descend from only eight partially related founders (3, 13). If the existing HWD does indeed represent the predecessor population of the famed NGSD, it is imperative that it be protected. Breeding programs should be established to infuse what is expected to be comparatively robust genetics of the free-breeding HWD into the genetically compromised captive NGSD population, for which inbreeding coefficients have been estimated at >0.50 (14). Finally, truly wild populations of dogs are rare, including the dingo and potentially the HWD, highlighting the need for rapid and well-organized conservation plans.

In 2016, an expedition led by the New Guinea Highland Wild Dog Foundation (NGHWDF) in collaboration with the University of Papua reported the existence of 15 HWD on the western side of the island of New Guinea near the open-cut Grasberg Mine (15). Photographs and fecal samples were collected but were insufficient for nuclear genome analysis. A subsequent 2018 field study led to the collection of blood samples from three putative HWD in their natural environment, as well as demographic, morphologic, and behavioral data. We have utilized those samples to produce a detailed analysis of the HWD nuclear genome, enabling us to determine the relationship between the captive NGSD and modern HWD and to answer the question of whether or not the NGSD is extinct in the wild.

Results

This research study was initiated as an affiliation and collaboration between the NGHWDF and the University of Cenderawasih, in Papua Province, Indonesia. Samples were collected during a 2018 expedition to the base of Puncak Jaya, within the Tembagapura district in the Mimika Regency of Papua, Indonesia. Wild dogs were conditioned to approach and enter cage traps over a period of 2 wk, during which time 18 dogs were observed, 10 of which were new to the study and 8 of which were observed in 2016 (15). Two captured dogs were fitted with Global Positioning System collars and released. All dogs fit the general description of an NGSD (Fig. 1) and the gross body measurements were within or very close to the expected range based on the small number of phenotyped, captive NGSD described previously () (10). All three samples were used for the genomic studies described below.

Fig. 1.

Oceanic dog populations compared to breed dogs and wild canids. (A) Images from left to right: HWD seen during the 2018 expedition, NGSD, dingo. (B) Neighbor-joining dendrogram of dogs from 161 breeds places the HWD (green) within the clade of East Asian (gold) and Arctic (yellow) breed dogs on a monophyletic clade with the other Oceanic dog populations (NGSD in red, dingoes in blue on three branches indicated with an asterisk). Branches with 100% bootstrap values are marked with a black dot. (C) Box plots indicating the distribution of total haplotype sharing between all pairs of dogs from different populations. Each graph represents one population, indicated above the graph, sharing haplotypes with all others. Each box represents a breed. Haplotype sharing with individuals from the same breed is not included. The three Oceanic breeds are between the two gray lines in order from top to bottom dingo, NGSD, and HWD. The rest of the populations are grouped and colored as follows: black, nondog canids; yellow, Asian/Arctic origin, orange, Nordic origin; brown, Mediterranean origin; pink, Western European origin; gray, mixed or unclassified. The full list of breeds in order can be found in .

Combining new and publicly available data, we created a dataset of 1,346 dogs from 161 breeds as well as 9 nondog canids, 16 captive NGSD, and 25 wild dingoes (14, 16, 17). Phylogenetic inference of this genotypic dataset places HWD as a monophyletic clade with 100% bootstrap support (Fig. 1). This clade is positioned adjacent to the captive NGSD and dingoes on a single branch comprising all of the dog populations we sampled that originate in Oceania. The branch is within the Asian dog clade, which also includes purebred dogs of East Asian and Arctic ancestry. Ninety-six percent of the branches (25/26) in this clade were positioned with 100% confidence. Dingoes representing multiple wild populations in Australia are not monophyletic.

To determine the genetic composition of HWD compared to all available dogs we assessed possible hybridization using identical-by-descent (IBD) haplotype sharing. Haplotypes were phased across the genome and those inferred to be IBD with a logarithm of the odds (LOD) score >3.0 were summed for every pair of dogs. The distribution of those sums was plotted for each breed dog or wild canid-to-query pair where the query is either HWD, NGSD, or DING (Fig. 1). HWDs display significant haplotype sharing with captive NGSDs and dingoes, but not with any modern breed dogs. However, the background level of sharing with non-Asian breeds is significantly higher in the HWD than in captive NGSD or dingoes (Wilcoxon P value < 2.2e-16) ().

We used the phased haplotypes to assign regions of the HWD genome to the most similar group representing an ancestral population using RFmix (Fig. 2). By assigning each allele to exactly one of these representative populations based on haplotype, we estimate that 72% of the HWD genome is most like captive NGSD or dingoes, the representative groups for Oceanic populations (Fig. 3). This is in stark contrast to the village dogs of New Guinea which share 87% of their genome with breed dogs and only 13% with the Oceanic dogs. An even weaker Oceanic signature is evident in village dogs from Vietnam (11%) and Namibia (<1%). The Australian cattle dog, a historically hybrid breed population, shows 1.5% Oceanic dog heritage using the same metrics. The HWD is predicted to share only 28% of its genome with breed dogs. The excess of European-specific haplotype signatures in all of the dogs from PNG suggests admixture as the cause. Alternatively, it is possible that any ancient variation found in the HWD could masquerade as a signature of recent admixture, as the dog populations representing the Oceanic branch have diverged in ways that would alter their allele frequencies from that of the HWD. The dingoes are as many as 1,200 generations removed from the original dogs of New Guinea, providing adequate time for drift and bottlenecks to allow for loss of diversity. In addition, the captive NGSD encapsulate a severely limited amount of the genomic variation of the original wild population. Thus, the allelic imbalances may, in part, be lost variation in the modern representatives of the ancestral line.

Fig. 2.

Assignment of HWD chromosomes to their ancestral geographic origins. Each HWD was phased and each marker was assigned an origin from the modern populations provided based on haplotype sharing with individuals from that population group: red, Oceanic; yellow, Asian; light blue, Western European; dark blue, Mediterranean; pink, Nordic; black, wolf. The graph displays the number of alleles among the six available that were assigned to each representative ancestral population at each site and is ordered by representative ancestral population, not by individual HWD. Only four chromosomes are represented on the X, as two of the HWD are male.

Fig. 3.

Assessment of admixture in HWD. (A) Bar plot of admixture in HWD compared to village dogs and a putatively admixed breed. Admixture is based on the assignment of SNPs to putative ancestral chromosomes by haplotype sharing with modern representatives and averaged across all individuals of the population: 3 HWD and 10 each of village dogs from New Guinea (VNG), Vietnam (VV), and Namibia (VN) and Australian Cattle dogs (ACD). Red, Oceanic; yellow, Asian; light blue, European; dark blue, Mediterranean; pink, Nordic; black, wolf. (B) PCA of 251 dogs and nine wild canids shows division between Oceanic and other dogs on PC1 and wild canids and domestic dogs on PC2. Geographic origin of the dogs is indicated by color. (C) PCA analysis of only Oceanic dogs reveals separation of pure Oceanic and hybrid village dogs along PC1 and separation of New Guinea dogs and Australian dingoes along PC2. (D) Maximum likelihood trees of geographic-based breed populations with migration. Migration is indicated by the arrows shaded from yellow to red according to their weight. Golden jackal was used to root the tree. (E) Maximum likelihood tree with New Guinea village dogs added and migration allowed.

Observed heterozygosity in the HWD is not significantly different from that of breed dogs from any region but is significantly higher than captive NGSD (pwilcox = 0.0021), or a mixed population of dingoes (pwilcox = 0.0006) (). It is also significantly lower than that of village dogs collected in New Guinea (pwilcox = 0.0023). Because the entire captive population of NGSDs is descended from only eight individuals, low genetic variation is expected. However, the dingoes derive from a varied, free-breeding population and would be expected to show much more heterozygosity. The low levels detected are likely an indication of ascertainment bias in the Illumina CanineHD single-nucleotide polymorphism (SNP) set, as the markers were originally chosen because they are polymorphic in European breed dogs and were not widely tested outside of that population. Similarly, the small number of wolves from different parts of the world analyzed with these SNPs show reduced heterozygosity (P = 9e10−6 compared to village dogs) and increased inbreeding over what is predicted in a free-breeding population but is not significantly different from that found in the HWD ().

To perform principal components analyses (PCA) we trimmed the dataset by removing SNPs in high linkage disequilibrium with other SNPs. To help balance the comparisons, we reduced the breed dogs to four groups of 46 to 48 dogs each, representing major geographic regions of breed development; East Asia and the Arctic (ASIA), northern Europe (NORD), western and central Europe (EURO), and the Mediterranean region (MEDI), including parts of southern Europe, northern Africa, and central Asia. Breeds displaying recent admixture across geographic regions were excluded (16). PCA shows distinct clustering of the Oceanic dog groups apart from the breed dogs and wild canines (Fig. 3) along the first PC, which accounts for 32% of the total variation in the dataset (pT-W = 2.03e-297). The village dogs from New Guinea cluster closer to the breed dogs than the pure Oceanic populations. Analysis of the Oceanic group alone reveals distinct, nonoverlapping groupings of captive NGSD, dingoes, HWD, and village dogs (Fig. 3). The first PC separates the village dogs from the NGSD, HWD, and dingoes (77% of variation, pT-W = 7.70e-15). The second PC separates the NGSD from the dingoes (23% of variation, pT-W = 5.83e-15), with the HWD trending toward the NGSD.

Using allele frequencies from the independent genotypes, we built maximum likelihood trees of the Oceanic populations and geographic breed groups to assess the likelihood of admixture using the program Treemix (18). This method predicted the migration of an ancestor of all European breed dogs into the HWD, as well as compensatory migration from the New Guinea dogs to breed dogs (Fig. 3). These same migration events are predicted if New Guinea village dogs are added to the tree, as well as a third migration of New Guinea naturalized dogs to the village dogs (Fig. 3).

We tested these predictions further by calculating f3 which, when negative, provides evidence of admixture from two source populations into a third. We observe such evidence in HWD when the source populations include captive NGSD and any of the four breed groups (Table 1).

Table 1.

Detection of possible admixture between Oceanic populations and breed dogs in the HWD using the F3 statistic

Source 1: Breed	Source 2: Oceanic	Target	Z-score
ASIA	NGSD	HWD	−4.708
EURO	NGSD	HWD	−15.289
MEDI	NGSD	HWD	−13.639
NORD	NGSD	HWD	−13.871
ASIA	Dingo	HWD	n.s.
EURO	Dingo	HWD	n.s.
MEDI	Dingo	HWD	n.s.
NORD	Dingo	HWD	n.s.
All	All	Dingo	n.s.
All	All	NGSD	n.s.

“All” indicates all possible combinations of populations that fit the column descriptions. n.s = not significant (Z > 0).

To further confirm the source and degree of possible admixture into the HWD we calculated D statistics and F4 ratios from four-branch tree arrangements including the HWD using either golden jackal or wolf as the outgroup. We find significant evidence for admixture (Z > 3.0) between all breed groups and the HWD (). F4 ratios of tree arrangements swapping any of the breed groups for the HWDs predicted admixture ranging from 15 to 50%; however, the z-scores were insignificant for all combinations. We next used qpWave to examine the number of likely populations contributing to the HWD and determined that a two-source model best fits the data provided. That number did not change, regardless of which breed groups or Oceanic populations were included in the analysis. We ran qpAdm and obtained a prediction of 29.6% breed dog contribution to the HWD. This is very similar to the 28.8% migration from the breed clade predicted by Treemix and the 28% nonoceanic contribution predicted from haplotype sharing.

Discussion

The 2018 expedition to Indonesian New Guinea was motivated by analysis of specimens collected in 2016 (15). Specifically, analysis of mitochondrial DNA isolated from 9 of 24 presumed HWD scat samples revealed the A29 haplotype, rather than the A79 haplotype. The A79 haplotype is unique to NGSD, while the A29 haplotype is observed in dingoes, NGSD, some Asian and Arctic breed dogs (19, 20), and village dogs. Thus, the results are not specific to the NGSD and can be generalized to dogs from East Asia and the Oceanic region (13, 21–23).

Three contemporary HWDs provided a look into the nuclear genome of this elusive, potential originator population that shares striking morphological similarity to both the NGSD and dingo. Distance measures place the HWD split basal to the dingo and NGSD division in a monophyletic clade with 100% bootstrap support. In addition, maximum likelihood trees, allowing for a small amount of migration, place the HWD, NGSD, and dingo on a branch excluding all breed dogs and diverging prior to the formation of breed dogs from four diverse geographic origins. These results suggest a distinct evolutionary lineage for the Oceanic populations. This is confirmed by the extensive haplotype sharing within dogs originating in Oceania to the exclusion of the remainder of canids. This shared genomic ancestry, along with estimates of genetic diversity, highlight the lack of an ancestral contribution from any lineage but Oceania, and specifically HWD, to the current captive NGSD. The lack of diversity renders it difficult to estimate the specific genetic origin of NGSD; however, the extant population of HWD is most similar of all populations thus far studied.

Using F3 statistics, admixture is detected in the HWD only when NGSD are part of the analysis. This suggests that the captive NGSD are derived solely from the HWD and not from another Oceanic population, that the HWD are not a new introduction of dingo to the islands, and that the presumptive admixture may represent early variation missing from the captive NGSD. The pattern of allele sharing is not found when HWD are analyzed with dingoes and breed dogs, suggesting that the derived alleles shared with breed dogs are similarly present in the dingo, and are missing only in the very restricted NGSD genome, likely as a result of reduced heterogeneity rather than recent migration.

As our knowledge of the genetic history of the dingo advances (14), it becomes increasingly necessary to understand the underlying contribution of phenotypically similar and closely related canids to the Oceanic dogs of the Canidae. Phylogenetic and haplotype analyses provide compelling evidence that the HWD studied here represent extant animals arising from the same historic stock as the captive NGSD. Analysis of admixture suggests that there is evidence for minimal introgression in the HWD, likely from a population ancestral to modern breeds. The fact that the New Guinea population of HWD possesses much higher genetic diversity than their largely inbred, captive counterparts argues that conservation efforts to revitalize the NGSD would benefit from the utilization of HWD in combination with the captive NGSD. Indeed, the inclusion of even a small number of individuals from a natural population can have a significant effect on the genetic diversity of a captive breeding population (24).

The size and extent of the HWD population, however, remains to be documented. Analyses in this work encompassed only three individuals. Sampling of dogs from more remote regions along the central range such as Puncak Mandala, Puncak Trikora, and equally remote areas on the PNG side of the Island of New Guinea would provide insight as to both the numbers of HWD in their natural habitat and their genetic diversity. Conservation efforts will benefit most from inclusion of the greatest number of specimens that best represent the original dogs, with the least amount of influence from outside sources, making it imperative that these studies be continued. In addition, studies of small pockets of wild living dogs will enable the search for the genetic variation that is required for free living, but irrelevant in a captive environment, and vice versa. There is urgency in this exercise, as each new generation of captive NGSD increases the risk for fixation of deleterious alleles, thus compromising subsequent breeding programs (25). This study provides insight into the genetic makeup of an isolated canid population and thus identifies a readily available genetic reservoir for the NGSD, previously thought to be extinct in the wild.

Materials and Methods

Sample Ascertainment.

This research study was sponsored by PT Freeport Indonesia and resulted from an established collaboration between the NGHWDF, under the sponsorship of PT Freeport Indonesia, and the University of Cenderawasih, in Papua Province, Indonesia. The Foundation has participated in two research trips, the first in 2016 (15) and the second in 2018, both to the base of Puncak Jaya, which is within the Tembagapura district in the Mimika Regency of Papua, Indonesia. Samples utilized for this study were collected in 2018 from locations at the periphery of the Grasberg Mine using a novel capture approach for HWD. Wild dogs were conditioned to approach and enter cage traps over a period of 2 wk following enticement with crushed coyote gland lure and skunk essence. Traps were baited inside and out of cages. Auditory North American coyote calls were played to attract wild dogs and elicit return howls. Trail cameras reviewed the progress. After 2 wk the traps were set and the capture campaign was initiated. Two wild dogs were captured, immobilized, and underwent extensive examinations including measurements () and collection of biologic samples. The dataset included DNA and information from those dogs as well as a third recently deceased dog that appeared to be the victim of a car collision.

DNA Preparation.

Three HWD samples were provided for processing. Two blood draws from live dogs were stored in acid-citrate-dextrose anticoagulant. DNA was isolated using described methods (26). The third sample, skin from the deceased dog, was stored in RNAlater (Qiagen) for transport. DNA was isolated using the DNeasy Blood and Tissue micro kit (Qiagen) and standard protocols following homogenization. All protocols completed at NIH and Texas A&M University followed Institutional Animal Care and Use–approved protocols.

SNP Genotypes and Phylogeny.

DNA samples from the HWD as well as 10 captive NGSD were genotyped on the Illumina canine HD SNP chip (Illumina) using standard protocols. NGSD blood samples were collected from sanctuaries and private owners who had rescued them from zoos and exotic animal breeders in seven US states and Canada. They were chosen to include as much diversity from the captive population as possible. Genotype calls were made using Genome Studio v2.0.4 with genotyping module v2.0.4 (Illumina). Data from captive NGSD and HWD (Gene Expression Omnibus [GEO] accession no. GSE143824) (27) were analyzed with 1,346 domestic dogs representing 161 breeds, 7 wolves, and 2 golden jackals genotyped at 150,112 SNPs spanning all autosomes and the X chromosome (GEO accession nos. GSE90441, GSE83160, GSE70454, and GSE96736) (16, 28–31). In addition, publicly available SNP chip data from 24 dingoes (Dryad submission DOI:10.5061/dryad.sq8d0) (14, 32) was incorporated into the analysis, as was data from an additional five NGSD and 60 village dogs collected in New Guinea, Vietnam, and Namibia (Dryad submission DOI:10.5061/dryad.v9t5h) (17, 33).

Linkage disequilibrium pruning was carried out with an r2 threshold > 0.5, as assessed by PLINK v1.9 using the –indep-pairwise option (34, 35). This reduced the dataset to 107,125 variants for subsequent principal component and admixture calculations.

To create a balanced dataset for PC and admixture analyses, purebred dogs were grouped according to ancestral geographic region of origin based on published phylogenetic results (16). Four groups were identified: Nordic, 46 dogs from breeds originating in Scandinavian countries; Mediterranean, 47 dogs from breeds originating in southern Europe, central Asia, and northern Africa; European, 46 dogs from breeds originating in Western European countries and the United Kingdom,; and Asian, 48 from breeds originating in East Asia and Arctic regions. Breeds that showed recent haplotype sharing across geographic groups were excluded. Each group was then trimmed to simulate the size of the smallest group, retaining as many different breeds as possible while randomly choosing representatives of the breeds. These were combined with the Oceanic populations, 25 dingoes, three HWD, and 16 NGSD, for the final balanced dataset.

Phylogenetic Analysis.

A pairwise identity-by-state distance matrix was computed using PLINK v1.9 and the –distance 1-ibs command (34, 35). Bootstrapped distance matrices were created by randomly resampling 150,112 markers with replacement 100 times and input into PHYLIP using neighbor and consense to construct neighbor-joining dendrograms (36). Dendrograms were visualized using FigTree v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/).

PCA.

PCs were calculated for the balanced dataset along with the wolves and jackals and village dogs from New Guinea. A second PCA was performed including only dogs from the Oceanic region: HWD, captive NGSD, Dingo, and New Guinea village dogs. The eigensoft package was used to calculate PCs with smartpca (37, 38). Twstats was used to determine the significant components using the Tracy–Widom distribution. Proportion of variance was calculated from only the significant PCs (PT-W < 0.05).

Haplotype Sharing.

Data were phased and haplotypes determined using the program Beagle v4.1 (39) with sliding windows of 1,000 SNPs and a 50-SNP overlap. Haplotypes that are identical-by-state and that are predicted to be IBD, based on excess length and sharing (40–42) with a LOD score greater than 3.0, were calculated between every pair of dogs using the ibd option of Beagle (40). The inferred IBD haplotypes for each pair were summed across the genome to assess the recent breed ancestry of the HWD, captive NGSD, and dingoes.

Chromosomal regions were mapped to the most likely geographic ancestor population by RFmix (43). Phased genotypes from the balanced dataset were used after collapsing haplotypes that were separated by only one SNP and reassigning phase to the original Beagle output as described in Browning et al. (44). Total estimated ancestry was calculated by averaging the ancestral assignments of all 150,000 markers.

Admixture and Migration.

The balanced dataset described above was used to assess the presence and extent of admixture in the HWD using Treemix (18) and AdmixTools (45). Wild canids were included in the dataset with the golden jackal serving as the outgroup for the maximum likelihood trees. Allele frequencies for each population were calculated in PLINK. Treemix was run with the bootstrap option and k = 1,000, allowing 0 to 10 migration events and repeated five times at each value. The optimal number of migrations was reached when the average variance came within one SD of 99.8% as calculated in the program OptM (https://rdrr.io/cran/OptM/) and by the Evanno method (46), the latter of which is calculated from the change in likelihood (∆m) with each additional migration event (). This analysis was repeated including village dogs collected in New Guinea (17). Maximum likelihood trees with predicted migration events and residuals were viewed using plotting_functions.R, provided with the Treemix package.

F3, D statistics, and F4 ratios were calculated from three and four branch tree arrangements to evaluate the ancestral components of HWD and geographical groups of breed dogs. F3 was calculated using threepop from Treemix (18). D-statistics and F4 ratios were calculated using Admixtools with golden jackal as the outgroup. Because all of the breed groups can be assumed to have some shared ancestry with one another, wolf was included as the unrelated branch in the F4 ratio calculations (45).

Also using the R package admixr (47) to run AdmixTools we calculated the most likely number of source populations contributing to the HWD and the extent of estimated hybridization. These were calculated using the qpWave and qpAdm commands. Because no one breed group was determined to be the source of admixture, all breeds were combined to obtain the percent hybridization. Therefore, when running qpAdm the target was set as HWD, the sources were Oceanic and breed dogs, and the outgroups were the wolves and jackals.