Population Characteristics of Feral Horses Impacted by Anthropogenic Factors and Their Management Implications

Feral horses form relatively stable harems over time that are characterized by long-lasting bonds among their members, a characteristic that makes them an exceptional case of a social system among terrestrial ungulates. Their social system has been described as uniform despite the wide differences in their environment and demography. Horse populations subjected to human interference often show higher levels of population instability that can ultimately compromise their reproductive success. In this article, we describe demographic and dynamic changes of a Portuguese population of Garranos in Serra d’Arga (SA), which is impacted by human and predation pressures, over six breeding seasons. Furthermore, we tested several hypotheses related to the impact of anthropogenic disturbance on the structure and dynamics of this population. Our results revealed that the SA population had relatively little human interference at the start of the project in 2016. This was supported by the natural composition of the herd (total number of individuals, 206), which consisted of several single- and multi-male harems (n = 17 and 7, respectively) and bachelor males (n = 9). However, from 2017 to 2021, SA’s Garrano population suffered a drastic decline. Approximately two-thirds of the individuals and all bachelor males disappeared, and 76% of adult female transfers occurred after the death or disappearance of the harem male. Predatory pressures and poor management of the population, which allowed illegal human interference, contributed to this population crisis. A low population growth rate, reduced birth and foal survival rates, in addition to a delayed primiparous age were observed in this population and exacerbated after its drastic decline; suggesting the viability and survival of this Garrano population were compromised. Investigating the population demographic changes and their causes and consequences can provide guidelines for managing populations and help fight the extinction of horse breeds.


INTRODUCTION
Feral horse populations, similar to their wild counterparts (i.e., zebras and wild asses), are shaped by births, deaths, and membership changes, including dispersals, transfers, and integration of new members in a group (Ransom et al., 2016). However, how these demographic factors influence equid social dynamics and ecology varies among species and populations subjected to different environmental pressures (Linklater, 2000;Ransom et al., 2016). In the case of horses in particular, management practices also shape population structure and dynamics (see for a review: Linklater, 2000). Horses are quite unique regarding their ecology; for example, they occupy a diverse range of habitats within all latitudes excluding the polar circle, from xeric to mesic environments (e.g., Klingel, 1972;Rubenstein, 1981;Hoffmann, 1983;Berger, 1986;Kaseda et al., 1997;Feh, 1999;Linklater, 2000).
Horses, including the Przewalski's horse (Equus ferus przewalski) and two out of three zebra species (plains zebra, Equus burchelli, and mountain zebra, Equus zebra), display a female defense polygynous mating system in which one or more males (stallions) can form a breeding group (or harem) with several unrelated adult females and their immature offspring (Berger, 1986;Rubenstein, 1994;Stanley et al., 2018). Multi-male groups have been mostly reported for populations with low levels of management (Tyler, 1972;Gates, 1979;Kaseda, 1981;Goodloe et al., 2000). Some studies have also reported the presence of all-female groups, and aggregations of young females and males; however, these groups are relatively scarce and of short duration (Keiper, 1976;Berger, 1977;Linklater et al., 2000; but see: Clutton-Brock et al., 1976). Furthermore, they also form all-male groups, known as bachelor groups, which are less stable and consist of young males, or colts, that have dispersed from their natal groups (Berger, 1986;Linklater, 2000;Pinto and Hirata, 2020). The presence of bachelor groups is often a sign of population stability, and several management practices involve controlling the surplus number of males, i.e., bachelors and pre-dispersal young males, which results in highly skewed female sex ratios (Linklater, 2000;Boyd and Keiper, 2005;Górecka-Bruzda et al., 2020). Young females, or fillies, also disperse from their natal groups before the onset of adulthood and can either integrate into a new harem or join bachelor males and form a new harem (Monard et al., 1996;Kaseda et al., 1997;Goodloe et al., 2000). The process of dispersal in horses is hypothesized to have evolved to avoid incest (Monard et al., 1996;Linklater and Cameron, 2009). Moreover, the age at which females disperse can be affected by several factors, including prolonged maternal investment (Kaseda et al., 1997; but see: Monard et al., 1996), onset of puberty (Rutberg and Keiper, 1993;Monard et al., 1996), and reproductive opportunities (i.e., number or existing harems or males; Monard et al., 1996).
Another unique aspect of horses relates to the stability and complexity of their social system. Feral or semi-feral horses form stable social groups (Linklater, 2000;Boyd et al., 2016;Maeda et al., 2021a). Year-round group stability has been reported for different populations under contrasting environmental, demographic, and management pressures (Klingel, 1982;Keiper and Sambraus, 1986;Linklater, 2000;Boyd and Keiper, 2005). As a result of their harem stability, horses develop long-term inter-sexual social bonds that have been proposed to promote cohesiveness to cope with outgroup male harassment and predation pressures (Feh, 1999;Linklater et al., 1999;Cameron et al., 2009); this ultimately contributes to increased female reproductive success ). Long-lasting bonds among females are maintained even after the death of the harem male; a new male often takes over the female unit (Klingel, 1982;Rubenstein, 1994;Mendonça et al., 2021). Despite the ubiquitous bond stability, adult females can voluntarily transfer to other harems (Berger, 1986;Rutberg, 1990;Monard et al., 1996) depending on factors such as stallion age, subordinate male presence, and harem size (Stevens, 1988).
Harem composition and size seem to vary both between and within populations, and is strongly linked to ecological factors such as resource abundance and predation pressures; however, more importantly, in managed populations, harem size is closely related to management practices (Boyd and Keiper, 2005). Populations with more males tend to form smaller groups as a result of a stronger competition between males for a smaller number of females (Berger, 1986;Keiper and Sambraus, 1986). Thus, populations subjected to greater levels of human interference (e.g., male-skewed removals) should comprise larger group sizes. Naturally, we would expect a female-skewed sex ratio in horse populations; similar to many polygynous ungulate species, adult male horses suffer greater mortality compared with females, and these disparities have often been attributed to costs associated with the evolution of male reproductive competition (Berger, 1983;Garrott, 1991). However, in populations subjected to predation pressures, the scenario can be exacerbated as males face more risks because of their predator-directed aggression (Berger, 1983;Garrott, 1991). Management practices have also been reported to limit population changes (e.g., Linklater et al., 2004). Nonetheless, ecological factors, such as resource abundance and predation pressures appear to be the major constraints of population growth (Keiper and Houpt, 1984;Goodloe et al., 2000;Dawson and Hone, 2012). In managed populations, human interference can often be a confounding factor and a hindrance to testing hypotheses related to animal ecology and social system evolution. Nonetheless, studying human impact can allow us to better elucidate the resilience and behavioral plasticity of feral horse populations.
The current study focuses on a Garrano population from Serra d'Arga (SA), northern Portugal. This population is quite unique because they are organized in a multi-level society that consists of multi-male, single-male, and bachelor groups in undefended territories (Maeda et al., 2021a); this differs from their Galician and other northern Portuguese counterparts, which are mostly managed under a traditional husbandry system (Lagos, 2013;Nuñez et al., 2016;Ringhofer et al., 2017). However, this population has been impacted by illicit and uncontrolled human intervention, particularly targeting the males, and it is subjected to predation pressures by the Iberian wolf (C. lupus signatus) (Primavera et al., 2015;Nakamura et al., 2021).
The aims of the current study were twofold: first, to describe a feral Garrano population in SA by examining its demography and membership changes over six breeding seasons using four measures of population stability: harem size, sex ratios, foaling rate, and female transfers; second, to use this population's demographic data to explore the effects of anthropogenic impacts on equid ecology, social dynamics, and ultimately on their reproductive success based on a set of hypotheses. First, we hypothesized that harem size may be affected by adult sex ratio (Boyd and Keiper, 2005). Populations with more males tend to form smaller groups as a result of a stronger competition between males for a smaller number of females (Berger, 1986;Keiper and Sambraus, 1986); therefore, we expected the group size to increase with the decrease of the number of males throughout the years. Second, we hypothesized that ratio differences tend to arise later during adulthood as an outcome of the costs of male-male competition for females, which results in greater male mortality (Berger, 1983). Although this hypothesis can only be partially tested given the original female sex biases of the population, we expected a decrease in the female-to-male ratio over the years and a female:male ratio to be 1:1 at birth. Third, we hypothesized that group transfers and integration, which are linked to higher stress levels in females, can compromise their body condition and ultimately affect the sex ratio at birth Nuñez et al., 2014). According to Trivers and Willard (1973) model horse females with better body condition should produce more sons; therefore, we could expect a variation in foal sex ratio after a period of more instability for females. However, body condition of the mothers only seems to affect the sex of the foal at time of conception . Because all the disturbances in the group compositions happened after the window of conception, we expected the sex ratio to be constant over the study period. Fourth, as harem stability seems to be linked to higher reproductive rates (Berger, 1983;Stevens, 1988), we predicted that this population would have relatively low foaling rates and low survival rates following the years of greater disturbance. Additionally, females that suffered fewer group changes should also be expected to have higher foaling ratios. Fifth, male tenure and female movements are good indicators of population stability (Rubenstein, 1981;Boyd et al., 2016). A short male tenure was expected for the SA population based on the management practices and predation pressures to which this population is subjected. An increased rate of female transfers would be expected in the years of less stability, which should coincide with the higher levels of male disappearance/removal. Sixth, if young females dispersed to avoid incest (Monard et al., 1996), we expected the replacement of the harem male to avert the process of dispersion. Seventh, dispersal age can be delayed by decrease in the number of males and groups (Monard et al., 1996). Thus, we expected a delay in dispersal age as the number of groups and males declined. Moreover, an overall delay in the age of primiparous females was also expected as a consequence of the delayed dispersal.

