High Degree of Multiple Paternity and Reproductive Skew in the Highly Fecund Live-Bearing Fish Poecilia gillii (Family Poeciliidae)

Multiple paternity is a common phenomenon within the live-bearing fish family Poeciliidae. There is a great variety in brood sizes of at least two orders-of-magnitude across the family. However, little is known about the ramifications of this remarkable variation for the incidence and degree of multiple paternity and reproductive skew. Mollies (subgenus Mollienesia, genus Poecilia) produce some of the largest broods in the family Poeciliidae, making them an excellent model to study the effects of intra-specific variation in brood size on patterns of multiple paternity. We collected samples of the live-bearing fish Poecilia gillii from 9 locations in Costa Rica. We measured body size of 159 adult females, of which 72 were pregnant. These samples had a mean brood size of 47.2 ± 3.0 embryos, ranging from 4 to 130 embryos. We genotyped 196 field-collected specimens with 5 microsatellite markers to obtain location-specific allele frequencies. In addition, we randomly selected 31 pregnant females, genotyped all their embryos (N = 1346) and calculated two different parameters of multiple paternity: i.e., the minimum number of sires per litter using an exclusion-based method (GERUD) and the estimated number of sires per litter using a maximum likelihood approach (COLONY). Based on these two approaches, multiple paternity was detected in 22 and 27 (out of the 31) females, respectively, with the minimum number of sires ranging from 1 to 4 (mean ± SE: 2.1 ± 0.16 sires per female) and the estimated number of fathers ranging from 1 to 9 (mean ± SE: 4.2 ± 0.35 sires per female). The number of fathers per brood was positively correlated with brood size, but not with female size. Next, we calculated the reproductive skew per brood using the estimated number of sires, and found that in 21 out of the 27 multiply sired broods sires did not contribute equally to the offspring. Skew was not correlated with either female size, brood size or the number of sires per brood. Finally, we discuss several biological mechanisms that may influence multiple paternity and reproductive skew in P. gillii as well as in the Poeciliidae in general.


