Genetic Structure and Gene Flow in Eastern Grey Kangaroos in an Isolated Conservation Reserve

: Dispersal is a key process for population persistence, particularly in fragmented landscapes. Connectivity between habitat fragments can be easily estimated by quantifying gene ﬂow among subpopulations. However, the focus in ecological research has been on endangered species, typically excluding species that are not of current conservation concern. Consequently, our current understanding of the behaviour and persistence of many species is incomplete. A case in point is the eastern grey kangaroo ( Macropus giganteus ), an Australian herbivore that is subjected to considerable harvesting and population control efforts. In this study, we used non-invasive genetic sampling of eastern grey kangaroos within and outside of the Mourachan Conservation Property to assess functional connectivity. In total, we genotyped 232 samples collected from 17 locations at 20 microsatellite loci. The clustering algorithm indicated the presence of two clusters, with some overlap between the groups within and outside of the reserve. This genetic assessment should be repeated in 10–15 years to observe changes in population structure and gene ﬂow over time, monitoring the potential impact of the planned exclusion fencing around the reserve.


Introduction
Habitat fragmentation due to land use change or exclusion fencing, or both, is a significant threat to biodiversity [1,2]. Isolated populations are prone to inbreeding, which, in turn, negatively affects reproductive fitness, the ability to adapt, and survival rate [3]. Connectivity and movement among wildlife populations and subpopulations are therefore essential for the long-term persistence of many species [4].
Currently, most of our knowledge about the movement of animals comes from observational data (however, these are often sparse and might provide poor estimates of population status and dynamics [5]), or else, from tracking technologies that require the capture of individuals and the process of fitting them with devices [6,7]. In contrast, several non-invasive landscape genetics methods can be employed to assess gene flow among populations or individuals, without the need to disturb the animals [8][9][10]. Furthermore, the genetic diversity among groups can be used to detect relatedness, population structure, inbreeding, and effective population size [11,12]. Estimation of these demographic parameters is important for the design and evaluation of conservation measures [13].
In Australia, the eastern grey kangaroo (Macropus giganteus) has been subjected to anthropogenic pressure since the first European settlement in 1788. Kangaroos were initially regarded as a sustainable food source [14], but towards the middle of the 19th century, kangaroos started to become "harvested" in great numbers [15]. Despite the scientific evidence to the contrary, kangaroos are still often perceived to compete with livestock by landholders and are subjected to considerable eradication efforts and a commercial killing industry [16][17][18]. Every year, approximately 3-4 million kangaroos are killed in Australia [19,20]. Consequently, eastern grey kangaroo populations in Queensland and other parts of Australia are disfavoured in agricultural settings but can find respite in conservation reserves. Therefore, it is crucial to monitor and assess the trajectories of kangaroo populations to understand the role of these conservation reserves on their long-term persistence.
The present study aimed to assess the fine-scale genetic variation and gene flow in a heavily fragmented landscape that contained a conservation reserve, the Mourachan Conservation Property, in which kangaroos are not persecuted. Recently, construction work of an exclusion fencing around the reserve commenced ( Figure 1). Exclusion fencing has been used in Australia for more than 100 years to exclude "unwanted" species, such as rabbits, dingoes, or emus [21,22]. Such a barrier will likely impact the gene flow among populations of the eastern grey kangaroo and other species within and outside of Mourachan, as well as influence the relationships between individuals and species within the reserve. Knowledge of the current genetic diversity and functional connectivity will help to establish a baseline for future genetic assessments. This information can be used both for increasing our understanding of eastern grey kangaroo ecology and for providing important tools for management authorities for the evaluation of the impact of exclusion fencing.  Table 1 for details on location ID). (C) Average assignment probabilities from STRUCTURE for K = 2 are represented by pie charts, proportionally sized to the number of samples: cluster 1 in blue, cluster 2 in pink. (D) Individual assignment coefficients are displayed in the bar plots (each individual is represented by a single vertical bar, with sampling location IDs shown on top. The map sections were produced using QGIS 3.16.4 [23]. Table 1. Locations ID and names; latitude; longitude; position (inside or outside of the Mourachan Conservation Property); number of samples (N); allelic richness (A R ); observed heterozygosity (H O ); expected heterozygosity (H E ); inbreeding coefficient (F IS ; statistically significant values in bold); and STRUCTURE Q-values assigning the population to cluster 1 (Q 1 ) or cluster 2 (Q 2 ; see Figure 1 for details).

Study Species
The eastern grey kangaroo (Macropus giganteus) is a grey-brown macropod reaching over 70 kg in weight, with a body length of up to 2.3 m in males and 1.9 in females [24,25]. Its distribution spreads from eastern Australia to south-eastern South Australia and northeastern Tasmania [26]. There have been significant differences in dispersal, movement, and home-range size patterns reported across different study areas of eastern grey kangaroo populations [27,28]. For instance, home-range size was recorded to range from 7.6 to 269 ha for male kangaroos [28]. Eastern grey kangaroos live in small communities, known as mobs, with a fission-fusion social structure [24,29]. They are crepuscular animals, with peak activity approximately two hours before sunset, retiring to shelter under shrubs or trees to rest two to three hours after sunrise [24].

Study Area
The Mourachan Conservation Property is a 48,000 ha reserve situated on the Brigalow Belt in Queensland, Australia, spanning approximately between 27 • 50 13.9 S 149 • 02 30.8 E and 27 • 37 56.1 S 149 • 06 05.2 E (Figure 1). It is one of the last conservation properties in the region, representing a stronghold of recovered vegetation and wildlife, despite a long history of clearing and agriculture [30]. The area around Mourachan is dominated by cotton and livestock farming, where large wildlife fauna are excluded. According to the SPEI drought index [31], the Brigalow Belt, along with much of Australia, experienced increasingly moderate-severe drought conditions after 2016, peaking as extreme conditions in the summer of 2018 and in late 2019 [32], when sampling took place.

Sample Collection
In this study, strictly non-invasive genetic sampling was implemented [33]. We collected eastern grey kangaroo scats at 17 sites in and around the Mourachan Conservation Reserve in November 2019 (Table 1 and Figure 1). To increase sampling efficiency, we used a targeted sampling approach of areas where the presence of kangaroos was known, i.e., around artificial dams or grass patches.
Samples were collected very early in the morning (ca. 5-7 a.m.) and late afternoon and night (ca. 5-10 p.m.), reflecting kangaroos' crepuscular activity. We placed 1-3 pellets from each fresh faeces pile in a plastic bag and stored them in a freezer at −20 • C.

DNA Extraction and Genotyping
We extracted genomic DNA from 246 scats using QIAamp Fast DNA Stool Mini Kit (Qiagen Pty Ltd, Doncaster, Australia) following the manufacturer's instructions, with a final DNA-elution in 50 µL of nuclease-free water (Thermo Fisher Scientific Pty Ltd, Scoresby, Australia). DNA concentrations were measured using a NANODROP 2000 (Thermo Fisher Scientific Pty Ltd, Scoresby, Australia) and the samples were then sent to the Australian Genome Research Facility Ltd. for sequencing. The samples were genotyped at 20 microsatellite markers ( Table 2) in three multiplex panels, using internal protocols. A repeat round of genotyping was performed if a sample failed for a given marker in the first round. Amplification products were sized by comparison to a GENESCAN™ 500 ROX™ Size Standard using the software GENEMAPPER 3.7 (Thermo Fisher Scientific Pty Ltd, Scoresby, Australia).

Data Analysis
The genetic diversity was calculated as the number of alleles per locus (N A ), observed (H O ) and expected (H E ) heterozygosity, potential deviations from Hardy-Weinberg equilibrium (HWE), and null alleles in CERVUS 3.0.7 [34]. We also used CERVUS to identify matching genotypes across all samples to exclude duplication of the same individuals. Allelic richness (A R ) and inbreeding coefficient (F IS ) were quantified in FSTAT 2.9.4 [35], H O , and H E per sampling location in GENALEX 6.5 [36]. Genetic differentiation between the pooled sampling locations inside and outside of the Mourachan Conservation Property (see Table 1 for details) was calculated as the pairwise genetic distance (F ST ) in FSTAT 2.9.4.
To assess the most likely population genetic structure of the eastern grey kangaroos within and around the Mourachan Conservation Property, we used the clustering algorithms implemented in STRUCTURE 2.3.4 [37]. STRUCTURE uses a Bayesian clustering method to assign individuals to one of K genetic clusters and to estimate the degree of admixture among these clusters. The settings were as follows: the number of genetic clusters (K) 1-10, 10 iterations for each K value, 10 5 burn-in and 10 5 Markov chain Monte Carlo replicates in each run (as recommended by [38]), using sampling locations as prior and with admixed and correlated allele frequencies. CLUMPP 1.1.2 [39] was used to re-evaluate the population assignment and generate an average output, which was plotted using DIS-TRUCT 1.1 [40]. To summarize major variation patterns in the dataset, we also conducted a principal coordinate analysis (PCoA) in GENALEX 6.5. To assess whether the populations underwent a recent genetic bottleneck, we compared the distribution of heterozygotes for each marker with the number of alleles in BOTTLENECK 1.2.02 [41]. We used the two-phase substitution model recommended for microsatellite markers, with default proportions of the infinite alleles and a single stepwise mutation model. The one-tailed Wilcoxon signed-rank test [42] was implemented to estimate the likelihood that the observed data deviated from what is expected under mutation-drift equilibrium. Since a minimum of thirty individuals are recommended to achieve power of at least 0.8 with the Wilcoxon signed-rank test [41], we pooled the sampling locations into two groups: one group inside and one group outside of the Mourachan Conservation Property (see Table 1 for details) for this analysis.
We estimated the current genetic effective population size (N e ) using the linkage disequilibrium (LDNE) method [43], implemented in NEESTIMATOR 2.1 [44], excluding singleton alleles (i.e., alleles that occur in one copy in one heterozygote, which can upward bias N e estimation). NEESTIMATOR estimates genetic effective population sizes (N e ) from genotype data and the recently improved method implemented in the new version reduces bias in estimating N e and increases reliability of confidence intervals. Table 2. Indices of the 20 microsatellite markers used in this study, genotyped across 232 samples (see Table 1 for more details).

Genetic Diversity
Out of the originally sampled and genotyped 246 samples, 14 were removed due to duplication, leaving 232 unique genotypes ( Table 1). The mean proportion of loci genotyped across the 232 samples was 95.26%. The number of alleles per locus (N A ) ranged from 13 (locus B123) to 57 (locus Me28), with an average of 26 (Table 2). H O across all 232 genotypes ranged from 0.526 at locus G16-2 to 0.916 at locus Me28 (Table 2). Across all loci and individuals, seven loci deviated significantly from HWE, probably due to the potential presence of null alleles (Table 2).
Overall, genetic diversity was relatively high across all sampling locations, with the average allelic richness (A R ) of 5.023 (Table 1) and observed heterozygosity (H O ) ranging from 0.631 (1000 Acre Dam) to 0.822 (End Dam; Table 1). F IS estimates, ranging from −0.043 (Nr 11 Dam) to 0.240 (Bob's Billy Dam) were statistically significant in 10 of the 17 sampling locations (Table 1). There was very low genetic differentiation between the pooled locations inside and outside of the reserve (F ST = 0.029).
We did not detect any evidence of genetic bottlenecks in the pooled inside (p = 0.559) or outside (p = 0.563) sampling locations. The effective population size N e for the inside subpopulation was estimated to be 47 (95% CI [46][47][48], for the outside subpopulation to be 40 (95% CI 38-42).

Genetic Structure and Gene Flow
The analysis with STRUCTURE showed high consistency among different runs and detected two clusters based on the mean ln P (k) and Delta K, determined by STRUCTURE HARVESTER (Figure 1). A slight trend in the grouping of individuals from sampling locations could be observed: individuals from all but one location outside of the reserve clustered together (Figure 1).
The PCoA showed three separate clusters: two with mostly individuals from inside the Mourachan Conservation Property and a third one with a mixture of individuals from inside and outside the property (Figure 2). Individuals from the same location clustered predominantly together, but different sampling locations overlapped. The first and second axis accounted for only 12.91% and 3.82% of the variation, respectively.  Table 1 and Figure 1 for details on location IDs). The percentage of the total variation in the dataset explained by each principal coordinate is given in parentheses.

Discussion
In this study, we implemented population genetics methods to assess the current genetic diversity gene flow among eastern grey kangaroo subpopulations within and outside of the Mourachan Conservation Property. The study was timed to coincide with the establishment of exclusion fencing around most of the property. The microsatellite markers proved to be highly polymorphic (Table 2), confirming the results of other studies [49]. We found relatively high levels of heterozygosity, similar to previous work on eastern grey kangaroos [27,45,48,49].
Despite the long history of harvesting and persecution of kangaroos outside of the reserve, we could not detect any sign of a bottleneck in the data. However, it is important to note that populations experiencing demographic decline may not necessarily suffer from a genetic bottleneck, and vice versa [50]. Only extreme cases of severe population decline will generate a detectable genetic footprint [51,52]. Kangaroos are known to move over large distances to pursue food and water resources [53]. The planned exclusion fencing will likely constitute a barrier to migration that might result in abrupt population decline, as already reported from other locations [54].
One of the most important concepts in conservation is the effective population size (N e ). Effective population size influences the rate of genetic variation loss, inbreeding, and gene flow [55]. We determined the N e to be only 47 for the subpopulation inside the reserve. This number is comparable to those reported for several populations of another macropod, Setonix brachyurus, in the study by Spencer et al. [56]. Effective population size is the size of an ideal population and is, in general, lower than the census population size. This might be the result of several factors, for instance, fluctuating population size [57], overlapping generations [58], or the pattern of spatial distribution of a population [59].
Individuals from the sampling locations within and outside of the Mourachan Conservation Property exhibited very limited differences in genetic diversity, with very low genetic differentiation (F ST ). Although the STRUCTURE and PCoA analyses revealed a trend in the individuals outside of the reserve clustering together (Figures 1 and 2), we cannot conclude that there has been a recent or current barrier to gene flow. Considering the large home range sizes reported for eastern grey kangaroos (7.6 ha-269 ha [28,49]), this finding is not surprising. Our results are also in line with previous studies indicating a limited genetic structure in kangaroo populations [28,60] and are analogous to patterns observed in other species with high dispersal ability such as ungulates [61] or carnivores [62]. Although eastern grey kangaroos are highly sedentary, some male-biased dispersal has been suggested [60,63] and may be sufficient to promote genetic mixing. Further investigations would be necessary to confirm this assumption.
Furthermore, because our results provide only a temporal snapshot, it is impossible to determine at this point whether the two clusters detected (Figure 1) are the beginning of a population differentiation process. We recommend expanding the sampling area to provide a more comprehensive overview of the species' genetic diversity, and repeated sampling over longer periods to capture changes over time. Exclusion fencing that has since been raised around the reserve will likely disrupt this connectivity: fences create a boundary that prevents migration of species [64]. It has been shown that even relatively recent barriers can significantly impact population structure. For instance, Epps et al. reported a marked decline in genetic diversity in the desert bighorn sheep population after the construction of highways within their habitat [65].
This study serves as a baseline of genetic diversity prior to the effects of the physical barrier can be felt; hence, it is crucial that this genetic assessment is conducted repeatedly every few years to assess the potential impact of the exclusion fencing. The knowledge of exclusion fencing effects can help in designing mitigation strategies that could be employed to preserve the natural processes and ecosystem services [66].  Table S1.