Has the introduction of two subspecies generated dispersal barriers among invasive possums in New Zealand?

The introduction of species into new environments provides the opportunity for the evolution of new forms through admixture and novel selection pressures. The common brushtail possum, Trichosurus vulpecula vulpecula from the Australian mainland and T.v.fuliginosus from Tasmania, were introduced multiple times to New Zealand from Australia to become one of New Zealand’s most significant pests. Although derived from two subspecies, possums in New Zealand are generally considered to be a single entity. In a previous analysis, we showed that possums in the Hawkes Bay region of New Zealand appeared to consist of at least two overlapping populations. Here, we extend that analysis using a genotype-by-sequencing approach to examine the origins and population structure of those possums and compare their genetic diversity to animals sampled from Australia. We identify two populations of each subspecies in Hawkes Bay and provide clear evidence of a contact zone between them in which a hybrid form is evident. Our analysis of private alleles shows higher rates of dispersal into the contact zone than away from it, suggesting that the contact zone functions as a sink (and hence as a barrier) between the two subspecies. Given the widespread and overlapping distribution of the two subspecies across both large islands in New Zealand, it is possible that many such contact zones exist. These results suggest an opportunity for a more targeted approach to controlling this pest by recognising sub-specific differences and identifying the contact zones that may form between them.


Introduction
Introductions of species into novel environments can provide the opportunity for the evolution of new forms as they naturalise (Sz} ucs et al. 2017). In some cases, alien species have been introduced to the same location on multiple occasions from different sources, creating genetic combinations of individuals that have not previously existed (Kolbe et al. 2004). Such admixture can enhance fitness through hybrid vigour or increased genetic diversity (Barker et al. 2019;Dlugosch, Parker 2008;Facon et al. 2008;Hirsch et al. 2017;Rius, Darling 2014) or reduce fitness through outbreeding depression (Rhymer, Simberloff 1996). It is unclear how frequently genetic factors are intrinsic to the success or otherwise of invasive species (Barker et al. 2019;Rius, Darling 2014) although propagule pressure (something that will affect the levels of genetic diversity within an invader) is one of the few acknowledged general characteristics of successful invaders (Bock et al. 2015). Identifying the role that such factors can play in successful establishment, may offer novel perspectives or pathways for the identification of high-risk invaders or provide opportunities for enhanced control of invasive pests (Sz} ucs et al. 2017).
The introduction of common brushtail possums (Trichosurus vulpecula; hereafter called possums) to New Zealand from Australia is one well documented instance where different forms of the same species have been brought together into a new environment. Introduced to establish a fur trade (Montague, 2000), possums were shipped to New Zealand from Tasmania and several other undocumented locations on mainland Australia, most likely in New South Wales and Victoria, between 1858 and the 1920s. All introductions ceased after that time but the transfer of possums within New Zealand continued into the 1940s. Possums are endemic to Australia with five currently recognised subspecies distributed across the continent and on several offshore islands (How, Hillcox 2000;Kerle et al. 1991). New Zealand has no naturally occurring marsupials and possums are now one of the most serious threats to biodiversity and agriculture in New Zealand damaging vegetation (Cowan 1991;Nugent et al. 2001); preying on birds and their eggs (Brown et al. 1993); and spreading bovine tuberculosis (Innes, Barker 1999). New Zealand management agencies have conducted extensive possum control for decades (Donnell 1995;Morgan 1990;Nugent, Cousins 2014;Nugent, Morriss 2013) and now aim, with the support of conservation groups, to eradicate brushtail possums as part of the ''Predator-free New Zealand'' campaign (Russell et al. 2015).
It is well known that the different colour forms seen in New Zealand possums are representative to some extent of their different ancestry (Kean 1971;Pracy 1974;Triggs, Green 1989), but it is generally assumed (although seldom stated) that New Zealand possums breed indiscriminately with respect to those origins (Montague 2000). However, a recent microsatellite DNA analysis of possums in Hawkes Bay, on the North Island of New Zealand, identified at least two genetically distinct groups of possums, separated by a zone of putatively introgressed animals (Sarre et al. 2014). That work implied that there were at least two breeding forms in the region, probably representing the Tasmanian subspecies and the south-eastern mainland subspecies (T.v.fuliginosus and T.v.vulpecula respectively), rather than one continuous population (Sarre et al. 2014). That proposition is supported by historical records which show that there were three recorded introductions of possums into the area ( Fig. 1): possums taken directly from Tasmania and released in the Lake Waikaremoana area in 1898; ''black furred'' possums (also presumed to be Tasmanian in origin but translocated from elsewhere in New Zealand) released on the Mahia Peninsula a short time after the first introduction; and possums from an unknown location on mainland Australia released in the Mohaka river area in the early 1900s (Pracy 1974;Sarre et al. 2014).
Resolving complex introduction histories and patterns of admixture in invasive species can be difficult, but access to population level genomic data can reveal how species integrate into new environments and whether diverse origins have enhanced or inhibited their invasion success. Here, we use genotype-bysequencing to examine the origins and contemporary population structure of possums in Hawkes Bay, and compare their genetic diversity to animals sampled from across their native range in Australia. We test the proposition that the sub-populations, as previously defined by the distribution of microsatellite DNA variation (Sarre et al. 2014), are representative of the ancestral mainland and ancestral Tasmanian origins, and that the Tasmanian and mainland forms have established a hybrid swarm in the Hawkes Bay region that affects the direction and rate of gene flow. We also use asymmetries in the number of private alleles between pairs of sub-populations to determine the direction of that geneflow. A symmetrical distribution of private alleles among sub-populations indicates symmetry in the rates of gene flow between subpopulations and subspecies and would argue against the formation of genetic barriers between the two subspecies. We demonstrate that one sub-population is genetically identifiable as the Tasmanian subspecies T.v. fuliginosus, that possums from that sub-population are distinct from animals that are most likely mainland in origin, and that they have maintained that distinctiveness even after more than 80 years of contact. We also show that the sub-population believed to originate from mainland Australia is indeed largely T.v.vulpecula in origin, although we could not identify the mainland region from which those possums were obtained. Our data support the proposition that brushtail possums in Hawkes Bay do not behave as a single panmictic population, but rather as two subspecies with a narrow zone of contact. Gene flow into and out of that contact zone is asymmetrical and indicative of a genetic barrier between the two subspecies. The formation of similar hybrid barriers elsewhere in New Zealand could have important implications for the control of this invasive species.