Population and Study Site
This study focused on the Garrano population living in SA (825 m a.s.l), a mountain range covering an area of 4,493 hectares in northwestern Portugal (41 • 48 W, 8 • 42 N) that is characterized by a Mediterranean climate with Atlantic influences. The highaltitude plateau (700-800 m a.s.l.) used by the horses during the breeding season mainly consisted of wetlands and shrublands dominated by heather species (Erica ciliaris, Erica tetralix, and Calluna vulgaris), gorse species (Ulex minor and Ulex europaeus), and broom species (Genista anglica and Genista micrantha) interspersed with extensive granite outcrops (Gonçalves et al., 2016). Currently, SA belongs to the Natura 2000 Network and Sites of Community Importance (SCI) and a place where local people still survive on traditional practices, such as agriculture and livestock husbandry (Primavera et al., 2015).
Garranos are an ancestral breed of ponies that inhabit the heath and moorlands of the northwestern Iberian Peninsula. They are small, sturdy, and robust; well adapted to mountainous environments; and have an average weight of 290 kg and average height of 1.30 m at the withers (Morais et al., 2005). The Garrano population from SA is not managed under a traditional husbandry system and it is not confined by any artificial or topographical barriers (Lagos, 2013;Nuñez et al., 2016;Ringhofer et al., 2017). A traditional system for Garranos basically consists of strict control over the number and pedigree of adult males, which results in an extremely skewed sex ratio population structure (sometimes with six males for 100 females). In populations where the traditional husbandry system prevails, horses are managed by an association that consists of local owners that organize yearly round-ups where they register and/or remove some foals to be sold. The Garrano owners receive annual subsidies for each horse from the European Union and are compensated for cases of predation by Iberian wolves. These regulations were established to encourage locals to preserve the breed, which was declared endangered by the European Union in 1994 (ACERG, 2021). Currently, there are 1,563 female and 208 male registered Garranos in Portugal (ACERG, 2021).
Horses in SA are not pure Garranos, which is common for most free-ranging horses found in northern Portugal and Galicia (Morais et al., 2005). Pure and crossbreed Garranos started to be released around the mid-70s in the SA mountain. Horses in SA coexist with other species of farm animals, such as endemic breeds of cows and goats. Initially, horses were introduced to the mountain as bait to protect other farm animals that yielded greater economic value from predation by Iberian wolves. Additionally, releasing horses yielded economic advantages to the farmers because they did not have to spend money caring for the horses and would be compensated if there was damage by wolves (Milheiras and Hodge, 2011).
However, over the years and changes in the laws regulating equids registration, the ownership of the horses in SA has become unclear. Very few horses in this population are privately owned by livestock farmers. Despite not being managed under the traditional husbandry system, removal of Garranos in SA takes place arbitrarily, and mostly illicitly, often between breeding seasons. Moreover, this Garrano population is subject to predation pressures by the Iberian wolf (Álvares, 2011). A wolf pack has been reproducing regularly on a yearly basis since 2014, and previous studies estimated a pack size of an average size of 10 ± 1.79 from 2014 to 2018 (Primavera et al., 2015;Freitas, 2019;Nakamura et al., 2021). Strong positive selection by wolves on horses in contrast to other prey species was also reported; horses represented more than 80% of the biomass consumed by wolves in 2018 and 2019 (Freitas, 2019).
Field research on this population has been ongoing since 2016; researchers have mainly focused on the horses' social system and socio-spatial behaviors using novel technologies and noninvasive methodologies (Ringhofer et al., 2017Inoue et al., 2018Inoue et al., , 2020Go et al., 2020;Mendonça et al., 2020Mendonça et al., , 2021Pinto and Hirata, 2020;Maeda et al., 2021a,b). During the study period, horses were not subjected to any anti-parasitic or other medical treatment and showed a "good" body condition (approximately 3 on a scale from 1 to 5; Pinto et al., under review). Moreover, they were not provided with additional food except during August-September 2016; this was after a fire in SA and we have reports of local people bringing hay and water to temporarily feed the horses.

