Effects of landscape matrix on population connectivity of an arboreal mammal, Petaurus breviceps.

Ongoing habitat loss and fragmentation is considered a threat to biodiversity as it can create small, isolated populations that are at increased risk of extinction. Tree-dependent species are predicted to be highly sensitive to forest and woodland loss and fragmentation, but few studies have tested the influence of different types of landscape matrix on gene flow and population structure of arboreal species. Here, we examine the effects of landscape matrix on population structure of the sugar glider (Petaurus breviceps) in a fragmented landscape in southeastern South Australia. We collected 250 individuals across 12 native Eucalyptus forest remnants surrounded by cleared agricultural land or exotic Pinus radiata plantations and a large continuous eucalypt forest. Fifteen microsatellite loci were genotyped and analyzed to infer levels of population differentiation and dispersal. Genetic differentiation among most forest patches was evident. We found evidence for female philopatry and restricted dispersal distances for females relative to males, suggesting there is male-biased dispersal. Among the environmental variables, spatial variables including geographic location, minimum distance to neighboring patch, and degree of isolation were the most important in explaining genetic variation. The permeability of a cleared agricultural matrix to dispersing gliders was significantly higher than that of a pine matrix, with the gliders dispersing shorter distances across the latter. Our results added to previous findings for other species of restricted dispersal and connectivity due to habitat fragmentation in the same region, providing valuable information for the development of strategies to improve the connectivity of populations in the future.

Introduction

The loss or fragmentation of natural habitats caused by human activities is one of the major threats to long‐term persistence of many species (Fahrig 2003; Foley et al. 2005). Remaining habitat patches are often small and isolated from each other by a matrix that may not be hospitable to resident fauna, such as agricultural and rural lands, plantations, settlements, and roads (Wilcove et al. 1986). The disruption to habitat continuity can lead to a reduction of population connectivity by preventing or reducing dispersal (e.g., Coulon et al. 2004; Banks et al. 2005; Lancaster et al. 2011), and potentially increasing the level of inbreeding and genetic drift within small isolated populations (Frankham 2005). The latter population processes result in a loss of genetic variability that provides the raw material for evolutionary change and is therefore crucial to the long‐term viability of isolated populations (Soule 1980; Reed and Frankham 2003; Frankham 2010).