INTRODUCTION
Female multiple mating, or polyandry, often leads to offspring being sired by multiple fathers, a phenomenon which is referred to as multiple paternity. Polyandry is prevalent in the animal kingdom, but it remains debated why it evolved in in most animal taxa (Jennions and Petrie, 2000;Zeh and Zeh, 2001;Simmons, 2005). Nonetheless, there are several hypotheses concerning how polyandry and multiple paternity can increase female fitness (Keller and Reeve, 1995;Zeh, 1996, 2001;Jennions and Petrie, 2000;Stockley, 2003;Simmons, 2005).
First, females might select multiple mates in order to receive direct benefits, such as nuptial gifts (Thornhill, 1976;Kondoh, 2001), or fertilization assurance (Sheldon, 1994). Second, polyandry can also arise because of indirect genetic benefits (Yasui, 1998;Jennions and Petrie, 2000), in three ways: (i) Females can mate with multiple males to improve the genetic quality of the offspring by 'trading-up' (Pitcher et al., 2003). This hypothesis predicts that a female that mates with a 'genetically inferior' male can increase her fitness by mating with another male of higher genetic quality (Jennions and Petrie, 2000;Simmons, 2001;Pitcher et al., 2003), (ii) Multiple paternity can also increase a female's fitness by increasing the genetic diversity of offspring (Yasui, 1998;Fox and Rauter, 2003). If the environment is unpredictable, the genetic diversity of the offspring could increase a female's reproductive success by ensuring that some of her offspring will survive when the environment changes, which is a form of bet-hedging (Loman et al., 1988;Yasui, 1998;Fox and Rauter, 2003;Garcia-Gonzalez et al., 2015). However, bet-hedging strategies only seem to outweigh the fitness costs of multiple mating in certain situations (Yasui, 1998;Fox and Rauter, 2003;Yasui and Garcia-Gonzalez, 2016), (iii) Multiple mating can further arise because it can promote post-copulatory sexual selection (e.g., by enhancing sperm competition; Parker, 1970). The enhanced sperm competition results in an increased fertilization chance by genetically fit males (Zeh and Zeh, 1997), thereby increasing the fitness of a females' offspring. Third, multiple paternity can also arise without female choice, as males force mating upon the females, i.e., coercive mating (Clutton-Brock and Parker, 1995). By doing so, males can increase their reproductive output. When more males mate coercively with the female, this can lead to polyandry and multiple paternity. In this way, multiple paternity can arise without increasing the fitness of a female.
In the live-bearing fish family Poeciliidae, both female mate choice (often in association with male courtship behavior) and coercive mating occur (Magurran, 2011;Rios-Cardenas and Morris, 2011). Moreover, multiple paternity is widespread in this family; multiple sired broods are found, for instance, in the genera Gambusia (Greene and Brown, 1991;Zane et al., 1999), Heterandria (Soucy and Travis, 2003), Poecilia (Travis et al., 1990;Kelly et al., 1999;Hain and Neff, 2007;Neff et al., 2008;Girndt et al., 2012), and Xiphophorus (Luo et al., 2005;Simmons et al., 2008;Tatarenkov et al., 2008;Paczolt et al., 2015). In multiple sired broods, it is common that sires do not contribute equally to the offspring, a phenomenon that is referred to as reproductive skew. Previous studies have shown that in natural populations of Poeciliidae often pronounced male reproductive skew is observed, this is for instance the case for swordtails (Luo et al., 2005;Simmons et al., 2008;Tatarenkov et al., 2008), guppies (Hain and Neff, 2007;Neff et al., 2008), and mollies (Girndt et al., 2012). There are several factors that can cause and enhance male reproductive skew in natural populations. For instance, when females mate with multiple males this will intensify sperm competition, which often results in one or a few male(s) being more prominent in siring offspring (Parker, 1970;Constantz, 1984;Evans, 2010). An increase in the number of sires may hereby lead to an increased chance of skew. Next to sperm competition, cryptic female choice can also cause paternity skew. This refers to a post-copulatory sexual selection process in which the female biases sperm use, thereby altering the paternal contribution (Eberhard, 1996;Firman et al., 2017).
In the Poeciliidae, there is a remarkable variation in brood sizes of at least two orders-of-magnitude, both within and among species. Most earlier multiple paternity studies have focused on species with small broods, such as Poecilia reticulata (Kelly et al., 1999;Hain and Neff, 2007;Neff et al., 2008), Heterandria formosa, (Kelly et al., 1999;Soucy and Travis, 2003), Xiphophorus multilineatus (Luo et al., 2005), X. nigrensis (Smith, 2014), and X. birchmanni (Paczolt et al., 2015). Only a few earlier studies focus on species with large broods, i.e., Poecilia latipinna and Gambusia affinis, but most of these only analyze a small part of the actual brood (Travis et al., 1990;Greene and Brown, 1991;Girndt et al., 2012), except for a study by Gao et al. (2019) that studied multiple paternity in whole broods of Gambusia affinis and G. holbrooki. Therefore, little is currently known about how the large variation in brood size affects multiple paternity, the number of sires per brood and male reproductive skew within species.
Prior studies on poeciliid fishes have found a positive relationship between multiple paternity and female traits such as brood size and female size (Neff et al., 2008;Girndt et al., 2012). Several hypotheses have been proposed to explain why multiple paternity may be correlated with brood size and female size. The number of offspring a female can carry is restricted by both the limited space available in the body cavity (physical limitation; Reznick and Miles, 1989;Pires et al., 2011) and the exacerbated negative consequences of an increase in reproductive allocation on a female's swimming ability (physiological limitation; Fleuren et al., 2019;Quicazan-Rubio et al., 2019). These factors explain (at least partly) why female size is often positively correlated with brood size (Cheong et al., 1984;Reznick et al., 1992Reznick et al., , 1993Neff et al., 2008;Schrader and Travis, 2009;Hagmayer et al., 2018Hagmayer et al., , 2020: larger females can physically carry larger broods than smaller females. Simple mathematics dictate that the larger a brood, the larger the absolute number of potential sires that can contribute to that brood (Avise and Liu, 2011): i.e., a brood of two offspring can be sired by a maximum of two fathers, while a brood of 30 offspring can be sired by a maximum of 30 different males. Several studies have furthermore shown that males prefer to mate with larger females (Bisazza et al., 1989;Ptacek and Travis, 1997;Dosen and Montgomerie, 2004;Herdman et al., 2004;Hoysak and Godin, 2007), presumably because they increase their chance of producing more offspring as these females generally carry more eggs (Schlupp, 2018). These above processes may explain the observed positive associations between female size, brood size and multiple paternity.
In this study, we investigate the degree of multiple paternity and reproductive skew in the live-bearing fish Poecilia gillii (Poeciliidae) across its South-Western range in Costa Rica. Members of the subgenus Mollienesia (genus Poecilia) are an excellent model to study the effects of brood size variation on patterns of multiple paternity at an inter-specific level, because they produce some of the largest broods in the family (Reznick and Miles, 1989;Winemiller, 1993;Johnson and Bagley, 2011). Specifically, we will test to what extent female size and brood size within P. gillii correlate with the (i) degree of multiple paternity (i.e., the number of sires within a brood) and (ii) male reproductive skew.