Data Collection
We usually started field observations at 09:30 and continued until 18:00 each day (see for more details on the methodology: Ringhofer et al., 2017;Maeda et al., 2021a;Mendonça et al., 2021). Horses were classified according to their sex and age-class as adult females, or mares, and adult males, or stallions (>4 years old); young females, or fillies, and young males, or colts (>1 year old and <4 years old);and foals (<1 year old). All individuals were given names and identified based on visible phenotypic traits such as body color, presence/absence and shape of white marks on the face and/or feet, and mane color and direction (Supplementary Figure 1).
Groups were named after the male of the group or, in the case of two-male groups, with both of their names. The number and composition of individuals in each group and membership changes were recorded daily during six breeding seasons (May-August) from 2016 to 2021, with the exception of 2020 when observations were only conducted during July. Because of weather constraints and the natural horse group movements, the number of groups that we were able to observe in a day varied. During the breeding seasons from April to July, the horse groups gather on the plateau of the mountain (Maeda et al., 2021a), which facilitates the simultaneous observation of several groups. A team of three to five researchers and volunteers followed the horses daily during the observation period, and conducted population censuses, registering the membership of each group (presence of individuals, dispersal, and integration), including the numbers of foals born. SI and MR followed one group per day (Inoue et al., 2018Ringhofer et al., 2020). RM and PP follow three to six groups per day (Mendonça et al., 2021). TM focused on inter-group dynamics and with the help of a drone could follow up to 20 groups per day (Maeda et al., 2021a). All researchers also checked for the presence and membership of the surrounding groups.
The number of horses reported in this study corresponded to the maximum number observed in each breeding season. Foaling rates were estimated by dividing the total number of females that foaled by the number of females of reproductive age for a given reproductive season. Fecundity rates were estimated by diving the total number of female offspring born by the number of females of reproductive age for a given reproductive season. Foal survival rates were estimated by dividing the number of foals that survived until the next breeding season by the number of foals that were born in the previous breeding season. Individuals were confirmed dead if their carcasses were found or if they were detected in wolf scat as a result of predation or scavenging. Carcasses were found by researchers or horse owners. Telemetry data from a collared wolf, between August and November 2018, provided information on wolves' feeding sites, including kill sites (Freitas, 2019). The horse DNA present in the wolves' scat was analyzed and compared with genotypes obtained through the analyses of feces of identified horses from 2017 to 2021 (n = 200) (Freitas, 2019;Mendonça et al., 2021). Horses found in wolf scat and identifiable carcasses were considered "dead." Identifiable carcasses or DNA from carcasses found in wolf kill sites by researchers were considered "depredated." We discussed horse removals with some horse owners that witnessed the removals or removed the horses themselves. However, there was no reliable and systematic record of horse removal. Horse removal could happen clandestinely, and in these cases we were unable to track it down or to obtain information regarding the removals. In the absence of a carcass, or if we could not confirm the horse was removed, we assumed that the individual disappearance cause was "unknown" if it had not been seen within a 1-year period, which corresponded to the interval between two breeding seasons. Moreover, horses that disappeared during the breeding seasons were considered "dead" because removals did not take place during those periods; therefore, the most likely cause of disappearance was a natural death.
We considered an adult female or a young female to have changed groups if they were found with a different harem from one breeding season to another. To assess group stability, we used the definition in our previous study (Mendonça et al., 2021); for a group to be considered stable, the proportion of individuals who transferred/disappeared or integrated into a new group had to be <0.5 (i.e., the number of individuals who left or integrated into the group had to be smaller than the number of individuals who composed the core group).
Male tenure was described as the years a male could hold a harem within the study period (2016)(2017)(2018)(2019)(2020)(2021). In the case of males that formed two or more groups, tenure was calculated by averaging the male's tenure in each group. Males in multi-male groups were analyzed individually to calculate tenure.
To determine the annual finite rate (λ) of population increase we used the Leslie age-structured model (Leslie, 1945). To build the model, we used age-specific fecundity rates, survival and horses' presence for each female age class at time t (corresponding to each breeding season). We built a matrix with three cohorts, foals (1 years old), young females (>1 year old and ≤3 years old), and adult females (>3 years old). Each age class at time t + 1 depended on: (1) counts of individuals in each age class at the time t; (2) per capita fecundity rates in each age class at the time t; and (3) probability of surviving of an individual in age class t surviving to the next age class t + 1.