The nature of the matrix surrounding remnant populations and the ability of individuals and species to use the matrix are important factors that influence the degree of connectivity across the broader landscape. If the matrix provides ecological requirements such as food and shelter, it may be utilized by a species (Kramer‐Schadt et al. 2004; Soule′ et al. 2004). Some matrices do not provide ecological requirements for species residency, but may be suitable for dispersal (Soule′ et al. 2004). Other matrices with no ecological value for a species may act as barriers to dispersal and prevent gene flow between neighboring populations. Before making general conclusions for conservation management of species within fragmented landscapes, individual species characteristics, including the degree of habitat specialization, dispersal potential through different matrix types, and behavioral responses to habitat fragmentation need to be considered (Weins 1997; O'Grady et al. 2004; Viveiros de Castro and Fernandez 2004; Meyer et al. 2009). Even species with similar life‐history characteristics may respond differently to the landscape matrix (Lindenmayer et al. 1999; Callens et al. 2011; Amos et al. 2014). As a consequence, we require species‐specific research to identify the influence of different surrounding matrices on dispersal and population connectivity (Debinski and Holt 2000; Cushman 2006; Callens et al. 2011).

Arboreal marsupials have long been recognized as a group of mammals potentially vulnerable to forest and woodland loss and fragmentation (McIllroy 1978; Bennett et al. 1991; Wormington et al. 2002; Laurance and Vasconcelos 2004; Lancaster et al. 2011; Taylor et al. 2011) due to their dependence on trees for nesting, foraging, and dispersal. Ecological and genetic consequences of fragmentation on arboreal marsupials have been infrequently reported (Taylor et al. 2007, 2011; Lancaster et al. 2011; Goldingay et al. 2013). In particular, few studies have tested the impact of fragmentation and the influence of different types of landscape matrix on gene flow and population structure in this group (e.g., Taylor et al. 2007; Lancaster et al. 2011; Goldingay et al. 2013).

In order to define the effect of landscape context on an arboreal species, traditional field methods such as radio‐tracking have been used (Dooley and Bowers 1998; Lindenmayer et al. 1999, 2000; Bladon et al. 2002). These methods are time‐consuming and limited in their ability to describe dispersal patterns at broad geographic scales. Without long‐term monitoring, it is difficult to infer the influence of different matrices on population connectivity and predict the viability of a species. In contrast, genetic markers can be used effectively to elucidate landscape‐scale dispersal and gene flow (e.g., Keller and Largiader 2003; Banks et al. 2005; Berry et al. 2005; Stow and Briscoe 2005; Taylor et al. 2007; Callens et al. 2011; Amos et al. 2014), providing indirect information about the response of a species to the surrounding matrix and the likely impacts of habitat fragmentation on population connectivity and persistence. Such information has important implications for conservation planning and understanding landscape effects on population structure.

Here, we use microsatellite markers to examine the effect of landscape matrix on population connectivity of an arboreal marsupial, the sugar glider Petaurus breviceps. Its distribution in Australia is confined to forests and woodlands of eastern and southern Australia. These forests and woodlands have been experienced widespread clearing since European settlement 200 years ago (Woodgate and Black 1988), resulting in species extinction and decline in this region (Woinarski et al. 2015). The remaining remnant native forests in some regions are now patchily distributed across the landscape and are isolated from one another by cleared agricultural land or pine (Pinus radiata) plantations. Petaurus breviceps is one of nine species of arboreal marsupials occurring in southeastern South Australia (Carthew 2004). Ecological connectivity for another arboreal species, the common ringtail possum Pseudocheirus peregrinus, in this fragmented landscape was investigated within seven patches surrounded by pine plantation, revealing that pine significantly impeded gene flow within the species compared to native forest (Lancaster et al. 2011). Although pine did not completely prevent movement of ringtail possums across the landscape, some consequences of isolation such as lower heterozygosity and genetic drift were evident in small patches. Pseudocheirus peregrinus is a generalist species that is not restricted to native forests for foraging and nesting (Lindenmayer et al. 2008). In contrast, P. breviceps is somewhat more specialized, in that it is not known to enter pine plantations, is more dependent on trees (including hollows for nesting) and is less likely to venture far along the ground for dispersal (Gibbons and Lindenmayer 2002). It is, therefore, expected to be more vulnerable to loss of habitat or replacement of native habitat with exotic plantations.

The effect of land clearing on the connectivity of populations of P. breviceps is difficult to predict as it is likely to be influenced by the density of eucalyptus trees that are remaining and the surrounding land use (e.g., cropped land or pastoral). Le Duff (2000) provided evidence that sugar gliders were more likely to be detected in patches which were surrounded by cleared agricultural lands rather than Pinus radiata plantations, suggesting that scattered large and old eucalyptus trees within cleared lands may assist P. breviceps to disperse between patches (Le Duff 2000). In a separate study of the arboreal marsupial P. peregrinus, in a nearby region in western Victoria, Lancaster et al. (In press) showed limited dispersal across an agricultural matrix compared to that within a continuous forest. Significant differentiation of populations and loss of genetic diversity within P. peregrinus were evident, and we predict that P. breviceps is likely to show similar patterns of population differentiation resulting from the agricultural matrix, albeit its ability to glide between trees over distances of ~90 m (Menkhorst and Knight 2004) may enable a higher level of dispersal compared to P. peregrinus.

This aimed to evaluate the effects of habitat fragmentation on population connectivity and genetic structure of P. breviceps in southeastern South Australia and explore mechanisms (such as sex‐biased dispersal and inbreeding avoidance) which may help sugar gliders maintain genetic diversity in this landscape. Specifically, this study aimed to (1) determine genetic diversity and population structure of the species, (2) test the influence of landscape features such as patch size, distance to nearest neighbor and degree of isolation on genetic diversity, (3) test the hypothesis that there will be isolation and restricted gene flow between forest patches, and (4) determine whether two different surrounding matrices, Pinus radiata plantations and cleared agricultural land, have a different effect on the dispersal capability of the species. Given previous results for P. peregrinus, and the habitat specialization of P. breviceps, we predict that gene flow will be significantly impeded across the pine matrix compared to cleared agricultural land.

Materials and Methods

Study sites and sampling

Samples of P. breviceps were collected from 12 remnant patches of native forest in southeastern South Australia (37°30′S, 140°25′E to 38°00′S, 141°00′E), and two sites (2.2 km apart) within a large continuous forest (Rennick State Forest, 5000 ha) in the adjacent southwest Victoria (37°55′S 140°58′E). Sampled patches vary in size (43–2216 ha) and the distance between them (1.4–19.2 km, Table 1), and are isolated from each other for more than 30 years by a matrix of cleared agricultural land or pine plantations (Fig. 1). Samples were collected between 2004 and 2009, using nest boxes or trapping. Two transects of five nest boxes were located in each of the patches commencing 50–100 m from access tracks, and had been installed as part of another project (Richardson and Carthew 2004). The distance between each transect was 100–500 m depending on the size of the patch. Nest boxes were checked monthly to sample gliders.

Table 1

Characteristics of sample sites including patch size, distance to neighboring patch, degree of isolation (mean distance to the closest three patches), and summary statistics of genetic diversity for Petaurus breviceps in the patches. Values are number of samples (N), allelic diversity (AD), allelic richness (AR), private alleles (PA), observed (H O), expected (H e) heterozygosity, inbreeding coefficient (F IS), probability (P) of F IS values, and mean relatedness (r). Numbers are given as mean ± SE

Patch Name	Patch abbreviation	Area (ha)	Distance to neighbor (km)	Degree of isolation (km)	N	AD	AR	PA	H o	H e	F IS	P	r
Casterton	CR	43	2.4	9.5	11	5.6 ± 0.7	3.7 ± 0.5	5	0.54 ± 0.07	0.63 ± 0.05	0.033 ± 0.06	0.12	0.08
Bourne	BN	80	4.7	15.7	16	6.1 ± 0.7	4.1 ± 0.5	4	0.69 ± 0.06	0.69 ± 0.03	−0.003 ± 0.08	0.13	0.06
Paltridges	PL	116	1.6	4	25	8.7 ± 1.3	5.3 ± 0.8	8	0.71 ± 0.05	0.73 ± 0.04	0.028 ± 0.051	0.15	0.04
Topperweins	TP	117	3.3	5.8	19	8.0 ± 0.9	4.8 ± 0.6	5	0.74 ± 0.04	0.72 ± 0.04	−0.028 ± 0.03	0.23	0.05
Penola	PN	139	4.7	16.2	14	6.2 ± 0.8	3.8 ± 0.5	1	0.65 ± 0.06	0.65 ± 0.05	−0.006 ± 0.05	0.10	0.09
Snowgum	SG	194	3.6	15.8	9	5.6 ± 0.6	3.7 ± 0.4	3	0.67 ± 0.06	0.65 ± 0.05	−0.028 ± 0.05	0.25	0.10
The Heath	TH	204	1.6	2.5	21	6.4 ± 0.7	4.0 ± 0.5	3	0.71 ± 0.05	0.69 ± 0.04	−0.02 ± 0.03	0.30	0.05
Mt Meredith	MM	250	6.8	8.6	10	4.9 ± 0.6	3.2 ± 0.5	5	0.53 ± 0.07	0.61 ± 0.04	−0.002 ± 0.07	0.10	0.08
Grundy's Lane	GL	260	8.2	13.8	43	9.8 ± 1.2	5.0 ± 0.6	5	0.71 ± 0.045	0.75 ± 0.03	0.044 ± 0.05	0.45	0.03
Yangery	YG	286	1.4	4	15	7.3 ± 0.9	4.4 ± 0.6	2	0.65 ± 0.07	0.69 ± 0.05	0.077 ± 0.07	0.04	0.08
Deadmans Swamp	DMS	525	19.2	22.7	48	9.5 ± 1.6	5.4 ± 0.9	12	0.69 ± 0.04	0.72 ± 0.04	0.042 ± 0.04	0.23	0.02
Nangwarry	NG	2216	8.6	9.1	5	3.5 ± 0.2	2.9 ± 0.2	3	0.74 ± 0.06	0.62 ± 0.04	−0.2 ± 0.07	0.87	0.13
Rennick State Forest	RSF	5200	NA	NA	14	7.3 ± 0.9	5.0 ± 0.7	3	0.76 ± 0.05	0.73 ± 0.03	−0.05 ± 0.06	0.33	0.07
Site 1	Ren1	NA	2.5	NA	8	6.1 ± 0.6	4.5 ± 0.5	2	0.82 ± 0.04	0.75 ± 0.04	−0.07 ± 0.05	0.23	0.03
Site 2	Ren2	NA	2.5	NA	6	5.0 ± 0.5	3.8 ± 0.4	1	0.77 ± 0.04	0.68 ± 0.05	−0.03 ± 0.04	0.40	0.05

Figure 1

Study area in the southeastern of South Australia and western Victoria, Australia, with remnant native forests embedded in pine plantation (light gray) or cleared agricultural lands (white). Sugar gliders were sampled from 12 labeled remnant patches (see Table 1 for full names of patches) and a large block of continuous forest (RSF).

Trapping was conducted to augment nest box samples. Wire cage traps were baited with creamed honey and were installed 3–6 m above the ground on metal brackets nailed to the trunk of trees (as per Carthew et al. 1999). Trap trees were also sprayed with a mixture of honey water, around and above the trap, as an attractant. Traps were checked at midnight and first light, and any captured gliders were checked for reproductive condition, sexed, and weighed. Reproductive condition in conjunction with Suckling's (1984) weight and tooth wear categories was used to allocate captured sugar gliders to three groups: adults, subadults, and juveniles. Small samples of ear tissue were taken from each individual, placed in a sterile vial of 50:50 ethanol:saline, and stored at room temperature before processing in the laboratory. Sugar gliders were individually tagged with uniquely numbered metal fingerling ear tags to avoid re‐sampling. In total, 250 tissue samples were collected from 12 native patches and two sites within the continuous forest (Table 1). Despite our attempt to sample gliders at five sites within the continuous forest, trapping was successful at only two sites, possibly due to the presence of the dominant yellow‐bellied glider, Petaurus australis, in this area. Further, difficulties associated with trapping gliders in tall eucalypt forests may also have resulted in a low sample size from the continuous forest. No nest boxes were available in the continuous forest prior to this study. Nest boxes were installed in this area in 2006; however, no animals were ever found in residence at the time of inspections.

DNA extraction and genotyping

Nuclear DNA was extracted from skin biopsies using the Gentra Puregene Extraction Kit and the manufacturer's procedure (Gentra Systems). Individuals were screened at 15 microsatellite loci, of which 13 had been previously isolated from the species (Brown et al. 2004; Malekian et al. 2013) and two originated from the squirrel glider, Petaurus norfolcensis (Millis 2000). The forward primer of each locus was fluorescently labeled, and PCR amplicons were visualized and scored as described in Malekian et al. (2013) (see supplementary table).

Genetic diversity analyses

Tests for departure from Hardy–Weinberg equilibrium (HWE) and linkage disequilibrium (LD) between loci were performed using the program GENPOP 3.4 (Raymond and Rousset 1995). A probability test based on a Markov chain algorithm (Guo and Thompson 1992) with 10,000 dememorizations, 100 batches, and 5000 iterations was conducted for each combination of locus and patch sample. The resulting P values were adjusted for multiple tests via the sequential Bonferroni method (Rice 1989). Observed (H o) and expected heterozygosity (H e) (Nei 1978) for each patch sample were calculated using Arlequin v 3.5.1.3 (Excoffier et al. 2005). Allelic diversity (average number of alleles per locus), allelic richness (allelic diversity corrected for sample size), private alleles (corrected for sample size), and the inbreeding coefficient (FIS) were estimated for each forest patch, using the program FSTAT 2.9.3.2 (Goudet 2001). The significance of FIS values was tested by permuting the alleles within samples over all loci in each sample using the program FSTAT 2.9.3.2 (Goudet 2001), with 1000 permutations. Bonferroni correction was applied to the resulting P values.

Relatedness analyses

Relatedness analysis was carried out to investigate the relationships among all sampled individuals using the Queller and Goodnight (1989) relatedness estimator with the software package GENALEX6 (Peakall and Smouse 2006). The age composition of individuals nesting together was used to categorize nesting groups as a putative family nest group (a mix of adults and juveniles/subadults), an adult nest group, or a subadult nest group. The differences in average relatedness among these categories were tested using a single‐factor analysis of variance (ANOVA). On the basis of these relatedness coefficients and prior to population genetics analysis, we removed one individual of each identified pair of relatives (full siblings or parent–offspring) to avoid any bias from sampling family groups in nest boxes. This process reduced the number of samples included in population structure analyses to 220 across the study area.

Genetic differentiation and population structure

Differences in allelic richness and heterozygosity (H e) among forest patches were assessed using one‐way analysis of variance (ANOVA) and post hoc Tukey tests in SPSS version 16.0.

We compared allele frequencies from different years (2005 and 2008) at two sites with the largest number of samples (Grundy's Lane, n = 43; Deadmans Swamp, n = 48). No significant differences (P > 0.05) were observed between years for these two sites, so at each site we pooled samples across years for further population genetic analyses. To determine whether there were differences in genetic variability among forest patches and to assess whether habitat fragmentation may have contributed to genetic differentiation, a hierarchical analysis of molecular variance (AMOVA) (Excoffier et al. 1992) was conducted using GENALEX6 (Peakall and Smouse 2006). We calculated pairwise F ST between patches and within continuous forest sites using ARLEQUIN v. 3.5 (Excoffier and Lischer 2010). Separate AMOVAs were conducted for patches surrounded by either pine plantations or cleared agricultural lands. F ST may underestimate differentiation among sites; therefore, we also calculated the estimator of actual differentiation Jost's D (Jost 2008) using the software program SMOGD (Crawford 2010). Statistical significance of these values was estimated using bootstrap analyses with 500 pseudo‐replicates.

To further investigate population structure, we employed several Bayesian approaches to identify genetic clusters across the landscape. Recent reviews and comparative tests on the use of Bayesian clustering software have highlighted the advantages of concurrently employing multiple programs to verify the number of clusters (K) within a dataset (Latch et al. 2006; Chen et al. 2007; Guillot et al. 2009; Francois and Durand 2010). Therefore, we implemented three Bayesian clustering packages to estimate K across our study region. Two of these analyses, STRUCTURE v. 2.2 (Pritchard et al. 2000) and BAPS v. 5.2 (Corander et al. 2003), were performed to infer clusters based on genotypic data alone, whereas TESS v. 2.3 (Chen et al. 2007) incorporated both genotypic and spatial (geographic coordinates of sampling locations) data to infer the most appropriate value of K.

The parameters used for STRUCTURE were admixed ancestry and correlated allele frequencies, and burn‐in and run lengths of 100,000 and 500,000, respectively. In both programs, K was investigated from 1 to 14 with 5 iterations of each K, as each patch could potentially represent one distinct population (12 patches plus two sites in continuous forest). Data from the continuous forest (RSF) were also analyzed separately to identify any structure between the two sites that may have not been detected in the whole dataset analysis.

The optimal value of K was assessed using the method described by Evanno et al.(2005), where the highest ΔK score represents the optimal number of populations. STRUCTURE HARVESTER v. 0.6.8 (Earl and vonHoldt 2012) was used to employ the Evanno method (Evanno et al. 2005) to select K from STRUCTURE results. We used CLUMPP 1.1.2 (Jakobsson and Rosenberg 2007) to average the membership probabilities of the most likely K.

BAPS was similarly run with five iterations of each K from 1 to 14 using the admixture model based on mixture clustering of individuals with 10,000 iterations, and true clusters were considered to be those that comprised three or more individuals (Corander and Marttinen 2006). Combining genotypic and spatial data, TESS was performed with 10,000 sweeps and a burn‐in of 5000 and we set multiple Ks from 2 to 14 with five iterations of each. We selected the conditional autoregression (CAR) admixture model and used the program default values for parameter settings. The optimal K value was chosen based on the stabilized value of the deviance information criterion (DIC).

Effects of landscape features on genetic diversity

We performed spatial autocorrelations within GENALEX6 (Peakall and Smouse 2006), using 0.5 and 1 km distance class sizes to explore the spatial scale of genetic variation. We assessed the dataset as a whole, and used two subsequent tests to independently assess patches separated by cleared agricultural land or pine. To test for sex‐biased dispersal, spatial autocorrelation was tested for adult female (86) and male (85) datasets separately, including adults and subadults. The significance of these analyses was statistically tested, using 95% confidence intervals defined by 1000 random permutations.

Finally, to further explore the relationship between genetic structure and environmental and spatial variables, we conducted a redundancy analysis (RDA) with the package vegan 2.2 in R (Oksanen et al. 2013). RDA is an analog of multivariate linear regression, using matrices of dependent and independent (explanatory) variables, and seems to have greater power than Mantel tests where there are multivariate species–environment relationships (Legendre and Fortin 2010). We used allele frequencies as dependent variables and environmental parameters as independent variables. The environmental parameters included spatial (X and Y coordinates, minimum distance to neighboring patch and degree of isolation) and habitat (patch size and percentage of suitable habitat around each patch) variables. For each patch, the degree of isolation was measured as the mean distance to the closest three patches (edge to edge, km) and the percentage of suitable habitat was estimated in a buffer zone of 1 km around each patch (Mapelli and Kittlein 2009). We then partitioned the variance components of the RDA by running three models – two separate partial RDAs in which genetic variance was conditioned on spatial and habitat variables, respectively, and a full model with all explanatory variables (Gugger 2012). These analyses allowed assessment of how much of the total genetic variance was explained by spatial factors, how much was explained by environmental factors, and how much was due to the joint effect of both factors.

Identification of dispersal events

To investigate whether dispersal is occurring across the matrices, we attempted to identify dispersal events by conducting first‐generation migrant detection in GeneClass v. 2 (Piry et al. 2004). Tests were performed according to the Bayesian method of Rannala and Mountain (1997) using the Monte Carlo re‐sampling approach of Paetkau et al. (2004) with 10,000 simulated individuals and a significance level of 0.05. Due to the size of the study region, several forest patches were not sampled so we implemented the appropriate model (“L=home”) for migrant detection that assumes that not all possible source forest patches have been sampled.

Results

Genetic diversity, Hardy–Weinberg equilibrium, and linkage disequilibrium

We used genotypes of all 250 individuals collected from 12 native forest patches and two sites within a larger tract of continuous forest to assess genetic diversity and departure from HWE and LD. A total of 278 alleles were scored across all forest patches and 15 loci, with the number of alleles at a locus ranging from 4 to 47. The average population allelic diversity ranged from 3.46 to 9.87. Moderate to high heterozygosity was found within each of the forest patches, with a mean heterozygosity of 0.68 across all loci (Table 1).

One locus (Petb15) showed a significant departure from HWE (P = 0.000), after Bonferroni correction for multiple tests (α = 0.003), suggesting that null alleles may be present. This locus was therefore removed from further analyses. No significant linkage disequilibrium was observed for pairwise locus combinations, as found in previous analyses of these microsatellite loci (Malekian et al. 2013).

Relatedness and family structure

Patch relatedness coefficients averaged between 0.02 in Deadmans Swamp and 0.13 in Nangwarry (Table 1), indicating that sampling was not biased toward highly related individuals. Average relatedness values were also calculated for nesting mates from 27 putative family (a mix of adults and subadults/juveniles), 14 adult, and 15 subadult nest groups. On average, putative family groups with an average group size of 4.4 individuals showed higher relatedness values (0.25 ± 0.04) than either adult (0.02 ± 0.07) or subadult (−0.06 ± 0.07) nest groups, and this difference was significant (single‐factor ANOVA; F = 14.36, P = 0.000). Pairwise relatedness values of nesting adult males, adult females, and adults of the opposite sex were also obtained. Results showed that, on average, adult females within nest boxes had higher relatedness values (0.19 ± 0.10) than did adult male pairs (0.10 ± 0.01) and adults' opposite sex pairs. Adult males and females that shared a nest box showed, on average, the lowest relatedness values (0.001 ± 0.06). Adult females with high relatedness values (r ≥ 0.5) were found in the same nest box (n = 27) with offspring.

Overall, there was a significant difference in allelic richness across the 13 sites (12 patches plus one continuous forest) (F = 3.49, P = 0.001). However, a post hoc Tukey test showed that only five pairwise comparisons were significant. Significant differences were between the four forest patches with the largest number of samples (Deadmans Swamp, Grundy's Lane, Paltridges, and Rennick State forest) and the patches with a small number of samples (e.g., Nangwarry and Mt. Meredith). Due to the lack of genetic differentiation within the continuous forest (see F ST results below), we pooled samples from the two sites. Heterozygosity levels were also significantly different among all sites (F = 1.99, P = 0.045). Two patches – Mt. Meredith (P = 0.04) and Casterton (P = 0.02) – showed significantly lower heterozygosity. No significant relationship was found between allelic richness and either patch size (r = 0.098, P = 0.85), distance to neighbor (r = 0.30, P = 0.2), or the degree of isolation (r = 0.2, P = 0.8). The relationships between heterozygosity and two of the three landscape features including distance to neighbor (r = 0.628, P = 0.044) and the degree of isolation (r = 0.48, P = 0.04) were significant. In addition, significant relationships were obtained between the number of private alleles and distance to neighbor (r = 0.758, P = 0.03) and the degree of isolation (r = 0.59, P = 0.027).

Despite an overall HW equilibrium in forest patches, positive values of FIS were recorded for five of them, although none were significant (Table 1) after corrections were made for multiple tests (k = 5 and P = 0.01); other forest samples, including the large continuous forest (Rennick State Forest), had slightly negative FIS values. Linear regression revealed no relationship between FIS and the three landscape features: patch size (r = 0.2, P = 0.09), distance to neighbor (r = 0.34, P = 0.08), and degree of isolation (r = 0.17, P = 0.5).

Overall, while the level of genetic subdivision (F ST) among all forest patches was not high, it was significantly different from zero (F ST = 0.1035 ± 0.005; P = 0.001). A large proportion (84%) of the genetic variance was explained by variation within patches, with 16% of variation among forest patches. Pairwise F ST values for forest patches ranged from 0.052 to 0.206 and were mostly significant (P < 0.00), with three exceptions: Bourne (BN) and Penola (PN); Paltriges (PL) and Yangery (YG); and Paltridges (PL) and The Heath (TH) (Table 2). Jost's D values were similar to F ST, ranging from 0.071 to 0.252 (Table 2). Pairwise F ST between the two sites within the continuous forest was nonsignificant (F ST = 0.02, P = 0.1).

Table 2

Pairwise F ST (below diagonal) and Jost's D values (above diagonal) for 12 native patches and a continuous forest (following sequential Bonferroni correction). Significant values are donated by *

	BN	CR	DMS	GL	MM	NG	PL	PN	RSF	SG	TH	TP	YG
BN	–	0.163*	0.105*	0.112*	0.183*	0.171*	0.093*	0.075	0.107*	0.155*	0.165*	0.106*	0.103*
CR	0.141*	–	0.193*	0.141*	0.252*	0.238*	0.123*	0.157*	0.106*	0.132*	0.213*	0.179*	0.143*
DMS	0.093*	0.127*	–	0.103*	0.181*	0.145*	0.071*	0.113*	0.091*	0.110*	0.113*	0.98*	0.097*
GL	0.094*	0.124*	0.081*	–	0.147*	0.125*	0.107*	0.123*	0.098*	0.117*	0.132*	0.107*	0.099*
MM	0.156*	0.206*	0.127*	0.116*	–	0.227*	0.138*	0.183*	0.132*	0.127*	0.165*	0.149*	0.132*
NG	0.140*	0.198*	0.130*	0.119*	0.194*	–	0.156*	0.147*	0.192*	0.211*	0.189*	0.171*	0.182*
PL	0.078*	0.117*	0.065*	0.088*	0.122*	0.145*	–	0.106*	0.099*	0.102*	0.071	0.093*	0.079
PN	0.052	0.130*	0.087*	0.092*	0.154*	0.138*	0.087*	–	0.147*	0.148*	0.147*	0.132*	0.126*
RSF	0.091*	0.091*	0.080*	0.067*	0.110*	0.130*	0.087*	0.090*	–	0.082*	0.131*	0.102*	0.112*
SG	0.118*	0.107*	0.093*	0.095*	0.112*	0.178*	0.085*	0.124*	0.061*	–	0.141*	0.125*	0.118*
TH	0.114*	0.172*	0.095*	0.110*	0.134*	0.171*	0.056	0.126*	0.112*	0.126*	–	0.105*	0.092*
TP	0.094*	0.156*	0.076*	0.099*	0.126*	0.144*	0.066*	0.100*	0.089*	0.111*	0.086*	–	0.082*
YG	0.092*	0.130*	0.076*	0.087*	0.120*	0.153*	0.053	0.104*	0.092*	0.097*	0.066*	0.086*	–

The effect of surrounding land use on genetic differentiation was assessed by dividing patches into two groups: those surrounded by cleared agricultural land and those surrounded by pine plantation. Significant and similar F ST values were obtained among patches surrounded by pine plantation (F ST = 0.08, P = 0.001) and cleared agricultural land (F ST = 0.10, P = 0.001).

STRUCTURE analysis identified nine well‐defined clusters (K = 9, Fig. 2). There was very little admixture between the clusters, which was reflected in the mean probability of membership (Q) for individuals assigned to each cluster (ranging from 0.60 to 0.95; Table 3). Some patches, including DMS, TP, MM, GL, NG, and SG, were each associated with unique clusters. We used a hierarchical cluster approach to further examine cluster 4, which occurred in more than two patches. Cluster 4 subdivided into one cluster associated with two patches (Paltridges and The Heath) and another with Yangery. BAPS found an optimal partition of eight clusters (K = 8). Proportional membership of individuals to each cluster based on BAPS and STRUCTURE was similar with one exception: the assignment of gliders from Snowgum, Rennick, and Casterton to one cluster. Population assignment in TESS was generally concordant with STRUCTURE results, with the highest DIC support for KMAX = 9 (Fig. 2B). No genetic structuring was identified within the continuous forest (K = 1, data not shown), but note the small total sample size (n = 14) from the continuous forest.

Figure 2

Genetic structure of the sugar glider in 12 patches and one continuous forest in southeastern South Australia: proportional membership (Q) of each individual to genetic clusters identified by STRUCTURE (A) and TESS (B). Each vertical bar represents the genotype of one individual glider, and the proportion of color in each bar represents the probability of membership in the relevant cluster. Forest patches abbreviated as given in Table 1.

Table 3

Mean probability of membership (Q) of gliders to each of the nine genetic clusters identified by STRUCTURE. Clusters where more than 60% of possums were strongly assigned to the same genetic cluster are bolded. Forest patches abbreviated according to Table 1

Cluster patch	1	2	3	4	5	6	7	8	9
BN	0.019	0.027	0.048	0.031	0.007	0.004	0.007	0.843	0.012
CR	0.006	0.890	0.030	0.005	0.005	0.009	0.027	0.014	0.009
DMS	0.89	0.0071	0.011	0.032	0.016	0.005	0.006	0.011	0.017
GL	0.035	0.026	0.009	0.063	0.028	0.801	0.010	0.015	0.008
MM	0.006	0.005	0.008	0.015	0.919	0.012	0.007	0.021	0.005
NG	0.006	0.007	0.004	0.003	0.005	0.006	0.003	0.010	0.952
PL	0.024	0.011	0.242	0.603	0.021	0.011	0.014	0.060	0.009
PN	0.018	0.007	0.021	0.015	0.019	0.194	0.018	0.699	0.016
RSF	0.007	0.652	0.008	0.008	0.244	0.012	0.033	0.016	0.017
SG	0.015	0.146	0.027	0.031	0.042	0.005	0.721	0.005	0.005
TH	0.036	0.012	0.017	0.865	0.017	0.006	0.008	0.006	0.029
TP	0.032	0.012	0.768	0.057	0.015	0.012	0.006	0.034	0.060
YG	0.016	0.006	0.020	0.722	0.019	0.184	0.012	0.009	0.007

Effects of landscape features on genetic differentiation

A significant relationship between genetic variation and the spatial variables was revealed by the full model of RDA (P = 0.001, Fig. 3A). When the analysis was controlled for habitat variables, we also found a significant association between genetic variation and the spatial variables (P = 0.015, Fig. 3B). Comparing the full model with the partial models indicated that the spatial variables explained 74.1% of the total explainable genetic variance. Habitat variables explained 20.3% of the total variance and spatial and habitat variables had a joint effect of 5.6% on genetic variance. Among the three spatial variables, geographic location showed the longest vector along each RDA axis, explaining 46% of genetic variation (Fig. 3B).

Figure 3

RDAs showing the contribution of spatial and habitat components to genetic structure in Peturus breviceps, for the (A) full model and (B) partial model controlled for habitat variables. Open circles are allele frequencies of each patch displayed in the RDA space, and the vectors show how explainable variables fall along that RDA space and crosses are centroids of environmental variables. DN, minimum distance to neighboring patch; DI, degree of isolation; PS, patch size; PSH, percentage of suitable habitat.

Spatial autocorrelation analysis of the whole dataset revealed a significantly positive coefficient at the distance size of 0.5 km (Fig. 4A). The signal of significant spatial autocorrelation was retained when the analysis was performed only on the females (Fig. 4B). In contrast, analysis of the males did not show a significant autocorrelation at 0.5 km (Fig. 4C) and did not depart from the assumption of random distribution of genotypes. Similar correlograms were obtained when a longer distance class of 1 km was selected (data not shown). Separate analysis of the five native patches surrounded by cleared agricultural lands showed significantly positive correlation values up to 3 km, and a positive, but not significant, trend up to 5.5 km (Fig. 4D). For the seven sites separated by pine, however, a positive and significant correlation was observed within 1 km distance, but not beyond that (Fig. 4E).

Figure 4

Correlograms of genetic correlation (r) plotted as a function of distance (0.5 km distance classes). The five plots represent data from (A) the whole dataset, (B) females, (C) males, (D) forest patches surrounded by cleared agricultural land, and (E) patches fragmented by pine plantation. The permuted 95% confidence intervals (dashed lines) and bootstrapped 95% error bars are shown.

Identification of migrants

Nineteen first‐generation migrant dispersal events were detected across the study region with the source population identified with high confidence for nine of these (Table 4). For the remaining events, the most likely source population was identified as the same site in which the individuals were captured, suggesting that the true source population was not represented in our sampled populations, but likely to be located near the patch the individual was sampled in. Of the nine putative dispersal events, four occurred across cleared agricultural land and five across pine plantations. The average dispersal distance spanning pine was about 2.5 km, while in the cleared agricultural matrix the dispersal distance averaged about 6.9 km. Most of the identified migrants (n = 7) were male, and two were female. Three of the dispersal events were between the continuous forest (RSF) and fragmented patches (SG and CR, Table 4).

Table 4

Summary of dispersal events detected by GENECLASS showing sampled and origin populations. Distance between origin patch and sampled patch was measured as edge to edge straight‐line in ArcGIS 9.3. Dispersal events were determined with a significance threshold of P < 0.05. M = male, F = female, A = cleared agricultural land, P = pine plantation. For full patch names, see Table 1

Sampled patch	Likely origin patch	Sex	Main surrounding matrix	Distance between origin patch and sampled patch (kms)
BN	PN	M	A	4.7
MM	GL	M	A	8.2
PL	TH	M	P	1.4
RSF	CR.	F	P	2.4
SG	RSF	F	P	3.6
SG	RSF	M	P	3.6
TH	YG	M	P	1.5
YG	MM	M	A	7.5
YG	MM	M	A	7.5

Discussion

Maintenance of connectivity among populations in fragmented landscapes is important to mitigate against the detrimental effects of inbreeding in small populations, which can in turn reduce fitness levels and the potential of a species to respond to environmental changes (e.g., the introduction of a virus or new pathogen) (Frankham 2005, 2010). However, information about genetic connectivity in fragmented landscapes, and how different matrices influence this connectivity, is lacking for many species. Here, we used 14 microsatellite markers to investigate genetic structure and connectivity of sugar glider populations in a fragmented landscape of southeastern South Australia.

Fragmentation effects on sugar glider population connectivity

Evidence for genetic structure across the fragmented landscape

Significant genetic structuring across the fragmented landscape was revealed from population structure analysis, with samples grouped into at least eight or nine distinct population clusters, of which six were each associated with a single patch. It is not possible to infer from the current data the exact time frame over which this differentiation arose. However, the concurrent lack of any such clustering over similar geographic distances in the continuous eucalypt forest of Rennick State Forest (although data were limited) tends to suggest that gene flow among patch populations is limited, and/or the effective population sizes within each patch may have reduced, leading to genetic differentiation by genetic drift. Similar results were obtained for a related gliding marsupial, Petaurus norfolcensis, for which significant genetic differentiation of populations was reported in each of two fragmented landscapes in central and southern Queensland (Goldingay et al. 2013) and in some coastal populations fragmented due to urbanization and agricultural practices (Taylor et al. 2011). Frequency‐based analyses also showed strong differentiation among patches, with those located closer together being less differentiated than those further apart, suggesting that dispersal might be occurring within short distances and between neighboring patches. This tendency for short dispersal distances and between proximate populations was also observed in the common ringtail possum (Lancaster et al. 2011) and southern brown bandicoot (Li et al. 2014) in the same region, and it has been observed in several other mammalian species in fragmented landscapes (Goossens et al. 2005; Bergl and Vigilant 2007; Taylor et al. 2007; Fitzgibbon et al. 2011).

Influence of the surrounding matrix on dispersal and gene flow among patches

Autocorrelation analyses revealed positive spatial genetic structure for gliders across pine and cleared agricultural matrices over distances of 1.0 and 3 km, respectively. It is likely that the pine matrix does not provide suitable habitat requirements, such as food and shelter, for sugar gliders, so that they do not utilize these matrices extensively, restricting their movement across pine and further suggesting that the presence of pine reflects a loss of habitat for sugar gliders. For the agricultural matrix, the longer distance for genetic spatial autocorrelation compared to the pine matrix may reflect a greater utility of this matrix for food, shelter, and dispersal; however, we cannot rule out the possibility that the spatial autocorrelation patterns in the correlograms (see Fig. 4) resulted from sampling differences. Overall, it appears that pine plantations have restricted gene flow of sugar gliders across the landscape, and as such, native patches remaining in the landscape are likely to represent most of the habitat available to support populations of the species. These findings are in line with other studies of marsupials in fragmented landscapes where pine plantations have restricted connectivity among native forest patches (Banks et al. 2005; Peakall et al. 2006; Taylor et al. 2007).

The migrant analyses based on GENECLASS provided further support for the greater permeability of cleared agricultural land compared to pine plantation for sugar gliders. Putative dispersal events across pine plantations were largely restricted to neighboring patches, with an average migration distance of 2.5 km compared with that for cleared agricultural land of 6.9 km. These results suggest that scattered “paddock trees” or thin corridors along fence lines and roadsides may provide species requirements, such as foraging and/or nesting sites, or assist P. breviceps to disperse between patches. Use of large trees as stepping stones for movement was also evident in the dasyurid marsupial Antechinus flavipes, where animals were tracked to large trees in otherwise open paddocks (Marchesan and Carthew 2008). Gliders' ability to move through the landscape is likely to be restricted by the extent of tree cover (van der Ree et al. 2004; Taylor and Goldingay 2009). van der Ree et al. (2004) showed that Petaurus species were most likely to occur in remnant sites within 75‐m proximity, which corresponds with the maximum gliding distance in a single movement between trees. Sugar gliders can move between trees over distances up to 90 m in a single leap, depending on the size of trees (Menkhorst and Knight 2004). Gaps in tree cover that exceed the gliding distance threshold can therefore act as barriers to movement for gliders (Ball and Goldingay 2008). Dispersal movements of young gliders is known to occur along forested roadside vegetation for distances up to 1.9 km and across treeless gaps up to 200 m (Suckling 1984). Radio‐tracking studies on the squirrel glider showed that movements decrease as canopy gaps increase beyond 50 m (van der Ree et al. 2004, 2010). In another study, loss of intervening tree cover led to genetic differentiation of squirrel glider populations inhabiting the fragments (Goldingay et al. 2013).

How gliders may mitigate the effects of fragmentation

Evidence for male‐biased dispersal

Spatial autocorrelation analyses of males and females separately suggested that male sugar gliders were largely unrelated even in the smallest distance class, in contrast to females. Additionally, the relatedness analysis of nest‐sharing animals showed elevated relatedness among adult females within a nest group compared to male adults. Analyses of migrants also detected a larger number of male than female migrants. Previous inferences of dispersal in sugar gliders have been based on limited direct observations and video camera data, with somewhat conflicting results. Although both male and female young dispersed from their natal population in eucalypt remnants in an agricultural matrix in southern Victoria (Suckling 1984), Sadler and Ward (1999) reported that nesting groups of sugar gliders had female‐biased sex ratios and inferred that the dispersing sex was male. Overall, our results are consistent with the species showing female philopatry and male‐biased dispersal.

Male‐biased dispersal may help gliders avoid inbreeding, an evolutionary process that might be perturbed by habitat fragmentation. Elevated relatedness among individuals inhabiting a remnant – presumably due to increased inbreeding – was evident in genetically isolated populations of squirrel gliders in fragmented landscapes (Goldingay et al. 2013). This effect was also observed in three lizard species and a bird over fragmented landscapes (Delaney et al. 2010). No significant positive inbreeding coefficients were recorded in this study, suggesting gliders may employ inbreeding avoidance mechanisms.

Evidence for inbreeding avoidance

Within nest boxes, adults of the opposite sex were less related than were pairs of adult males or females. The presence of unrelated males and females as potential sexual partners within nest boxes could result from a natural tendency of the species to choose unrelated partners to avoid inbreeding and may be a consequence of sex‐biased dispersal. The lack of significantly positive FIS values at the patch level can also be interpreted as a sign of inbreeding avoidance.

In a fragmented landscape with reduced capacity for dispersal, kin recognition may act as a mechanism for inbreeding avoidance, thereby mitigating the problem of related individuals breeding in small populations. Mechanisms of social recognition in natural populations of P. breviceps remain unclear (Mallick et al. 1994; Klettenheimer et al. 1997; Sadler and Ward 1999). Male P. breviceps possess scent glands (e.g., frontal, sternal, and urogenital) that secrete pheromones, and Schultze‐Westrum (1965, 1969; cited in Suckling 1984) suggested that scent marking is important in determining the social organization of sugar glider captive groups. Schultze‐Westrum (1969) suggested that pheromones may be transferred to group members by one or two of the dominant males. These dominant males also perform most of the other social activities such as mating, territory maintenance, territory patrolling, and aggression against outside individuals. However, little is known about scent markings in natural populations and the role of pheromones in the social structure of P. breviceps.

Evidence for family structure within nest boxes

The analysis of genetic relatedness between individuals within nest boxes showed that groups of gliders sharing nest boxes were generally comprised of related individuals, suggesting P. breviceps preferentially shared nests with kin. The occurrence of co‐nesting by related adult females may imply that they live and rear their offspring together. Adult females with high genetic relatedness were found in the same nest box with juveniles. Although data are limited, the presence of related females with juveniles in the same box may imply cooperative rearing of offspring of P. breviceps. Potential benefits of this behavior include protection of offspring from infanticide, improved thermoregulation, and adoption of young whose mother dies (Hayes 2000).

Conclusion

Within the fragmented landscape of southeastern South Australia, research on multiple species with different life strategies (sugar glider, common ringtail possum, and southern brown bandicoot has provided evidence of restricted dispersal and connectivity of habitat patches due to fragmentation. This raises concerns about the long‐term viability of native mammal species in this region and further suggests there is a need to develop long‐term conservation management plans that mitigate the effects of fragmentation. A strategy to develop corridors between native forest patches is being implemented (Horn 2003), and our analyses here will provide an important basis for assessing whether this strategy has been successful in improving gene flow among patches in the future.

Conflict of Interest

None declared.

Table S1. Fifteen developed polymorphic loci isolated from Petaurus breviceps and P. norfolcensis including locus name, primer sequence (F, forward; R, reverse complement), core repeat motif, GenBank Accession numbers and the reference.

Click here for additional data file.