Sample selection
In order to determine the Australian origins of the invasive possums at our study site in Hawkes Bay, New Zealand, we obtained samples from a broad range of locations in eastern mainland Australia (n = 19) and Tasmania (n = 11). Specifically, we sampled tissue (mainly ear tissue) from museum specimens and road killed animals, from six states and territories (Northern Territory, Australian Capital Territory, Queensland, New South Wales, Victoria and Tasmania) and Flinders Island (located in the Bass Strait between Tasmania and mainland Australia; Fig. 2). These samples encompass the four commonly recognised subspecies of Trichosurus vulpecula in eastern and northern Australia (Kerle et al. 1991): T.v vulpecula (south-eastern Australia), T.v. arnhemensis (Northern Territory), T.v. johnstonii (Queensland) and T.v.fuliginosus (Tasmania and Flinders Island) and represent the subspecies most likely to include the source populations for New Zealand possums (Fig. 2). Possums from the Hawkes Bay region of New Zealand were sampled from 24 sites (n = 253) and were selected as a subset of the samples previously genotyped using microsatellite DNA (Sarre et al. 2014).
Specifically, we applied New Hybrid analysis (Anderson, Thompson 2002) to the microsatellite data previously reported by Sarre et al. (2014) to estimate the percentage of ancestry for each individual possum and indicate whether an individual is an F1 or F2 cross between subspecies. We then selected individual possums that represented all possible combinations of F1, F2 and 'pure' ancestry, as well as an even number of males and females, and used the NewHybrid analysis to sort all 1,605 individuals into three genetic groups: (1) largely Tasmanian in origin, (2) largely Australian mainland in origin, and (3) of mixed ancestry. We then stratified for collection site within those three groups, focussing most on samples from the sites closest to the three known sites of introduction and the previously identified contact zone (Fig. 2), and then randomly chose individuals across those strata for inclusion in the present study using genotyping-by-sequencing.
Genotype-by-sequencing DNA was extracted from each 2mm 2 piece of ear tissue sample by Diversity Arrays Technology Pty Ltd (DArT) using a salting out method as reported by (Couch et al. 2016). We used the DArTseq TM genotype-by-sequencing approach, which uses a combination of restriction enzyme complexity reduction and high throughput sequencing (Courtois et al. 2013;Cruz et al. 2013;Kilian et al. 2012) to simultaneously identify and genotype SNP markers in the absence of a reference genome. A double digest approach with the restriction enzymes PstI and SphI was used for complexity reduction before the addition of custom adapters (Kilian et al. 2012) to restriction site overhangs. Fragments that included both a PstI and an SphI adapter were amplified using primers complementary to the adapters. These also incorporated molecular identifier barcode tags, to allow multiplexing of up to 96 samples per sequencing run. PCR conditions consisted of: denaturation at 94 0 C for 1 min; 30 cycles of 94 0 C for 20 s, 58 0 C for 30 s and 72 0 C for 45 s; and a final extension period of 72 0 C for 7 min. PCR products were pooled for sequencing on an Illumina HiSeq 2000 platform using 77 cycles of single end sequencing. Raw sequence reads were filtered and processed using a proprietary DArT analytical pipeline. Poor quality sequences were removed, with more stringent criteria being placed upon the barcode region than the rest of the sequence. Approximately 2,000,000 sequences per barcode were identified and used in marker calling, during which identical sequences were collapsed and filtered before screening to identify variable markers using DArT proprietary SNP and SilicoDArT algorithms (DArTsoft14).