Genetic Relatedness Analyses
During the breeding seasons of 2017-2021, between May and August, 290 fecal samples were collected from visually identified horses. Individual identification was achieved using a set of 11 horse microsatellite markers (AHT4, ASB2, ASB17, ASB23, HMS1, HMS2, HMS3, HMS6, HMS7, HTG6, and VHL20), which matched the recommended core or extra panels for individual genotyping of horses by the International Society for Animal Genetics (ISAG, 2017). The protocol for DNA extraction, markers, and genotyping is described in Mendonça et al. (2021). Genetic relatedness (r) was calculated between 200 individuals using the Triadic likelihood estimator (TrioML) implemented in COANCESTRY 1.0 (Wang, 2011), with 10,000 bootstrap replicates and allele frequencies of 84 individuals from the SA Garrano population. The methodology to select the best relatedness estimator from 11 microsatellite markers genotyped from SA Garranos is described in Mendonça et al. (2021). The empirical TrioML estimate values, ranging from 0 to 1 for each dyad where, for example, the probability that two siblings share the same gene by descent is 0.5 and for cousins is 0.125 (Nowak, 2006).

Demographic Changes Over the Years
We conducted binomial tests to assess the skewness of the sex ratios in adult individuals and foals for each breeding season (2016)(2017). To examine the population stability throughout the years, we examined the effect of breeding seasons on the following stability measures: harem size, sex ratio, foaling rates, and female transfers, we built a generalized linear mixed model (GLMM) using the R package 'glmmTMB' (R × 64 3.5.0 1 ; Brooks et al., 2017). The model that examined the effect of breeding season on harem size was built with a Poisson error structure and logit link function. For each group, we considered the count of adult females over six breeding seasons. The identity (ID) of the group was included as a random effect; we examined 42 groups and 105 data points. The remaining models were built with binomial structure and logit link function. In the model of both adult and foal sex ratio, females were scored as 1 and males as 0. For the adult sex ratio, we examined 178 individuals and 545 data points, and included female ID as a random factor. For the foal sex ratio, we examined 140 individuals born throughout the study period, and included foal ID as a random factor. For the foaling rate model, each adult female that had a foal was scored 1, and those who did not have a foal were scored 0. We examined 125 females and a total of 389 data points. Each data point corresponded to a female's reproductive outcome in a given breeding season. In the model investigating female transfers, each female that transferred was scored 1, and those who did not transfer were scored 0. In this model, 109 adult females and a total of 268 data points were examined. We did not include young females that dispersed for the first time or females that disappeared between the breeding seasons.

Best Fit Models for Hypotheses Testing
We tested the effect of group disruption, which was inferred based on the proportion of males disappearing each year, on the stability measures mentioned in the previous section. Because different breeding seasons were characterized by different demographic factors (i.e., number of males, females, and groups), we evaluated the impact of these factors to test the proposed hypothesis. Therefore, GLMMs with logit link function were built to investigate the best explanatory models for harem size (Poisson structure, number of females of reproductive age per group); sex ratios (binomial structure, 0 -male; 1 -female), foaling rates (binomial structure, 0 -did not foal; 1 -foaled), and female transfer (binomial structure, 0 -did not transfer; 1 -transferred). Because the models did not converge when the breeding season and the proportion of males' that disappeared were simultaneously added as predictors, we did not test for breeding season in these models.
For the adult sex ratio model, the total number of males and females each year and the proportion of males that disappeared before the breeding season were tested as fixed effects; we examined 161 individuals and a total of 399 data points. For the foal sex ratio model, we tested the effect of the proportion of males that disappeared, total number of females, and total number of foals; we examined 94 foals. For the model examining female transfers, we tested the effect of the proportion of males that disappeared, the number of groups in the population before transfer, and harem size before transfer. We examined 109 females and a total of 268 data points. For the foaling rate model, we tested the harem size, number of males per group, and proportion of males that disappeared between breeding seasons. We examined 110 females and a total of 284 data points. After a stepwise selection procedure, we reported the models with the best fit, i.e., the lowest AIC value.
Finally, we tested the effect of female group changes on foaling rate (binomial structure, 0 -did not foal; 1 -foaled). We tested female transfer prior to giving birth as a categorical predictor (yes or no). We analyzed 16 adult females that remained present during the six breeding seasons and 80 data points.

