Female-biased introductions produce higher predicted population size and genetic diversity in simulations of a small, isolated tiger (Panthera tigris) population

Isolation of wildlife populations represents a key conservation challenge in the twenty-first century. This may necessitate consideration of translocations to ensure population viability. We investigated the potential population and genetic trajectory of a small, isolated tiger (Panthera tigris) population in Thailand’s Dong Phayayen-Khao Yai forest complex across a range of scenarios. Using an individual-based, spatially-explicit population modelling approach, we simulate population and genetic trajectories and evaluate the relative impact of translocations from a related population. Population and genetic trajectories in our study were most sensitive to sex and number of individuals translocated and translocation frequency. Translocation of females produced consistently higher population, allelic richness, and heterozygosity compared to equal numbers of males. Despite population increases, declines in allelic richness and heterozygosity across simulations were stark, with simulations predicting a mean decline of allelic richness and heterozygosity of 46.5% and 53.5% without intervention, respectively. Translocations of four females every generation or every other generation were required to prevent substantial heterozygosity declines. While translocations could increase population size, they may fail to prevent long-term loss of genetic diversity in small populations unless applied frequently. This reinforces the importance of incorporating realistic processes of genetic inheritance and gene flow in modelling small populations.


Supplementary Materials 1 -Model Overview
In this section, we provide a description of our individual-based model and approach via the OOD (Overview, Design concepts, Detail) protocol, as described by Grimm et al. 1,2 .

Purpose
The purpose of our modelling approach is to (1) quantify the likely population and genetic trajectory of the tiger population of the Dong Phayayen-Khao Yai Forest Complex (DPKY) in its current state and (2) to determine the relative impact of translocations of tigers from elsewhere in Thailand on these trajectories.
In addition, we aimed to evaluate the sensitivity of predictions to starting genetic variation (APL), genetic relatedness between source and destination populations (SA%), variation in the number and sex of individuals translocated into DPKY (TSEX), frequency of translocations (TFRQ), and mortality risk in translocated individuals (TMOR).