Animal Ethics
The study was assessed by the institutional Animal Welfare Body (AWB) of Wageningen University (the Netherlands). The AWB judged that the study complied with the Dutch Act on Animal Experiments (AAE), which complies to European Directive 2010/63/EU. All fish were collected and exported following the regulations of the Costa Rican government under permit number SINAC-CUS-PI-R-005-2017, SINAC-PNI-ACLAP-004-2020, and SINAC-ACOSA-DT-PI-INV-003-2020.

Field Collection
We collected individuals of the live-bearing fish species Poecilia gillii (family Poeciliidae) in Costa Rica using seine and cast nets during the dry season from 27 February to 12 March 2017. To minimize the potentially confounding effect of location on multiple paternity, we randomly sampled 9 different locations in the Rio Terraba General and the Rio Coto drainages in the Province of Puntarenas, Costa Rica (Table 1). To obtain sufficient samples to estimate population allele frequencies per location, additional samples were collected from two of these locations (specifically; N = 9 individuals from Pedegroso, and N = 8 from Tigre) from 18 February to 16 March 2020, which resulted in a total of 196 specimens for molecular analysis (range of specimens per location: 15-26). The collected specimens were euthanized in the field with an overdose of MS-222 (Tricaine methanesulfonate) and preserved in 96% ethanol. The fish were subsequently transported to the Pollux lab at Wageningen University (the Netherlands) and stored in a freezer at −20 • C until further processing.

Female Life History Traits
We measured the total length (TL; i.e., length from snout to the distal end of tail), standard length (SL; i.e., total length excluding length of caudal fin) and total wet mass of 159 females, of which 72 were pregnant (Supplementary Table S1). For all pregnant females, we measured the wet mass of the reproductive tissue (ovary with embryos) and counted and staged their developing embryos, following Haynes (1995) staging classification for  Table S1). Based on these measurements, we calculated the reproductive allotment (RA) by dividing wet mass of the ovary by the total wet mass of the female (Supplementary Table S1). For molecular paternity analysis, we randomly selected 32 Poecilia gillii females that were pregnant with broods that were at developmental stage 20 or later (early eyed;Haynes, 1995) to ensure that sufficient DNA could be extracted from the embryos (resulting in a total of 1371 embryos for further analysis).

DNA Extraction
DNA was extracted from the 32 gravid females of Poecilia gillii and all their embryos (N = 1371) for multiple paternity analysis, and an additional samples (males, females and juveniles) to estimate population allele frequencies. For the DNA extraction, we used the following adapted protocol from the Wizard Genomic DNA purification kit (Promega). Entire embryos and tailfin clips from adults were taken and dried for 10 min at room temperature (20-25 • C) and stored at -20 • C before DNA extraction. Lysis was performed using 300 µL Nuclei Lysis Solution/EDTA mix (by preparing a total of 310 µL mix: 60 µL of 0.5 M EDTA and 250 µL of Nuclei Lysis Solution) and 9 µL of Proteinase K (20 mg ml −1 ) and by incubating at 55 • C for 2-2.5 h (vortexed every 30 min). Protein precipitation was performed by adding 100 µL Protein Precipitation solution, after which the samples were vortexed for 20 s, cooled on ice for 5 min, and finally centrifuged for 4 min at 16000 × g. DNA was precipitated from the supernatant by adding 300 µL isopropanol, mixing, and leaving the samples at room temperature (20-25 • C) for 30-60 min, then centrifuging at 16000 × g for 1 min to form a pellet. Next, the supernatant was decanted and the pellet (DNA) was washed with 300 µL 70% cooled ethanol (0 • C), mixed, centrifuged at 16000 × g for 2 min and dried at room temperature for 20-30 min. Finally, the pellet was dissolved in 25 µl Tris (10 mM) for embryonic samples and in 50 µl Tris (10 mM) for maternal samples. The stock DNA was diluted to 20 ng/µL for further analyses.