Changes in the Population Structure, Size, and Survival Rates
During the 5 years of our study, which included six breeding seasons, the population decreased from a maximum of 159 individuals (>1 year old) to 49 individuals ( Table 1). The average annual finite rate of population increase (λ) was 0.82 and the exponential rate of increase (r) was −0.23. From August 2018 to the breeding season of 2019, the population suffered a severe decline; approximately two-thirds of the population (including all bachelor males) and more than half of the harems disappeared (Figure 1 and Table 1). The causes of the disappearance known are reported in Figure 1. From 2020 to 2021, this trend seemed to be reversed; the population size showed a 25% increase ( Table 1).
Between the breeding seasons of 2016 and 2021, the number of observed harems in SA varied from 7 to 30 ( Table 1). The number of harems increased from 24 to 30 in the 2017 BS, which resulted from individuals finding new groups (n = 4) and one group splitting into two smaller harems (n = 2). In 2017, the proportion of multi-male groups increased to 37% (n = 11). In the 2018 breeding season, the number of groups decreased to 25 harems (24% were two-male groups) as a result of the disbandment of some groups; however, the number of bachelor males increased (n = 13). From 2016 to 2018, bachelors banded in two separate all-male groups or remained solitary. Consequently, membership changes in group composition were observed ( Table 2). Despite the decline in the number of harems over the years, the number of adult females per group (harem size) remained constant (GLMM ,  Table 3). In the best fit model, the only predictor that increased the fit of the model was the proportion of males that disappeared, which had no significant effect on harem size (GLMM, Z = −0.91, P = 0.36, Table 4).
The adult Garrano population (>3 years old) was characterized by a stable female-biased sex ratio, with an average of 2.45 ± 0.28 (mean ± SD) females per male ( Table 5). Adult sex ratio remained constant over the years (GLMM ,  Table 3) and neither the proportion of males that nor the total number of males in the population affected the adult sex ratio (GLMM , Table 4). Foal sex ratio showed a very mild sex bias toward females, 1.04 ± 0.53, and no effect was found for breeding seasons or males' disappearance (GLMMs, Tables 3, 4).
Our results on foaling rates showed that there was a significant decrease of the foaling rates in the 2019 breeding season (GLMM, Z = −3.42, P < 0.001, Table 3 and Supplementary Table 1), which coincided with the population decline that occurred after the 2018 breeding season. In the best fit model, the proportion of males that disappeared negatively affected the foaling rate (GLMM, Z = −3.63, P < 0.001, Table 4). Foal survival within the first year of life averaged 0.18 ± 0.24 and fluctuated between 0.00 and 0.23 during the 2016-2021 period (Figure 2).

Male Tenure
Over the study period (2016-2021) and out of the 40 groups observed, only two groups (Takaoka&Uozu and Gozen&Nagaoka) remained stable throughout the whole period (Supplementary Table 3). The average male tenure for all males that formed harems was 1.9 ± 1.2 years during the six breeding season period. Out of 44 males, only 11% (n = 5) had a tenure of 5 years; four belonged to two-male groups (Takaoka&Uozu and Gozen&Nagaoka) and one belonged to a single-male group (Toki) that was part of a three-male group in 2016. In six multimale groups, one of the males disappeared or became solitary, changing the group status from multi-male to single-male. Five bachelor males formed their own harem during our study period; however, one group, Seitsu's, was artificially formed because the male was kept in captivity with the females, which were also captured from SA, for 2 years.
During our study period, three of the males lost their harem at least once and formed a new harem in the subsequent year. Of all 45 harem males, only 20% remained alive in the 2021 breeding season. The cause of their disappearance could be predation, natural causes leading to scavenging, or removal by humans. Of all 36 harem males that disappeared, we did not identify the 1 Females and males between 1 and 4 years old. 2 Weighted average. Group size includes males and females older than 1 year. Females/harem includes only adult females (older than 3 years old) and corresponds to the harem size.x refers to the yearly mean.
FIGURE 1 | Population changes by age class for females and males (>1 year old) and foals (<1 year old), representing the individuals that survived and disappeared in SA from 2016 to 2021. Individuals were categorized as "Dead" if they were identified in the wolves' scats or disappeared during the breeding season when removals did not take place; "Depredated" if the carcasses were found in a kill site (Freitas, 2019); "Unknown" if the horse had disappeared outside the breeding season period and its death had not been confirmed; and "Removed" if we could confirm they were taken from the mountain by owners or locals.
cause of disappearance for 50% of the males (n = 19); 17% were considered dead (n = 6) and 31% (n = 11) were removed by humans; one of the human-removed males was released back into the population at the end of 2020.

Female Transfers
During the study period, 43% (n = 47) of SA population females underwent a transfer event, but only 14% (n = 15) of the females voluntarily transferred. From 2016 to 2021, we recorded 63 inter-breeding season transfer events, 76% (n = 48) of which occurred after the disappearance or death of the male where females either sought to integrate into another harem or joined a solitary male (i.e., passive transfers). Females that voluntarily transferred in the presence of the harem males were always transferred alone. On passive transfers, females could either transfer alone or in groups (collective transfers, Table 6). Out of 48 females that passively transferred, 49% (n = 31) transferred in groups were collective transfers. From 2016 to 2021, 14 clusters of females that included 2.6 ± 0.61 females either transferred to another harem together or were taken by a new male.  Table 3). The interaction between the proportion of males that disappeared and harem size showed a significant effect on the female transfer rates (GLMM: Z = 2.92, P = 0.004, Table 4). During the years of more disturbance, females from larger groups were more likely to transfer; this was likely due to the group destruction affecting groups with larger harem size.
Out of the 16 adult females that were present during the six BS, three had four foals and two of them were never observed with a foal. The average foaling rate per female was 0.33 ± 0.22 (Supplementary Table 2). We did not find any relationship between the number of changes suffered by females and the foaling rate (GLMM: Z = −1.10, P = 0.27).