Analysis of genotype-by-sequencing data
The R package dartR (Gruber et al. 2018) was used for population genomic analyses of the DArTseq SNP data, including F-statistics, principal coordinates analysis (PCoA) and fixed difference analysis. SNPs were filtered using the RepAvg function (Reproducibility C 0.95), CallRate (0.95), with lower and upper read depth thresholds of 5 and 90 respectively, and to exclude monomorphic loci and secondary SNPs.
We estimated ancestry proportions using SNP data and geographic location in tess3r (Caye et al. 2016), for K values of between two and eight with the default settings (code supplied on request). Fixed differences were identified between the New Zealand and Australian (mainland and Tasmanian combined) possum populations and through a second analysis in which we considered only the sub-populations at the Hawkes Bay study site identified using PCoA and tess3r (Caye et al. 2016).
To further define and visualise the population genetics of animals in Hawkes Bay, we ran a STRUCTURE analysis (Pritchard et al. 2000) using all loci with identical filters (n = 502) for each value of K (assumed number of ancestral populations) between 2 and 5 with 10 repeats each of 100,000 iterations. The first 50,000 iterations were disregarded as burnin and then each of the 10 repeats were summarised using the programs Evanno (Evanno et al. 2005) and Clumpp (Jakobsson, Rosenberg 2007).
We ran STRUCTURE using both ancestral models (No Admixture and Admixture) because they differ in their underlying assumptions about the origins of the population under study (Porras-Hurtado et al. 2013). The No Admixture scenario assumes no prior knowledge about the origin of populations under study. In our case, we know something about the origins of the populations drawn from the introduction history and we have an hypothesis arising from a previous microsatellite DNA study (mainland, Tasmanian and Hybrid) so for this analysis we ignored that knowledge and checked to see if the three groups still emerged. We then ran the Admixture model using our prior knowledge of the origin of the populations to provide insight into the proportion of individuals that can be attributed to each of the presumed ancestral populations. In order to do this, we combined individuals from those sites closest to the mainland introduction sites, those from the previously identified hybrid zone (Sarre et al. 2014) and those closest to the sites populated by possums from Tasmania into subpopulations A, B and C respectively. We expected that A and C (Mainland and Tasmanian origin) would each have a unique signature, while individuals from population B should show a mixture of both A and C.
We also used a NewHybrids analysis to examine population structure among Hawkes Bay possums. While broadly similar to that provided by STRUC-TURE, NewHybrids analysis differs in that it uses a detailed inheritance model to check to see if individuals are indeed hybrids when two potential source populations mainly comprising of pure individuals can be identified (Anderson, Thompson 2002). As NewHybrids can only be applied to a maximum of 200 loci, we reduced the available loci to 15,523 by filtering on call rate (C 0.95), and then selected the most informative (based on their polymorphic information content as implemented in dartR) 200 SNPs for NewHybrids analysis, using 30,000 sweeps with a burn in of approximately 100 sweeps, and using the Jefferies prior for both Pi and Theta. To further validate the result and its independence from the 200 loci identified as most informative, we re-ran the analysis with 200 randomly selected loci 30 times.