Entities, state variables, and scales
Our models utilize the program CDPOP (Cost Distance POPulations 3 ) which uses an individual-based, spatially-explicit framework to model mating, genetic inheritance, dispersal, and mortality through space and time. Individuals are characterized by several state variables (see Table 1), including: a unique alphanumeric identifier, sex (M/F), location (coordinates), genotype data for user specified number of loci and alleles, and territory unit. Females mate with a single male and males can mate with multiple females in a given timestep. The number of offspring is drawn probabilistically from a normal distribution (mean of 3 and standard deviation of 2, with a 1:1 sex-ratio (Thatte et al. 4 ; Tables 1 and 2). Generations are nonoverlapping and, therefore, we do not distinguish individuals by age class.
All individuals in our model are assigned unique genotypes based on 14 theoretical loci with an average number of alleles per locus determined by one of several potential values. The number of alleles for each locus was generated via a Poisson probability draw with a mean corresponding to either two (APL2), three (APL3), or four (APL4) mean alleles per locus. Allele frequency probabilities were determined by a generating random number between 0 and 1 for each allele, and by dividing the result by the sum of all random numbers generated at that locus. The final allele frequency probability value was then used by CDPOP to generate unique genotypes for each individual in the simulation (intgenesans function in CDPOP). The same process was repeated to simulate genotypes for translocated individuals from WEFCOM with 14 loci and an average of 6 alleles per locus, corresponding to values from Klinsawat 5 . Due to the unknown degree of genetic relatedness between DPKY and WEFCOM populations at identical loci, when generating alleles for each locus between populations, we varied the percentage of DPKY alleles present in simulated WEFCOM individuals (SA%) at three levels: 25% (SA25), 50% (SA50), and 75% (SA75).
Occupiable locations in our study extent were determined by resampling a cost-distance kernel of 45,000 cost-distance units, originating from DPKY's contiguous forest boundary, to a resolution of 10km. These territorial units were then converted to points representing potential tiger home range centroids at a theoretical population density of 1/100 km 2 . For simplicity, maximum carrying capacity density is assumed to be constant over time.
Spatial patterns of the study landscape, with which individuals in simulations would interact, were defined in a resistance surface. A resistance surface reflects step-wise cost to individual movement, with high pixel values representing high resistance to movement and low values imparting little resistance. In our study, resistance surfaces were developed at a 250m resolution. We utilized a base resistance surface developed by Ash et al. 6 and utilized in Ash et al. 7 as the most representative of likely patterns of resistance among four candidate models. This expert-based surface was based on 2018 SERVIR-Mekong Regional Land Cover Monitoring System (RLCMS 8 ) land cover data, classified into dense forest (resistance value of 1), scrub forest (resistance of 20), agriculture/village matrix (resistance of 50), reservoirs/surface water (resistance of 80), and urban areas (resistance of 100). Minor (resistance of 30) and major roads (resistance of 100; 9 ) were also included.
Individual movement is defined by sex-dependent movement parameters, differentiated for mating and dispersal behaviour. These and other simulation parameters in our models were based on other studies on tigers and other wide-ranging felids 4,10,11 . The CDPOP model simulates mating among individuals via probabilistic selection of mating pairs based on the cost distance between their locations across the resistance surface. Individuals that are in close proximity in cost distance are more likely to mate than those that are farther, with the probability defined by the mating distance function and maximum distance (see Table 1). Dispersal of offspring produced in mating is also probabilistic based on cost distance, with locations in close cost-proximity more likely to be selected than those that are far, governed by the specified dispersal function and maximum cost threshold.
Movement to find mates was simulated as an inverse square of the cost distance with a threshold adjusted by multiplying 12000m (as per Thatte et al. 4 ) by the mean resistance per unit area. Movement of dispersing individuals in our model is sex-specific and uses a negative exponential function with two maximum dispersal distances of 200,000 cost units for females and 600,000 cost units for males.
Simulated movement across the landscape in this approach was limited by landscape resistance, defined by individual movement and the cost-distance between occupiable locations. For this purpose, a costdistance matrix was generated from the resistance surface (250m resolution) and all occupiable points (defined above based on a maximum density of 1 individual per 100km 2 ). Individuals in our models are included in our simulations in two ways. First, an initial DPKY population (timestep 1) was generated via a sample of 30 points, representing the upper range of the current tiger population estimates within DPKY boundaries (Ash et al. 12  (TMOR50), and 75% (TMOR75). In R 13 , each translocated individual in CDPOP input files was assigned a randomly generated number between 0 and 1. If this randomly generated number fell below a given threshold (i.e., <0.25 for TMOR25, <0.5 for TMOR50, and <0.75 for TMOR75) the individual was removed.
One time step in our simulation represents one non-overlapping generation within the population and simulations were run for 20 generations (approximately 100 years, as per Thatte et al. 4 ).

Process overview and scheduling
Our simulations were carried out in CDPOP in single, distinct time-steps, representing a non-overlapping generation. The CDPOP simulation works in several discrete steps within each simulated generation: 1. Mating pairs are selected probabilistically according to the mating function described above.
2. These mating events produce offspring based on the demographic function described above, and with Mendelian recombination of genotypes from their parents.
3. The adults die, as this is a non-overlapping generation model. If a given point is occupied at the time of translocation, placement of the translocated individual was random among available grid cells in the landscape. If a population reached carrying capacity during a generation in which a translocation was to occur, individuals would not be introduced into the DPKY population for that generation.
At the end of these steps a new set of adults is distributed in the landscape and these steps are then iterated across 20 generations. At each generation we track the number of individuals in the population (N), their locations, allelic richness (AR), and observed heterozygosity (Ho). The entire process was repeated for each factor combination (APL, SA%, TSEX, TFRQ, and TMOR) across 100 Monte Carlo replicates.

Design concepts
Basic Principles: Demographic patterns governing species population dynamics and genetic structure are dependent on spatially-dynamic processes across heterogeneous landscapes. Spatially-explicit models can enable understanding of the effect of translocation of individuals into a population on population and genetic structure in relation to the specific spatial configuration of a real-world landscape. Interaction: Interactions between individuals are of two kinds. First, mating occurs between individuals of opposite sexes that are selected based on cost proximity and reproductive state (e.g., adults and females that have not already mated in that time step). Second, dispersal is limited to territories not already occupied, limiting population density and leading to density-dependent mortality if no available territories are accessible within the dispersal capability of the individual.
Stochasticity: Individual life-history and behavioural characteristics emerge stochastically in simulations, governed by probabilistic draws from offspring, sex, movement, and mortality parameters.

Initialization
An initial population (timestep 1) was generated via a sample of 30 points, representing the upper range of the current tiger population estimates within DPKY boundaries 12 . This was generated probabilistically and proportional to predicted tiger habitat suitability. Source point locations were determined by subtracting a random raster, with values from 0-1, from a rescaled multi-scale-and shape-optimized tiger

Resistance Surface
The process for developing the resistance surface was originally described by Ash et al. 1  Authors of the study determined that the expert map more broadly captured likely resistance to movement in the study landscape and utilized the expert-based resistance surface in their evaluation of landscape connectivity for tigers. This expert-based resistance layer was used in Ash et al. 8 to represent current estimates of landscape resistance to tiger movement and to act as a baseline for comparison of the effect of various landscape change scenarios.

Mortality Function
In population connectivity studies, mortality may be implicitly incorporated through resistance surfaces 9 , though the relationship between resistance surface and mortality may be unclear. The explicit incorporation of spatially-differential mortality risk in landscape-scale population modelling studies has been shown to dramatically affect conclusions pertaining to landscape connectivity and population persistence 1,8,[10][11][12][13] .
In order to account for heterogeneity in mortality risk across the study landscape, we incorporated a habitat-based mortality function (MH10) which emerged as the most plausible function among a series of functions assessed in Ash et al. 8 . This mortality function was developed as a product of predicted tiger habitat suitability, landscape resistance, and protected area status with values representing the probability of a tiger that has dispersed to a given pixel dying before reproducing. Forested areas of low resistance (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20) were adjusted by a habitat suitability prediction surface from layer from Ash et al. 2 , inverted and rescaled between 0-20 percent mortality probability. In effect, the highest quality habitat would have a low probability of mortality (~2%-8%) with a max mortality probability of 20% (corresponding with a mean ~80% annual survival documented in well-protected forest in Western Thailand; 14 ). Outside PAs, the mortality risk in forest areas was doubled, with most unprotected forest corresponding to roughly 40% mortality probability before further adjustments. Similar to Kaszta et al. 13 , long-distance effects of anthropogenic structures, such as roads, dams/reservoirs, and urban areas were included as buffers, with mortality probability declining logarithmically up to the distance of 17,400m for roads and dams, and over 5,000m for urban areas.

Supplementary Results
Here, we discuss the relative effect of translocations in terms of the number of simulations in which population, allelic richness, and heterozygosity values were (1) stable/increased, (2) declined, and (3) zero (extinction) across various scenarios.

Effect of translocations on genetic trajectory
We evaluate the simulated genetic trajectory of the DPKY population in two ways, based on allelic richness (AR) and observed heterozygosity (Ho). In simulations in which no translocations occurred AR and Ho values declined or were zero due to extinction in 99.1% and 99.7% of simulations, respectively (Table 2).
Comparing adjustment of starting mean alleles per locus for DPKY individuals (APL), declines in AR and Ho were more frequent in simulations where DPKY individuals had a higher initial average number of alleles per locus. AR was lower in 51.9% of simulations in APL2 compared to 82.1% in APL4 and Ho was lower in 75.8% of APL2 simulations compared to 93.2% in APL4 (Table 2). (-29.1%) of simulations in these respective scenarios.

Variance of Population and Genetic Trajectories
Overall, population trajectories even in simulations using the same parameters, were highly variable,