Young Female Dispersal and Primiparous Age
In the SA population, 12/22 young females (55%, >1 year old) were observed dispersing from their natal groups between 2016 and 2021 ( Table 7). The average age of dispersal was approximately 1.83 ± 0.55 years (range: 1-3 years, n = 12). Out of the 12 young females, 58% (n = 7) were closely related to the harem male (r > 0.50), and 25% (n = 2) of females were not closely related to the male (r = 0 and r = 0.16). The other two young females were from a two-male group and not related at least to one of the males (r = 0 and r = 0.13).
Three 5-year-old females that were followed since they were born did not disperse. One (Imizu) was born in a two-male group and was not related to one of the males ( Table 7); Mizuho was born in a three-male group in which two of the males left the group (one of them was not related to her, the other we could not determine) and the remaining one was not genetically related to her. Sue was born in a two-male group and related to one of the males. After both males disappeared, a non-related male took over.
The average primiparous age during our study period in SA was 4.4 ± 0.5 (range: 4-5, n = 5). Of all the females that were born between 2016 and 2021 and gave birth, 40% (n = 2) dispersed from their natal group, whereas the remaining 60% (n = 3) either stayed in their natal group, which was originally a multi-male group (Imizu, Mizuho) or suffered a change of male (Sue), and none of them birthed offspring with their potential father.
Frontiers in Ecology and Evolution | www.frontiersin.org The name of the group is given by the name of the male(s). Groups in bold are multi-male groups. Group size include all individuals older than 1 year. (+), females integrating into the harem. Hyuga became a single-male group, and the females remained with him. A young male dispersed and joined the bachelor males. 10 Kujin disappeared. The female remained with Kitakami, and six more females joined them. 11 Natori disappeared and Suwa joined the bachelor males. 12 Hyuga lost his females and was observed with the bachelor males during the 2018 breeding season. 13 Hirosaki followed Takaoka&Uozu during the 2018 breeding season. 14 Enzan disappeared. Chuo&Enzan became a single-male group, and the females remained with Chuo. 15 Tennoji disappeared together with the females from his group. Nanba formed a new single-male group. 16 Two foals became young individuals and were therefore considered in the group size. 17 Souja was observed following Unnan group and ended up replacing Unnan. Unnan was found injured and alone in June 2021. Significance was set at P < 0.05 and highlighted in bold. We reported the models with the lowest AIC for the predictors tested. Significance was set at P < 0.05 and highlighted in bold.
Frontiers in Ecology and Evolution | www.frontiersin.org FIGURE 2 | Foaling rate and foal survival over the study period. Foaling rate was calculated as the number of foals born divided by the number of females of reproductive age. Survival rate was calculated as the number of foals that survived the first year within a year divided by the number of foals born in each breeding season. Percentage of transferred females were calculated by diving the number of females transferring between breeding seasons, by the total number of females. The number of females disappearing between breeding seasons were excluded from the total. Highlighted in bold is the period when more than half of the harem males disappeared. 1 Transfer with the harem male present. 2 Transfer after disappearance of the harem male. 3 Annual average.