Microsatellite Analysis
All obtained DNA samples (N = 1567) were analyzed using 5 polymorphic microsatellite loci (  We show the number of alleles (N A ), the allelic range, the observed (H O ) and expected heterozygosity (H E ), and the exclusion probability (E 1 ) per locus and for the 5 loci combined. *This location was excluded from the multiple paternity analyses, due to the low exclusion probability.
We checked the PCR products by gel electrophoresis and diluted samples based on band thickness. We analyzed the diluted PCR samples with capillary gel electrophoresis using an ABI3730 Sequencer (Applied Biosystems). For the analysis, we used a total volume of 10 µL per sample containing 1 µL of diluted PCR samples (all markers combined) and 9 µL ladder mix (0.5% (v/v) Liz500 size ladder in formamide). We genotyped the samples using GeneMapper software version 4.0. We scored peaks with automatic binning in GeneMapper but checked all samples manually and rescored when automatic scoring errors occurred. Allele count, range and observed and expected heterozygosity were calculated per sampling site ( Table 2), using GenAlex version 6.5 Smouse, 2006, 2012). Exclusion probabilities were furthermore calculated per sampling site, using GERUD 2.0 (Dodds et al., 1996;Jones, 2005). As exclusion probabilities are determined by both the number and frequency of alleles in a population (Dodds et al., 1996), we included location-specific allele frequencies in our analyses that were based on the 196 field-collected individuals. The exclusion probability represents the probability of a marker to exclude a random unrelated male from paternity for a motheroffspring pair (Dodds et al., 1996;Wang, 2007) and provides a comparative measure of marker information. However, it has been argued that the exclusion probability should not be used as a measure of confidence in parentage analyses, because: (i) exclusion probability calculations do not account for marker errors, and (ii) the probabilities are derived for single-offspring parentage tests and therefore do not correct for experimentwide error that arises from the hundreds of comparisons in a typical parentage study . Nevertheless, the low exclusion probability in the Incendio (<0,5; Table 2), prompted the exclusion of the single female from this location from our subsequent paternity analyses.

Paternity Analysis
Multiple paternity was quantified for the remaining 31 pregnant females and their 1346 embryos (using two software programs: GERUD version 2.0 (Jones, 2005) and COLONY version 2.0 (Jones and Wang, 2010). GERUD is an exclusion-based method that calculates the minimum number of sires for multi-locus data using co-dominant markers, such as microsatellites (Jones, 2005). It excludes maternal alleles in the offspring to calculate the number of alleles that offspring received from their fathers (Jones, 2005). Based on this data, it reconstructs the paternal genotypes to calculate the minimum number of sires (Jones, 2005). As missing data is not accepted in GERUD, we excluded samples that lacked data for one or more markers (see Table 3). COLONY calculates the estimated number of sires based on maximum likelihood methods, using simulated annealing to search for a global optimum (Jones and Wang, 2010). The program uses multi-locus data and considers all mothers and embryos jointly for paternity assignments (Jones and Wang, 2010). When marker diversity is limited (e.g., low expected heterozygosities and/or few microsatellite markers with a limited number of alleles), COLONY tends to overestimate the true number of sires (Sefc and Koblmüller, 2009). We performed a full-likelihood analysis in COLONY per sampling location (i.e., population), using the maximum run length and the highest precision. Population allele frequencies were estimated with GenAlEx version 6.5 Smouse, 2006, 2012), utilizing the extra samples for each location ( Table 2). Allele frequencies for each locus and sampling location were then applied to run the COLONY analyses. The expected genotyping error rate was set at 0.025, as suggested by Wang (2004). Mating systems for males and females were set on polygamous, and, as we do not have any information on sibship size in Poecilia gillii, no prior was set for sibship size. COLONY can run with missing data, but with limited genetic info paternity estimates become less accurate. Therefore, we excluded samples from the analysis that had missing genotypes for more than 2 out of the 5 loci (see Table 3). Since location-specific allele frequencies are used in these analyses to determine the number of compatible fathers, we calculated the number of alleles (N A ), observed heterozygosity (H O ), and expected heterozygosity (H E ) using GenAlex version 6.5 Smouse, 2006, 2012), and the exclusion probabilities (E 1 ) using GERUD 2.0 (Dodds et al., 1996;Jones, 2005) per locus per sampling location ( Table 2).