Directionality of gene flow through asymmetry in private alleles
To further determine the direction of interactions within the Hawkes Bay study site, we estimated the levels of asymmetry of gene flow between subpopulations using a novel approach based on the numbers of private alleles present in each subpopulation. Specifically, we used the idea that a subpopulation that has a net loss of individuals towards another sub-population through emigration will lose private alleles over time as it shares them with the recipient sub-population. Under this process, the subpopulation that has fewer emigrants will lose private alleles (through sharing) at a slower rate that its paired sub-population while also accumulating private alleles through genetic drift. We used this logic to explore the direction of gene flow by counting the number of private alleles between sub-populations across the study area while accounting for differences in sample size. This approach differs from those developed for microsatellite DNA markers (Kennington et al. 2003;Sundqvist et al. 2016) by using the raw counts of private alleles rather than comparisons between observed and simulated asymmetries. Under our approach, we consider that a detectable asymmetry in the number of private alleles between neighbouring sub-populations is indicative of asymmetric geneflow between those two sub-populations. Conversely, symmetry in the numbers of private alleles is indicative of equal rates of gene flow in both directions between two sub-populations.
To study the gene flow into and out of the contact zone using private alleles we clustered adjacent sample sites with the aim of obtaining an approximately even number of individuals per cluster. We used the proximity of individuals as the main criterion for clustering, while also considering the genetic similarity of individuals as indicated by the NewHybrid approach (see above) using microsatellite DNA data gathered previously by Sarre et al. (2014). Using this approach, we generated six clusters (A1, B1, B2, B3, C1, C2) within the three sub-populations and confirmed their validity by running a Discriminant Analysis of Principal Component (Jombart et al. 2010). We then counted the number of private alleles for each pairwise comparison among the 6 subpopulations. Sample sizes of individuals in each subpopulation varied between 20 and 36. To account for the uneven number of individuals, we sampled all subpopulations to 20 individuals (with replacement) and counted the number of private alleles in each pair of sub-samples. To reduce noise in the data set, we filtered the identified markers by minor allele frequency with a threshold of [ 0.1 (based on 40*0.1 [ 4 alleles) so that a cluster contained at least 5 alleles before an allele could be called private. This resulted in a data set of 506 SNPs, with a minimum minor allele frequency of 0.1 and mean minor allele frequency of 0.23. The mean number of private alleles between each resampled pair of sub-populations was then calculated from a thousand repetitions of the procedure. Finally, we calculated the mean and 95% confidence intervals of the number of private alleles and the pairwise ratio of private alleles. Sub-population pairs that had a private allele ratio larger than one, in more than 975 out of 1000 cases, were considered significantly asymmetric and hence regarded as exhibiting asymmetric gene flow.

Results
We genotyped 283 individual possums from Australia and New Zealand using the DArTseq approach, generating 27,339 SNPs, 6,596 of which were retained after filtering. Ten individuals were also removed by filtering because of high numbers of missing data leaving a total of 273 possum genotypes for analysis. We identified far fewer private alleles in the T.v.fuliginosus samples from Tasmania (mean = 25.3; n = 11) than in T.v.vulpecula sampled from mainland Australia (mean = 296.6; n = 12). Possums sampled in Tasmania also exhibited lower levels of expected (0.04) and observed (0.03) heterozygosity relative to their mainland counterparts (0.13 and 0.11 respectively). There was strong support (17.3% weight for PCoA Axis 1) for differentiation between Tasmanian individuals of T.v.fuliginosus (including those from Flinders Island) and all other individuals sampled from mainland Australia (Fig. 3). These findings are indicative of the Tasmanian history of geographical isolation, following separation of the island from mainland Australia by rising sea levels more than 9,000 years ago (Veevers 1991). There was also good support for differentiation between T. v. vulpecula and two other mainland subspecies (T. v. arnhemensis and T. v. johnsoni) although there was overlap between the latter two subspecies (Fig. 3).
Within the Hawkes Bay study area, we identified the three distinct sub-populations (A, B, and C) seen previously using microsatellite DNA analysis (Sarre et al. 2014). Sub-population A, genetically the closest to the introduction site of animals from mainland Australia, is discrete from sub-populations B and C, but does not cluster with any mainland Australian animals sampled here (Fig. 4). Sub-populations B and C cluster most closely together and have a lower pairwise F ST (0.03; Table 1) Table 1). No fixed allelic differences were observed between any pairwise combination of the Hawkes Bay sub-  (Table 2). Private alleles were observed in all pairwise combinations of groups. Overall, these data suggest that that both B and C are composed largely of descendants of Tasmanian animals. All three Hawkes Bay genetic sub-populations exhibit private alleles when compared to each other (Table 2) suggesting that gene flow is restricted between them. The smallest number of private alleles were seen between sub-populations B and C, suggesting that gene flow is higher between those two than between sub-populations A and B or between A and C.
The PCoA analysis including all mainland Australian, Tasmanian, and Hawkes Bay samples (Fig. 4) shows that the majority of mainland Australian possums are genetically different from both New Zealand and Tasmanian individuals. In contrast, the native Tasmanian possums cluster very closely with the Hawkes Bay possums of putative Tasmanian ancestry (sub-populations B and C), providing additional genetic evidence for Tasmanian ancestry for those Hawkes Bay animals. The possums from subpopulation A appear to be distinct from both Hawkes Bay sub-populations B and C and from the native Tasmanian and mainland Australian samples. We identified two mainland Australian individuals, one sampled from around the Sydney area of New South Wales and the other from the Australian Capital Territory, that cluster closest to the Hawkes Bay individuals of putatively mainland ancestry. We also used SNP genotype data combined with geographic information on each sample location to generate a visual representation of population genetic structure among the Hawkes Bay possums only using tess3r ( Fig. 5; Fig. S1). Here, we examined three scenarios for the ancestral makeup of the possums. (1) Two ancestral sub-populations (mainland Australia and Tasmania origins; k = 2); (2) Three ancestral subpopulations (T.v.vulpecula-mainland, T. v. fuliginosus-Tasmania, and the contact zone population; k = 3); and (3) Four ancestral sub-populations (T.v.vulpecul-mainland, two different introductions of T.v.fuliginosus-Tasmania, and the contact zone population; k = 4). Our estimates of ancestry proportions and ancestral allele frequencies show a clear disjunction between the two subspecies at their contact zone ( Fig. 5; Panel 1), with the contact zone being genetically closer to the populations descended from T.v.fuliginosus than to those descended from T.v.vulpecula. The contact zone itself seems to run to the west of the Mohaka River that divides the two release points (Fig. 5; Panel 2). We find no evidence of two distinct ancestral groups of T.v.fuliginosus arising from the dual introductions of possums from Tasmania ( Fig. 5; Panel 3) and it appears that this subspecies is relatively homogenous at its source within Tasmania (Fig. 4). These results were similar to those obtained using the STRUCTURE analysis (Fig. 6) which showed that the most appropriate number of groups to emerge from both No Admixture and Admixture models was three (Fig. S2) and, as predicted, demonstrated that the sub-population designated as ''B'' most closely resembled a combination of sub-populations A and C. Similarly, the NewHybrid plots ( Fig. S3; S4) show that sub-population B consists largely of F2 hybrid individuals along.
with small numbers of both parental genotypes while sub-populations A and C are comprised largely of parental individuals with mainland and Tasmanian origins respectively. There is little evidence of F1 individuals in any of the Hawkes Bay animals genotyped.
Finally, asymmetric geneflow was identified between several pairs of the six genetic clusters using our private allele approach (Fig. 7, Table S1). The number of private alleles was consistently lower in sub-population A1 (mainland origin individuals) than in the three sub-populations within the contact zone (B1, B2 and B3). The same was true for subpopulation C1 (Tasmanian origin individuals) relative to sub-populations within the contact zone (B1, B2 and B3). We also identified a significant asymmetry in the geneflow between sub-populations C1 and C2. There was no apparent effect of distance on the asymmetries between clusters (Fig. S6). In summary, the asymmetry in private alleles indicates that the direction of Table 2 Fixed differences, number of private alleles, and corrected number of private alleles between Hawkes Bay, New Zealand groups of possums (groups A, B and C) and Australian common brushtail possums sampled from Tasmania (T.v. fuliginosus) and the south eastern mainland (T.v. vulpecula).
The accumulation of fixed differences between populations indicates low levels of gene flow. The ''corrected'' number of private alleles has been divided by the larger of the two samples sizes to account for the difference between groups geneflow is into the contact zone, with only minor geneflow outwards.

Discussion
New Zealand possums are the subject of enormous national pest animal control effort and considerable scientific endeavour (Byrom et al. 2016;Clout, Sarre 1997;Eason et al. 2017;Forsyth et al. 2018;Goldson et al. 2015;Livingstone et al. 2015;Montague 2000;Rouco et al. 2018;Rouco et al. 2017). While the history of their introduction from Tasmania and the Australian mainland is extremely well documented, New Zealand possums are considered almost universally to be representatives of a single lineage, albeit one with colour variants. This view is consistent with the prevailing taxonomy that discriminates between them only at the subspecific level (Kerle et al. 1991).
Here, we demonstrate that this expectation of possums mixing indiscriminately with respect to subspecies, does not hold in practice in the introduced New Zealand population, at least within the Hawkes Bay study area. Specifically, possums in Hawkes Bay that are derived from two Australian subspecies are more closely resembling their two ancestral entities separated by a contact zone than one integrated population. This is despite having ample opportunity to disperse, mingle, and interbreed for over 80 years (between 40 and 80 generations). We found that new genetic forms, representing individuals of mixed ancestry (F2 ?) between the two subspecies, occur mainly in the intermediate contact zone between the sites of the original introductions. We also observe that gene flow tends to move inwards, into that contact zone more commonly than outwards and found no evidence of parental genotypes moving across the contact zone and into the opposite parental sub-population. As such, the contact zone that has formed appears to act as a sink and limits dispersal.
Our results provide clear evidence that possums in Hawkes Bay behave as three distinct groups rather than a single panmictic population-a result that reflects previous findings of genetic structure at the site (Etherington et al. 2014;Sarre et al. 2014). We show that Hawkes Bay possums with putative Tasmanian ancestry have high genetic similarity to modern Tasmanian populations (T.v.fuliginosus), with  only modest levels of introgression of alleles from T.v.vulpecula. In contrast, Hawkes Bay possums with putatively Australian mainland origin are genetically distinct from all Australian mainland individuals genotyped. As such, we were unable to identify the mainland Australian source for this group. These genome-level observations from the Hawkes Bay region mirror genetic patterns found more broadly in New Zealand possums using allozyme electrophoresis (Triggs, Green 1989). In that study, populations consisting predominantly of black possums were shown to be genetically similar to possums in Tasmania, while predominantly grey populations were genetically closer to Victorian and New South Wales possums.
There are several possible explanations for the genetic differences observed between the putative T.v.vulpecula animals in Hawkes Bay and the T.v.vulpecula animals sampled from mainland Australia. One possibility is that hybridisation between individuals of mainland Australian and Tasmanian origin masks the genetic signature of the mainland origins of New Zealand possums. Alternatively, the possums that derive from mainland Australia may be descended from a mixture of animals from several mainland sources, that were brought together in Hawkes Bay, and therefore do not resemble any single extant Australian population. Another possibility is that genetic drift-through founder effects and 80 years of isolation and/or selection imposed by a New Zealand environment-may have altered the genetic makeup of the Hawkes Bay population of T.v.vulpecula sufficiently to diverge from mainland Australian possums. This may even have resulted in the retention of genetic variants that are now rare (or even lost) in mainland Australia-a phenomenon seen in stoats introduced to New Zealand from England (Veale et al. 2015). Each of these three factors may contribute to our observations, although the fact that possums from sub-population A are genetically more similar to T.v.fuliginosus and to sub-populations B and C (Fig. 4), than they are to T.v.vulpecula, suggests that there has been leakage of T.v.fuliginosus alleles into sub-population A. A more detailed genomic haplotypic study, combined with mitochondrial DNA analysis, may resolve this issue. Additional sampling closer to the original release site for T.v.vulpecula in the Hawkes Bay may also help to identify the mainland location of origin for these possums.
Our results provide strong evidence that animals in sub-population B are a hybrid form of T.v. vulpecula and T.v. fuliginosus-the product of interbreeding between sub-populations A and C. The formation of a stable contact zone between two colonisers can arise when competition between the two is stronger than intraspecific competition or when hybrids are inviable or less fit than the parent taxa (Goldberg, Lande 2007). To our knowledge, interactions between the two subspecies or the viability of their offspring in nature have not been reported. We cannot therefore distinguish between these two possibilities. The positioning of the contact zone between two colonisers is likely to be attracted to, but not necessarily centred upon, a partial dispersal barrier (Goldberg, Lande 2007). In the case of the possums of Hawkes Bay, the presence of the Mohaka River between the release sites of the mainland Australian origin T.v.vulpecula and the Tasmanian origin T.v. fuliginosus provides such a partial barrier, and may have influenced the positioning of the contact zone near the eastern edge of this river.
The distribution of private alleles among the possum sites suggests that there is more immigration into the contact zone than out of it, from both the T.v.fuliginosus dominated population in the east and north of the zone, and from the T.v.vulpecula dominated population to the south west. The net result is that very few individuals with predominantly T.v.vulpecula ancestry (sub-population A) migrate completely through the contact zone (sub-population B) and integrate into the predominantly T.v.fuliginosus population (sub-population C), or vice versa. The contact zone thus appears to act as a sink for the two subspecies and, as such, represents a nontopographical (invisible) barrier to dispersal.
The existence of non-topographical or invisible barriers have implications for management, in this case particularly for the national program aimed at a predator free New Zealand (Russell et al. 2015). Populations with a mix of both grey (T.v.vulpecula) and black (T.v.fuliginosus) individuals occur in many parts of New Zealand (Cowan 2001;Kean 1971;Triggs, Green 1989) and hybrid animals have been reported previously in the Orongorongo valley on the North Island (Kean 1971). It is therefore possible that many such contact zones occur on the main islands of New Zealand. These zones could inhibit local dispersal and provide barriers that were previously invisible, but which could now be exploited for control. Conversely, possums can disperse quickly into depopulated habitats (Efford et al. 2000;Ji et al. 2001), so poisoning or trapping across contact zones could inadvertently break down existing subspecific barriers and increase local rates of dispersal. Documentation of the existence of additional contact zones on the main islands of New Zealand would provide greater opportunities to take advantage of these invisible barriers, with potential to produce faster localised eradication results.
Possums are the principal wildlife host of bovine tuberculosis in New Zealand (Morgan 1990), and the main cause of Tb persistence in cattle and deer herds in New Zealand, although the disease has not been recorded in Australian possums (Nugent, Cousins 2014;Tweddle, Livingstone 1994). If hybrid possums are less fit than either of the parental subspecies then the risk that possums are susceptible to bovine tuberculosis may be exacerbated. Moreover, there is some evidence that black coated animals, presumably T.v.fuliginosus, exhibit higher resistance to 1080 poison (Henderson et al. 1999;McIlroy 1982) than T.v.vulpecula-a phenomenon that has previously been hypothesised for marsupials more generally (Deakin et al. 2013). Thus, baiting with 1080 (a common approach used to control possums in New Zealand) could affect disproportionately the viability of one subspecies over the other, or even the viability of certain hybrids over others. Further research in this area, particularly through genomic analyses (Eldridge et al. 2020) will be particularly important if total eradication of possums is to be a key management objective.
In conclusion, we have identified clear evidence of a contact zone or border between two subspecies of introduced brushtail possums in Hawkes Bay, New Zealand, that produces a hybrid form. More work is required to determine whether this zone has formed because of inter-sub-specific competition or reduced hybrid fitness. The contact zone functions as a sink, with an asymmetry of gene flow tending towards the contact zone and away from both subspecies. It is therefore likely that the zone functions as a barrier to dispersal between the subspecies. Greater sampling would be required to determine the extent of the contact zone. Given the widespread overlapping distribution of the two subspecies in New Zealand, it is possible that there are many such hybrid zones that could function as barriers to dispersal across New Zealand. These hybrid zones could provide opportunities for managers to adopt a more nuanced approach to control of this introduced pest species, by accounting for differences between the subspecies and by taking advantage of previously unrecognised borders at contact zones. the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.