Population Structure and Demography
Feral horse populations have been considered uniform with regard to their social and spatial structure; nevertheless, group composition and size tends to vary across different populations depending on the management level and ecological pressures to which they are subjected ( Table 8). In unmanaged populations of the United States and Argentina, 83-86% of the groups included a single male, 11-17% of the groups included two males, and 0-3% included only females with their offspring (Berger, 1986;Scorolli, 1999;Goodloe et al., 2000). By contrast, other feral populations showed higher percentages (67%) of multi-male groups (Linnartz and Linnartz-Nieuwdorp, 2017). In populations where non-skewed sex removals occurred, up to 58% of groups included multiple males (with a maximum of five males) Linnartz and Linnartz-Nieuwdorp, 2017).
In SA, despite the skewed sex ratio, the percentage of singlemale groups varied between 66 and 82% during a period of six breeding seasons, a feature similar to what is observed in unmanaged populations. A maximum of four males was observed in a group in the 2016 breeding season. However, groups with more than two males were not stable; in the following breeding season, they split up into smaller harems. It is possible that those groups were formed by young males that were on the verge of dispersing because sometimes it was difficult to discern between immature (around 3 years old) and adult horses. Nonetheless, most of the two-male groups were relatively stable during the observation period, with males in two-male groups showing the longest tenure in our study-5 years.
Harem size has been reported to fluctuate between 2.2 and 4 individuals in unmanaged feral horse populations (Klingel, 1982;Keiper and Sambraus, 1986;Boyd and Keiper, 2005;Scorolli, 2007;Roelle et al., 2010) and 1.8-5.7 in female-skewed populations (Kaseda, 1981;Keiper and Sambraus, 1986). Some managed horse populations can have a maximum of 11-22 adult females (Pacheco and Herrera, 1997;Linklater et al., 2000). The average harem size of the SA population across the six breeding season period ranged from 4.5 to 6.7 females and the maximum number of females per group was 12, which was similar to that of managed and sex-skewed populations. Despite selective pressure suffered by males and disbandment of some harems, harem size remained relatively constant over the years. This indicates that, in response to sudden changes in their group composition, horses are prone to reorganizing their social composition. After the loss Updated and modified version of Table 2 from Ringhofer et al. (2017), and Table 4.1 of Mills and McDonnell (2005); we used the published report that mentioned both group sizes and the observed number of mixed-sex groups from Mills and McDonnell (2005).
of the harem male, we observed three different scenarios: (1) another male, often a bachelor, tried to integrate the females' group, (2) a female group moved together roaming around and tried to join several groups until they settled for a group, and (3) a group of females splitting and each subgroup integrated into a new harem (Mendonça, personal observations). Because females tend to join other harems or be joined by a male, we did not observe stable all-female groups as have been reported by some other studies (Nelson, 1980;Kaseda, 1981;Goodloe et al., 2000). Considering the high predation pressure during our study period, integration in a new harem is likely to be important for female survival and therefore reproductive fitness (Freitas, 2019;Mendonça et al., 2020). Factors such as familiarity between females or their ecological needs may have played roles in making this decision because familiarity was shown to be a crucial factor that predicted affiliation in population studies (Monard et al., 1996;van Dierendonck et al., 2004;Mendonça et al., 2021). It has been hypothesized that group size may be affected by adult sex ratio (Boyd and Keiper, 2005). Populations with more males tend to form smaller groups as a result of a stronger competition between males for a smaller number of females (Berger, 1986;Keiper and Sambraus, 1986). Our results showed that the decrease in the number of males did not influence the harem size, which remained relatively constant over the years. However, our results corroborate a previous study on Misaki horses in Japan (Kaseda and Khalil, 1996), which reported an average harem size of 3.8 in a 3:1 female-to-male ratio. In SA, with a similar female-to-male ratio (2.45:1), the average number of females per harem was 3.7.
Adult sex ratio did not decrease over the years, as predicted. Therefore, our results do not support the hypothesis stating that sex ratio differences tend to arise later during adulthood as an outcome of the costs of male-male competition for females (Berger, 1983). Nonetheless, as expected, the female-to-male sex ratio at birth was close to 1:1 and did not show significant changes over time, neither after the decline of two-thirds of the population triggering substantial changes in the groups. However, because these changes occurred after the breeding season and far from the conception window of females, these results were to be expected . Nonetheless, a foal sex ratio of 1:1 suggests that approximately two-thirds of the males in the SA population have disappeared or died. Arbitrary male-biased removals and increased predation by wolves might have been the main causes for the skewed sex ratio. In 2017 and 2018, horses in SA represented 85-93% of wolves' consumed biomass, and horse consumption seemed to be biased toward adult horses (Freitas, 2019). However, in light of recent events and given the very good body condition of the horses in this population (Pinto et al., 2022), we think that predation on adult male horses might have been overestimated in a previous report and may be a result of scavenging rather than predation (Freitas, 2019). Previous studies reported that feral horse populations usually have relatively high adult survival rates with low annual variation (Duncan, 1992;Cameron et al., 2001;Dawson and Hone, 2012;Zabek et al., 2016) and foals are more likely to be predated by wolves because of their vulnerability in the first months of age compared with healthy adults (Vos, 2000;Lagos, 2013). Recently, an adult male was confirmed to be illegally shot before being consumed by wolves, and that was not the first report of male horses being shot in SA (Mendonça personal observations). Therefore, we can attribute the higher rates of adults disappearing to removals rather than predation, and the high consumption of adult horses by wolves to scavenging rather than predation.
Foal survival in this population, average, 0.18 ± 0.24 (range, 0.00-0.72), was severely lower compared with the mean foal survival rates for feral horse populations, average 0.84 ± 0.002 (range, 0.09-0.97) (see for a review: Ransom et al., 2016). Foal survival plummeted from 2016, and predation by the Iberian wolf might have strongly contributed to this (Freitas, 2019). In populations where feral horses coexist with predators, a high percentage of foal mortality is often attributed to predation in contrast with disease, accidents, or other natural causes. Previous studies on predation by the Iberian wolves on Garranos reported that wolves have been responsible for 89-96% of horse mortality (Gomes and Oom, 2000;Lagos, 2013). During 2018 and 2019, a strong positive selection by wolves was observed in SA in contrast with other prey species, including both livestock and wild ungulates (Freitas, 2019). The impact of wolf predation in our study population was estimated to be 6 horses/wolf/year, which represented a loss of approximately 41% of individuals each year (Freitas, 2019). From 2017 to 2019, all the foals disappeared during our observation period, when removal by humans did not occur. This suggests that human interference was not likely to have had a major impact on foal mortality or disappearance during those years. However, foal carcasses were hard to detect because they can be entirely consumed within hours to a few days (Lagos, 2013;Mendonça et al., 2020); therefore, our reports might underestimate foal depredation. Nonetheless, we found some evidence of wolf attacks, such as foal carcasses and foals with open wounds on their rump, and evidence of wolf consumption on foals through analysis of their scats (Freitas, 2019;Mendonça et al., 2020).
The horse removals SA population underwent allied with high predation pressures greatly limited the population growth throughout the years. The average finite growth of SA population (λ = 0.82) reflected an overall declining trend, much lower than the ones observed for other feral horses' populations (λ = 1.03-1.49), even compared to populations constraint by ranging space and compromised nutrition (Keiper and Houpt, 1984;Goodloe et al., 2000;Cameron et al., 2001;Dawson and Hone, 2012;Zabek et al., 2016). The highest values have been reported for populations that experienced food abundance and low predation pressures (Eberhardt et al., 1982;Garrott et al., 1991;Duncan, 1992;Grange et al., 2009).
From 2020, foal survival considerably increased, which coincided with an increased foaling rate, a decrease in horse removals, and fewer passive female transfers. Moreover, an apparent decrease of Iberian wolf presence in SA was also reported (Álvares, personal communication). However, the foal survival reached 0% in an adjacent Garrano population managed under a traditional husbandry system during the 2021 breeding season. Further studies on populations that coexist with predators are needed to address how stability and management practices affect or are affected by variations in predation pressures.