Reproductive Skew
We calculated the male reproductive skew per mother using the estimated number of sires obtained from COLONY. The effective number (N eff ) of sires was estimated by N eff = 1/ (O i /BS) 2 , in which BS is the total brood size of a female, and O i is the number of offspring assigned to each sire i. We then calculated the reproductive skew (R skew ) by dividing the effective number of sires by the actual number of sires, and subtracting the outcome from 1 (Neff et al., 2008).

Data Analysis
All the analyses were carried out in R v 3.5 (R Core Team, 2019). Mixed models were fitted using the lme4 package (Bates et al., 2015), and significance tests for the fixed effects were performed with lmerTest (Kuznetsova et al., 2017). Poisson models were checked for over-and under-dispersion using functions implemented in the in the DHARMa package (Hartig, 2018). In the case of data dispersion, we switched to a Conway-Maxwell-Poisson (CMP) distribution, which is a viable count distribution that generalizes the Poisson distribution in light of data dispersion (Shmueli et al., 2005). CMP models can be fitted using the glmmTMB package (Brooks et al., 2017). Model equations of the analyses discussed below can be found in the Supporting methods.
To evaluate the effect of maternal standard length on the probability of being pregnant and brood size, we fitted two separate generalized linear mixed effect models that each included maternal standard length as an explanatory variable and sampling location as a random effect (a categorical variable with 9 levels). In the first model (model 1), the probability of being pregnant (a dichotomous variable: 1 = pregnant, 0 = not pregnant; N = 159 females) was fitted using a binomial response distribution and a logit link function. In the second model (model 2), brood size (a discrete variable based on the embryo counts in a brood) was analyzed for pregnant females only (N = 72) using a Poisson distribution and a log link function.
The potential effects of brood size and female size on multiple paternity were assessed by fitting the (i) minimum number of sires inferred from GERUD (model 3), and (ii) estimated number of sires inferred from COLONY (model 4) in generalized linear mixed effect models with a log link for the Conway-Maxwell-Poisson-distributed response to handle the under-dispersed data in (i), and the Poisson-distributed response in (ii). These models included female standard length, brood size and their interaction term as fixed effects and sampling location as a random effect (the latter to account for between-location variation in multiple paternity that is not accounted for by the fixed effects). The correspondence of the number of sires obtained with the two different programs (GERUD and COLONY) was evaluated by calculating the Spearman rank correlation coefficient between the minimum (GERUD) and estimated number of sires (COLONY).
To study the paternity skew in broods, we first tested whether the observed contribution of sires in a brood significantly deviates from the expected contribution if there was no skew (i.e., equal  For each female the standard length, total wet mass, reproductive allotment (RA) and brood size is given. For each brood the developmental stage (according to Haynes, 1995) is given. The minimum number of sires was estimated with GERUD, the estimated number of sires and reproductive skew were calculated in COLONY. The embryos analysed in these programs can be less than the actual brood size, as some samples were omitted because of lacking genotype data for one or more markers (for details see "Materials and Methods").
contribution) by means of χ 2 -tests for goodness of fit (Yue and Chang, 2010;Green et al., 2017). To identify potential sources of variation in paternity skew, we then fitted a series of generalized linear models with a logit link function for the quasibinomial-distributed response, as paternity skew values range between 0 and 1. To prevent model over-fitting, we did not include a random sampling location effect. The full model included brood size, female standard length, and the estimated number of sires inferred from COLONY as fixed effects. In addition, we tested all the two-way interactions by adding them to the model case-by-case. Because none of the interaction terms were significant (brood size × female standard length:

Variation in Female Life History Traits
A total of 159 adult Poecilia gillii females collected from 9 locations were studied to determine their life history traits (Supplementary Table S1

Multiple Paternity
Molecular characteristics of the 5 microsatellite loci used to detect multiple paternity in Poecilia gillii are summarized in Table 2. Loci were on average 93.33% (SE = 4.71%) polymorphic. Marker diversity was high with the total number of alleles per marker ranging from 14 (for GA-IV42) to 23 (Pvm16) ( Table 2,  Supplementary Table S2).
The minimum number of sires in the broods of the 31 females was determined with GERUD (Table 3). Multiple paternity was observed in 22 of the 31 studied females, with the minimum number of sires ranging from 1 to 4 (mean ± SE: 2.1 ± 0.16 sires per female). There was a significant positive correlation between brood size and the minimum number of sires (Figure 2B; model 3: Z 25 = 2.393, p = 0.017), but not between female SL and the minimum number of sires (Figure 2A; model 3: Z 25 = 1.937, p = 0.053). In addition, there was no significant interaction effect between SL and brood size (model 3: Z 25 = 0.375, p = 0.708). Next, the estimated number of sires contributing to the litters of the of the 31 P. gillii females was calculated with COLONY (Table 3). Multiple paternity was observed in 27 of the 31 females, with the estimated number of sires ranging from 1 to 9 (mean ± SE: 4.2 ± 0.35 sires per female). The estimated number of sires was positively associated with brood size (Figure 2D; model 4: Z 26 = 2.198, p = 0.028), but not with female standard length (Figure 2C; model 4: Z 26 = 0.936, p = 0.349). There was no significant interaction between brood size and female standard length on the estimated number of sires (Z 26 = -0.951, p = 0.342). We found a strong positive correlation between the minimum number of sires calculated in GERUD and the estimated number of sires calculated in COLONY (Figure 3; Spearman's rank correlation: r = 0.698 p < 0.001).

Paternity Skew
Broods with multiple paternity (determined with COLONY) were associated with considerable paternity skew, ranging from 0.08 to 0.53 (Table 3) with a mean (± SE) of 0.27 ± 0.03. A value of zero implies that all sires contribute equally to the brood (no skew), and a value of one implies sires differ maximally in their contribution to the brood (maximal skew). In most broods (21 out of 27 multiply sired broods), the observed paternal contribution among offspring (i.e., skew) significantly deviated from the expected contribution ( Figure 4A; goodness-of-fit χ 2tests p < 0.05). There was no significant correlation between skew and the number of sires inferred from COLONY (Figure 4B

Multiple Paternity
Local environmental conditions (e.g., light intensity, predation risk, food availability) might play an important role in determining (the incidence of) multiple paternity within species (Kelly et al., 1999;Soucy and Travis, 2003). Moreover, local 'logistic constraints' (sensu Avise and Liu, 2011; e.g., population density, sex ratio, mating strategy) are also likely to influence geographic patterns of multiple paternity. To minimize the potentially confounding effect of local conditions on multiple paternity, we randomly sampled 31 pregnant females from 8 different locations to characterize multiple paternity in P. gillii across its South-Western range in Costa Rica. Moreover, many studies on multiple paternity in species with large broods have analyzed only a fraction, or a set number of embryos, per brood (Travis et al., 1990;Greene and Brown, 1991;Girndt et al., 2012 but see Gao et al., 2019), potentially underestimating the absolute number of sires that may have contributed to the broods. To avoid this, we genotyped and analyzed all embryos (N = 1346) to quantify multiple paternity in our 31 P. gillii females.
We report a high incidence of multiple paternity in P. gillii; ranging from 22 to 27 out of 31 broods (i.e., 71%-87%, based on estimates of the minimum and estimated number of fathers per litter, respectively). Such high incidences of multiple paternity are common within this genus, e.g., in Poecilia reticulata multiple paternity was observed in 95% of the pregnant females (Neff et al., 2008). However, in other genera, lower rates of multiple paternity are also observed, sometimes as low as 23% in Poeciliopsis monacha (Lesie and Vrijenhoek, 1977) or 46% in Heterandria formosa (Soucy and Travis, 2003). In some genera, there is a FIGURE 4 | A Paternity skew in 31 female Poecilia gillii. The bar graph shows the relative contribution in percentages of sires to the broods of 31 females (different shades represent different sires within each brood). Asterisks indicate that observed paternal contribution among offspring differed significantly from the expected (equal) contribution (goodness-of-fit χ 2 -tests, p < 0.05), meaning that paternity was significantly skewed for those broods. The paternity skew and number of sires that contribute to each brood are shown above the graph. Panels considerable variation in the incidence of multiple paternity, ranging from 28% in Xiphophorus multilineatus (Luo et al., 2005) to 84% in X. birchmanni (Paczolt et al., 2015), possibly due to different mating strategies among the species. For instance, in the genus Xiphophorus, males of sword-carrying species show more dominant behavior and females seem to have much lower rates of multiple paternity than species that do not carry swords (Luo et al., 2005;Paczolt et al., 2015). Environmental factors, like predation, can also affect mating strategies and influence the rate of multiple paternity. For example, in guppies, the number of sneak mating attempts and the degree of multiple paternity were significantly higher for high predation populations (Matthews et al., 1997;Kelly et al., 1999).

Effect of Female Traits on Multiple Paternity
We found a significant positive correlation between brood size and the number of sires in Poecilia gillii (GERUD; model 3: Z 25 = 2.393, p = 0.017, COLONY; model 4: Z 26 = 2.198, p = 0.028). This correlation can arise simply because the larger a brood, the larger the number of potential sires that can contribute to that brood (Avise and Liu, 2011), but it could also arise because females with larger broods mate with more males (e.g., possibly to ensure fertilization of all eggs, or to increase the genetic variability of her offspring in bigger broods). To date, correlations between brood size and multiple paternity have been studied in species from four genera of the family Poeciliidae: Gambusia, Heterandria, Poecilia, Xiphophorus. Similar positive relationships have been reported in the closely related Poecilia latipinna and P. reticulata (based on estimated number of sires: Travis et al., 1990;Greene and Brown, 1991;Neff et al., 2008;Girndt et al., 2012;Zeng et al., 2017), and the more distantly related Xiphophorus birchmanni (based on minimum number of sires: Paczolt et al., 2015) and Gambusia holbrooki (based on minimum and estimated number of sires: Zeng et al., 2017). Some studies did, however, not find a positive relationship between brood size and multiple paternity; i.e., no significant correlation was found in Gambusia affinis, Heterandria formosa and Xiphophorus helleri (Greene and Brown, 1991;Soucy and Travis, 2003;Tatarenkov et al., 2008).
We did not find a significant relationship between the number of sires and female size in P. gillii. Most previous studies in the Poeciliidae do not report a significant correlation between female size and multiple paternity (Trexler, 1997;Zane et al., 1999;Neff et al., 2008;Girndt et al., 2012;Paczolt et al., 2015;Zeng et al., 2017), although a significant relationship was found in Gambusia affinis and Poecilia latipinna (Travis et al., 1990;Greene and Brown, 1991). We expected to find a positive relationship between female size and multiple paternity, as many earlier studies have shown a male preference for larger females (including the closely related P. reticulata and P. latipinna; Bisazza et al., 1989;Ptacek and Travis, 1997;Dosen and Montgomerie, 2004;Herdman et al., 2004;Hoysak and Godin, 2007;Schlupp, 2018). But this preference, usually assessed using choice-experiments in the lab, might not reflect actual mating opportunities in natural situations. In mating systems that include sexual harassment, larger females may be better than small females at avoiding unwanted sexual copulation attempts; i.e., earlier studies have shown that larger females are more likely to swim away from approaching males by moving to deeper and faster flowing water to avoid sexual harassment by males (Brewster and Houde, 2003;Croft et al., 2006;Magellan and Magurran, 2006). Laboratory and field observations have shown that Poecilia gillii lacks courtship behavior and relies solely on coercive mating, with males relentlessly chasing females and females continuously fending of obtrusive and unwanted suitors (Ptacek, 2005;Goldberg et al., 2019;Furness, Hagmayer and Pollux, unpublished data). Female poeciliid fish seem to avoid these coercive mating attempts in two ways: by chasing away males or by moving to a habitat with a lower density of males (Croft et al., 2006;Magellan and Magurran, 2006;Magurran, 2011). In both strategies, larger females are likely to be more successful in avoiding unwanted mating attempts than smaller females. In other cases, when some males are dominant and deny female access to other males (e.g., in swordtails: Luo et al., 2005), larger females might actually have lower rates of multiple paternity. This may be the case in Poecilia gillii, as some males are dominant and guard a harem of females/a territory, continuously chasing away competitors (Goldberg et al., 2019;Furness, Hagmayer and Pollux, unpublished data). Further, male preferences could also be more complicated, as was for instance shown in Brachyrhaphis rhabdophora: large males of this species had a preference for large females, but small males had a preference for small females (Basolo, 2004). Because we do not have information on the body sizes of sires in this study, we cannot assess whether this was the case for Poecilia gillii as well.

Pre-copulatory Mechanisms
There are several pre-copulatory mechanisms that might cause male reproductive skew. First, female mating behavior can lead to reproductive skew. A compelling example can be found in the guppy (P. reticulata), where females can increase the duration of copulation when mating with an attractive male, thereby increasing the amount of sperm inseminated (Pilastro et al., 2004(Pilastro et al., , 2007. Second, females may 'trade-up, ' meaning that they mate with additional males if they encounter males of better genetic quality than those previously mated with. This strategy tends to increase reproductive skew, particularly if the female preferentially uses the sperm from the most recent copulation event to fertilize her eggs, a phenomenon that is referred to as last-male sperm precedence (Pitcher et al., 2003). Third, paternity skew can be caused when males in a population vary in their reproductive success. This might occur when there are different male mating strategies in the population. In Xiphophorus nigrensis, for example, different mating strategies exist within a single population: sneaky and displaying males. For this species, reproductive success was significantly lower for sneaky males than for the dominant displaying males (Zimmerer, 1989). In guppies, males also have these two mating strategies, and Pilastro and Bisazza (1999) showed that courting males delivered three times higher numbers of sperm to the females than the sneaky ones, which significantly increased their probability of insemination. This could also be the case in Poecilia gillii, as males show different mating strategies with large and colorful dominant males guarding a harem of females and continuously chasing away competitors, whereas small and dull males mimic females in appearance and sneak mate when the opportunity presents itself (Goldberg et al., 2019;Furness, Hagmayer and Pollux, unpublished data).

Post-copulatory Mechanisms
The high degree of skew could also be an indication that there is strong post-copulatory selection. For instance, sperm competition can lead to increased fertilization chances for fastswimming sperm, leading to a higher contribution of males with fast sperm compared to males with slower swimming sperm (Constantz, 1984;Evans and Pilastro, 2011). Sperm competition is enhanced by polyandry (Constantz, 1984;Evans and Pilastro, 2011), therefore, the high degree of skew found in P. gillii could be caused by the high occurrence of multiple mating in P. gillii. Another factor that influences post-copulatory selection is sperm storage, which is common among poeciliids (Lopez-Sepulcre et al., 2013). When sperm is stored there is a prolonged time for sperm competition (Birkhead and Parker, 1997), and more opportunity for the females to influence paternity by cryptic female choice (Birkhead, 1998;Evans and Pilastro, 2011).

CONCLUSION AND FUTURE PERSPECTIVES
In this study, we assessed the influence of brood size and female length on the degree of multiple paternity and paternity skew in the live-bearing fish Poecilia gillii. Our results show that there is a high degree of multiple paternity and male reproductive skew in P. gillli. Since both the degree of multiple paternity and male mating strategies tend to vary greatly between populations/species of Poeciliidae, future studies should focus on elucidating the link between observed male mating strategies in natural populations (e.g., sneak mating, dominance and courtship), and the degree of multiple paternity and paternal skew. In addition, since many factors can potentially affect multiple paternity in the field, controlled laboratory studies should be performed to help identify the proximate causes of multiple paternity. Finally, the relative contribution of pre-(i.e., behavior) and post-copulatory sexual selection processes (i.e., sperm competition or cryptic female choice) to male reproductive skew can be assessed in controlled laboratory studies.

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

ETHICS STATEMENT
The study was assessed by the institutional Animal Welfare Body (AWB) of Wageningen University (the Netherlands). The AWB judged that the study complied with the Dutch Act on Animal Experiments (AAE), which complies to European Directive 2010/63/EU. All fish were collected and exported following the regulations of the Costa Rican government under permit number SINAC-CUS-PI-R-005-2017, SINAC-PNI-ACLAP-004-2020, and SINAC-ACOSA-DT-PI-INV-003-2020.

AUTHOR CONTRIBUTIONS
MD and BP conceived the project idea. MD, AH, AF, and BP planned the fieldwork and collected the samples. MD and KL carried out the dissections and performed the molecular analyses. MD analyzed the data and wrote the first draft of the manuscript supervised by BP. All authors critically reviewed the initial manuscript, provided helpful input and approved the final manuscript.

ACKNOWLEDGMENTS
We are grateful to Pam Heijmans and Nick Koning for their assistance in the lab. We are also grateful for the possibility to collaborate with the Organization for Tropical Studies and the Las Cruces Biological Field Research Station (San Vito, Puntarenas, Costa Rica). We thank Zak Zahawi (Las Cruces Research Station and Wilson Botanical Garden), Minor Porras, Francisco Campos Rivera and Rebecca Cole (Organization for Tropical Studies, OTS) for their help with the research and export permits, and Ana Rosa Ramírez (Museo de Zoología, Escuela de Biología, Universidad de Costa Rica) for the use of the research facilities.