Harem Stability and Membership Changes
Harems are often relatively stable over time. However, the SA population was characterized by substantial instability during the six breeding season period, particularly between 2017 and 2019. Over the study period, 43% of the females transferred at least once; however, only 14% of these females voluntarily transferred. The total percentage of adult females that voluntarily transferred between groups varied largely across populations, from 2 to 30% (Berger, 1986;Stevens, 1988;Rutberg and Keiper, 1993;Monard et al., 1996;Goodloe et al., 2000). In the SA population, the rates of voluntary transfers (0-13%) fall within the interval reported by several studies with different levels of management. Nonetheless, when considering both voluntary and passive transfers, our results far exceeded those observed for other managed horse populations; this indicates that the major cause of the excessive female transfer could be caused by the disappearance of males.
Male tenure and the proportion of groups that remained stable are good indicators of population stability (Boyd et al., 2016).
Studies on feral horses have reported relatively high population stability, with less than a third of the original groups being disbanded and approximately 30-37% of the males reaching maximum tenure during 3-8-year study periods (Berger, 1986;Linklater et al., 1999;Goodloe et al., 2000;Cameron et al., 2001;Scorolli and Cazorla, 2010). Those results highly contrast with ours because only 11% of the males in our population reached maximum tenure and only two groups remained stable for the whole study period.
Harem instability seems to be linked to lower reproductive rates (Berger, 1983). In the SA population, the average foaling rate per female was approximately 0.33 foals/year. According to previous studies, it is common for females to foal every other year, or even two out of every 3 years (Tyler, 1972;Keiper and Houpt, 1984;Cameron and Linklater, 2000); thus, SA foaling rates were lower than what was observed for other populations. Overall, the number of foals born per adult female in this population varied from 2016 (44%) to 2019 (9%). The values obtained for the SA population were much lower than the ones reported for either stable (71%) or unstable (44%) harems (Goodloe et al., 2000). Removal of putatively pregnant females after the breeding season might have contributed to the extremely low foaling rates in 2019. However, it could be possible that abortions may have happened in the SA females that were triggered by the replacement of a male and interaction with a new male or high stress levels associated with transfers (Nuñez et al., 2014). Previous studies reported the occurrence of abortions in different populations after replacement of the harem male by an outgroup male (Berger, 1983;Kirkpatrick and Turner, 1986). Moreover, preliminary results indicate that female fecal glucocorticoid metabolite levels were higher in the 2017 and 2018 breeding seasons (Mendonça, unpublished data), which could have affected fetal loss and, therefore, reproductive success in the following year (Keiper and Houpt, 1984;Romero et al., 1998). However, a recent study reported harem male takeovers or changes did not have any influence in the pregnancy rate in female Konik horses . Despite the small data set, our results also show that, among the 16 females that remained throughout the whole study period, there was no tendency for a correlation between changes in the group for females and foaling rate.
Young female dispersal age in SA was similar to what was reported in previous studies on feral horses, from 1 to 3.5 years (Monard et al., 1996;Kaseda et al., 1997;Goodloe et al., 2000). Contrary to our predictions, we did not observe a delay in the age of dispersal; however, the change in the harem male might have prevented the dispersal of young females. Four females that did not disperse from their natal group had suffered a change in the harem male or belonged to a multi-male group where at least one of the males was not genetically related to them. Of all the females that have dispersed, 58% were closely related to the male, one of the females was born after the mother migrated to a new group, and we could not confirm the paternity for two young females.
Nonetheless, the primiparous age for SA females was slightly delayed, at 4-5 years, compared with other studies; 3 years is often reported, and around 12-36% of 2-year-old females were observed to have offspring (Berger, 1986;Kirkpatrick and Turner, 1986;Zabek, 2015). However, other studies reported that females around 2 years old are reproductively inactive (Keiper and Houpt, 1984;Cameron and Linklater, 2000;Goodloe et al., 2000). This delay in the reproductive age of young females may also be a consequence of the disturbances this population underwent; either because females did not survive past puberty, or the survivors suffered changes in their group composition after dispersal, or both. However, given our small sample size and because most of the young females disappeared, we have to be cautious when drawing conclusions. To better understand the effects that these disturbances had on the age of first reproduction of females, research should continue beyond periods of great instability.

CONCLUSION AND MANAGEMENT IMPLICATIONS
Our results show that, except for the female-biased sex ratio, the SA population had relatively low levels of human interference from 2016 to 2018, despite the female-skewed sex ratio. This is supported by the existence of several multi-male groups and bachelor males organized in a multi-level society (Maeda et al., 2021a). Starting from 2017 and peaking at the end of the 2018 breeding season, the SA population suffered a drastic decline that was mainly due to bad management policies aligned with predation pressures. This sudden population decrease had severe impacts on the foaling rate and could compromise the viability and survival of this population. Predation seems to be the main cause of foal mortality, but human interference also plays a critical role in the decline of the adult population. Nonetheless, in the last 2 years (2020-2021), the population's decreasing trend appeared to be curbed as foaling rates and foal survival increased and the removals ceased. To better understand and project the population growth and therefore its recovery, more data post-disturbance (from 2019) are needed. We are aware that observing the horses only during breeding seasons drastically reduces the ability to follow the fate of the individuals. Consequently, efforts are being made to improve this; for example, researchers are conducting year-round monitoring and censusing of the population.
To improve the SA population and populations facing a similar scenario, we first have to investigate the real motivations behind male removals, which might have been for economic or selfish reasons, and offer viable alternatives. Second, the reinforcement of existing laws that protect horse welfare and their owners is paramount. Illegal removals and shootings should be further investigated and horse owners duly compensated. Finally, a coalition among horse owners, researchers, and local government representatives should be created to develop and smoothly execute consensual collective decisions on management of the population by prioritizing the horses' welfare and the livelihoods of the coalition's representatives.
Research on equine behavior and ecology can provide essential guidelines to aid population management and fight the extinction of Garranos and other endangered horse breeds. We suggest minimizing human interference with these feral populations. This should start by not arbitrarily and irresponsibly removing the adult males because they are the key for group protection, both against predation or harassment from other males and for genetic variability. Working closely with researchers is of utmost importance because it will allow horse owners to better understand horses' behaviors and needs, which is crucial for implementing adequate management policies adjusted to the needs of the respective populations. If we enable these circumstances to be met, we will be on the right track toward regenerating horse populations and mitigating the impacts of predation pressures.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because the current study used only non-invasive methodologies.

AUTHOR CONTRIBUTIONS
RM, PP, and MR collected the data. RM analyzed the data and wrote the first draft of the manuscript. All authors contributed to the final version of the manuscript.

FUNDING
This study was financially supported by JSPS Core-to-Core CCSN, JSPS-LGPU04, JSPS KAKENHI (Nos. 16H06283 to Tetsuro Matsuzawa, 18H05524 to SH, 17H05862 and 19H00629 to SY, and 18K18342 to MR) and Kyoto University SPIRITS to SY.