Effective population size and the genetic consequences of commercial whaling on the humpback whales (Megaptera novaeangliae) from Southwestern Atlantic Ocean

Abstract Genotypes of 10 microsatellite loci of 420 humpback whales from the Southwestern Atlantic Ocean population were used to estimate for the first time its contemporary effective (N e) and census (N c) population sizes and to test the genetic effect of commercial whaling. The results are in agreement with our previous studies that found high genetic diversity for this breeding population. Using an approximate Bayesian computation approach, the scenario of constant N e was significantly supported over scenarios with moderate to strong size changes during the commercial whaling period. The previous generation N c (N e multiplied by 3.6), which should corresponds to the years between around 1980 and 1990, was estimated between ~2,600 and 6,800 whales (point estimate ~4,000), and is broadly compatible with the recent abundance surveys extrapolated to the past using a growth rate of 7.4% per annum. The long-term N c in the constant scenario (point estimate ~15,000) was broadly compatible (considering the confidence interval) with pre-whaling catch records estimates (point estimate ~25,000). Overall, our results shown that the Southwestern Atlantic Ocean humpback whale population is genetically very diverse and resisted well to the strong population reduction during commercial whaling.


Introduction
Estimation of historic (pre-whaling) and contemporary population sizes are important to offer guidelines for managing and restoring populations that suffered overexploitation, such as those of baleen whales (Baker and Clapham, 2004;Jackson et al., 2008). Although abundance is obviously important in the very short term, effective population size (N e ) is a key parameter in the long term, and essential for population genetics, evolutionary biology and conservation biology (Schwartz et al., 1998;Charlesworth, 2009). It is directly related to evolutionary processes, such as rates of genetic drift and loss of genetic variability, levels of inbreeding, and effectiveness of selection (Frankham et al., 2002). Generally, N e is lower than the census population size (N c ) since individuals do not contribute genes equally to the next generation. One important advance has been the development of methods to estimate N e from ge-netic data (Schwartz et al., 1998;Leberg, 2005;Waples, 2005;Palstra and Ruzzante, 2008;Luikart et al., 2010). This approach has provided important information to investigate whale population dynamics, although N e estimates are only available for some species and in a few areas (e.g. Rooney et al., 1999Rooney et al., , 2001Waldick et al., 2002;Roman and Palumbi, 2003;Alter et al., 2007Alter et al., , 2012Ruegg et al., 2010Ruegg et al., , 2013. Another relevant information for management is knowing the impact induced by the extreme reduction of the abundance during 20 th century whaling activities on the genetic diversity of the species, since the loss of genetic variation can impact population viability leading to its premature extinction. However, although conservation actions should prioritize the populations with minor genetic diversity, it is not known which populations went extinct following such bottlenecks.
Humpback whales, Megaptera novaeangliae (Borowski 1781), were among the most exploited baleen whale species by commercial whaling. The populations are found throughout the world's ocean basins, undertaking annual migrations between the low latitude waters, where they breed and calve during the winter-spring months, and the high latitude waters, where they feed during the summer (Dawbin, 1966). The International Whaling Commission (IWC) recognizes seven humpback whale breeding stocks (termed A-G) in the Southern Hemisphere (International Whaling Commission, 2015). In the Southwestern Atlantic Ocean, the humpback whale population wintering along the Brazilian coast (~5°to 23°S) (Andriolo et al., 2010) is recognized as the Breeding Stock A (BSA). The main mating and calving area for this population is in the Abrolhos Bank (16°40'-19°30' S and 37°25'-39°45' W), where most whales concentrate (about 85% of the density) during the breeding season (Andriolo et al., 2010;Bortolotto et al., 2016a). However, in recent years, the number of sightings and strandings has increased beyond the BSA range, indicating recovery of this population and likely expansion of its distribution range (Pretto et al., 2009;Wedekin et al., 2014;Bortolotto et al., 2016b). This population migrates to summer feeding grounds around South Georgia and South Sandwich islands in the Southern Ocean (Engel and Martin, 2009;Zerbini et al., 2006Zerbini et al., , 2011a. Commercial whaling during the 20 th century reduced the worldwide humpback whale population to a small fraction of its pre-exploitation abundance (Tønnessen and Johnsen, 1982). In the Southern Hemisphere, approximately 200,000 humpback whales were hunted from 1904 to 1972, after accounting for the illegal Soviet whaling, mainly by whaling operations around Antarctica feeding areas (Findlay, 2001;Clapham et al., 2009;Allison, 2010). In Brazil, pre-modern whaling began in the early 1600s, ending in the 1830s in southern of the country, with the collapse of the southern right whale population, but lasted until the 1920s in northeastern region due to the high density of humpback whales, mainly in Caravelas, Bahia. It was estimated that between 11,000 and 32,000 humpback whales were captured from 1830 to 1924 (Morais et al., 2017). Modern whaling operations that began in the 20 th century expanded the activities of the whaling stations mainly for the coasts of Costinha (7°S) (between 1910 and 1967) and Cabo Frio (23°S) (between 1960 and 1963), where 352 humpbacks were caught in 1913, but only around 13 whales in 1967, already indicating a significant population size reduction Grangeiro, 1965, 1970;Williamson, 1975). In addition, modern whaling activities in high-density areas in the Antarctic and Sub-Antarctic feeding grounds increased the annual catch to several thousand whales (Findlay, 2001). Only in the surroundings of the South Georgia Island about 22,717 humpback whales were killed between 1904 and 1915, when the exploitation of this stock was most extensive (Edmundson and Hart, 2014). Although the species have been protected since 1966, the former Soviet Union fleet took humpback whales illegally off the central coast of Brazil until 1973(Yablokov et al., 1998. It was estimated that 48,477 humpbacks were caught by Soviet whaling in the Southern Hemisphere be-tween 1948 and 1973, of which 1,407 were caught in the South Atlantic Ocean between 1960 and 1967 (Berzin, 2008).
The BSA population size before the exploitation by modern whaling was estimated using catch records to nearly 24,700 individuals, and it reached its lowest numbers in the late 1950s, when there were less than 500 individuals (Zerbini et al., 2011b). Presently this population is recovering fast (growth rate of 7.4% per annum, Ward et al., 2011) and the abundance in 2015 was estimated around 12,123 individuals (Pavanato et al., 2017). Interestingly, despite these well-documented census size changes in the BSA, recent studies have not detected a genetic bottleneck, that is, a significant reduction in the effective population size of this population (Engel et al., 2008;Cypriano-Souza et al., 2010). The absence of a signal for a genetic bottleneck was explained by the low intensity (in terms of duration and minimum population size) of this bottleneck (Engel et al., 2008;Cypriano-Souza et al., 2010). However, these studies used standard methods (heterozygosity excess, mode-shift and M-ratio tests), that have reduced power to detect moderate bottlenecks (Peery et al., 2012). Besides, none of the genetic studies so far has provided estimates of N e for the Brazilian humpback whale population.
The present study aims to estimate the effective and census population sizes of the Southwestern Atlantic Ocean humpback population, and investigate the effects of commercial whaling in the 20 th century on its genetic diversity, based on the analysis of genotypes constructed from 10 microsatellite loci for 420 individuals sampled off the Brazilian coast.

Sample collection and DNA extraction
Between 1999 and 2007, 379 tissue samples of humpback whales were collected by the biopsy dart procedure (Lambertsen, 1987) at two geographic locations off the Brazilian coast, the Abrolhos Bank (n = 332), in southern Bahia and northern Espírito Santo, and Praia do Forte (n = 47) in northern Bahia ( Figure 1). Only adult animals were sampled within the social groups, which showed a body size longer than 11 meters. Additional samples (n = 41) resulted from individuals adults stranded along the coasts of both states. Samples were preserved in 70% ethanol and were stored at -20°C until processed. Genomic DNA was extracted using proteinase K digestion followed by phenol/chloroform extraction method (Palsbøll et al., 1995) or using the DNeasy Blood and Tissue kit (QIAGEN).
tides: GATA028, GATA053, GATA417, (Palsbøll et al., 1997)]. Genotypes of 268 of the individuals used here were described previously in Cypriano-Souza et al. (2010), and genotyping of the additional samples was conducted exactly as described in that study, including in the same machine.

Genetic variation
The program MICRO-CHECKER v.2.2.3 (Van Oosterhout et al., 2004) was used to identify possible null alleles, large allele dropout, and scoring errors due to stutter peaks. While in the previous study (Cypriano-Souza et al., 2010) locus 417/418 showed a weak sign of null alleles, this was not detected here, since homozygous excess was insufficient to suggest the presence of null alleles.
Genetic diversity was measured as the number of alleles per locus (K), observed and expected heterozygosities (H O and H E , respectively) under Hardy-Weinberg equilibrium (HWE) (Nei, 1978), using FSTAT v.2.9.3 (Goudet, 2002). FSTAT was also used to calculate the measure of F IS (Weir and Cockerham, 1984). Deviations from HWE for each locus (Guo and Thompson, 1992) and linkage disequilibrium between loci were tested using ARLEQUIN 3.5 (Excoffier and Lischer, 2010), corrected for simultaneous comparisons with the sequential Bonferroni test (Rice, 1989).

N e estimation
Two methods were used to estimate effective population size (N e ), both assuming a closed population with discrete generations and random variance in reproductive success. The program NeEstimator v.2.01 (Do et al., 2014) was used to estimate the parental generation N e from genotypic data based on the linkage disequilibrium (LD) method, which calculates separate estimates using different criteria for excluding rare alleles. We used the random mating model and the following critical values (P crit ): 0.05, 0.02 and 0.01, as suggested in the program manual.
An approximate Bayesian computation (ABC) approach implemented in the program DIYABC v.1.0.46 (Cornuet et al., 2008(Cornuet et al., , 2010 was also used to test four different scenarios based in the possible demographic history of this humpback whale population during commercial whaling in the 20 th century. Scenario 1 is a constant size population (no bottleneck), scenario 2 consisted of a population that is still experiencing a bottleneck, scenario 3 is a population that expanded recently from a bottleneck, and scenario 4 is a population that experienced a transitory bottleneck ( Figure 2a). The priors for all parameters were uniformly distributed between specified minimum and maximum values (Table 1 and Figure 2a), which were based on the available information of the whaling history of this population and its present day census data (see Introduction). The demographic parameters were: Scenario 1: N e (longterm historical N e ); Scenario 2: N e 2 (current N e ), N a 2 (pre-Humpback whale genetic diversity 255 bottleneck N e ); Scenario 3: N e 3 (current N e ), N a 3 (N e during bottleneck); t (time since demographic change in scenarios 2 and 3); Scenario 4: N e 4 (current N e ), N b (N e during bottleneck), N a 4 (pre-bottleneck N e ), t2 and t1 (time since the beginning and end of the bottleneck, respectively). All times are in number of generations [generation time of 18 years taking into account the range between 12 and 24 years estimated for the humpback whales (Chittleborough, 1965;Roman and Palumbi, 2003)] from the present, with t2 > t1.
The 10 microsatellites loci were assumed to evolve under the generalized stepwise mutation model (GSM) (Estoup et al., 2002) with the widely used mutation rate (m) range for mammals, from 10 -4 to 10 -3 per generation (Ellegren, 1995, Whittaker et al., 2003, Hoffman et al., 2011 and the coefficients of geometric distribution (P) from 0.1 to 0.7. Motif sizes and alleles ranges followed the empirical data of each locus. The summary statistics were the mean number of alleles (A), genetic diversity (H E ), allelic size range (AR), and Garza-Williamson's index (M). A total of 3,000,000 simu-lations were performed to generate the reference table, using the four scenarios according to their prior probability and their parameter values drawn from the prior distributions. The posterior probability of each scenario was assessed using both direct estimate and logistic regression approaches using between 500 and 30,000 best simulations. Under an ABC approach, the best scenario is the one with the simulated data set closest to observed data set. For the best scenario, the posterior distribution of the parameters was estimated using logit transformation for the 8,000 best simulations.
We converted the effective sizes (N e ) obtained in each method to the census sizes (N c ) as follows. First, the ratio of mature adults to the effective number of adults (N T :N e ) approaches 2 for most populations, and was based on the equation N e = N/(2 -T -1) from Nunney and Elam (1994), where T is the generation length. Second, the proportion of juveniles in the population (number of adults + juveniles)/(number of adults), was estimated between 1.6 to 2.0 256 Cypriano-Souza et al. for humpback whales based on catch and survey data (Chittleborough, 1965, Roman andPalumbi, 2003). Therefore, multiplying the two ratios, the average ratio of census population size to effective population size was 3.6, with a variation from 3.2 to 4.0, which has also been used in previous studies (Roman and Palumbi, 2003;Alter et al., 2007Alter et al., , 2012Ruegg et al., 2010Ruegg et al., , 2013.

Genetic variability
Individual multilocus genotypes were on average 98.5% complete. Summary statistics are presented in Table  2. The number of alleles identified at the 10 microsatellite loci ranged from five (EV1) to 18 (GATA417), with a mean of 12.6. The mean observed (H O ) heterozygosity was 0.736, ranging from 0.553 (EV1) to 0.923 (GATA417), and the mean expected (H E ) heterozygosity was 0.746, ranging from 0.532 (EV1) to 0.923 (EV37). Population-wide F IS values were low for all loci (below 0.05), except for the locus GATA053 (F IS = 0.053) and the locus 417/418 (F IS = 0.068), but these values were not significant. Moreover, no evidence of null alleles and no significant deviation from HWE expectations were seen at any of the loci. Pairwise comparison of allele frequencies revealed no significant linkage disequilibrium after Bonferroni correction.
In the comparison of the four scenarios (constant population, bottlenecked, expanded, and transitory bottleneck) using the ABC approach implemented in DIYABC, the constant population (no demographic changes) scenario was highly supported (posterior probability > 0.99) in relation to the other scenarios with demographic changes, in both the direct estimate and logistic regression approaches (Figure 2b,c). In the constant scenario, the mode of the posterior distribution for the long term N e was 4,170 (95% CI = 2,330 -26,600) ( Figure 3). The mode for the N c was then 15,012 (95% CI = 8,388 -95,760), using the 3.6 census/effective ratio.
Humpback whale genetic diversity 257

Discussion
Our extended sampling confirms our previous results on the high nuclear DNA diversity of the humpback whale population that winters off the Brazilian coast (BSA), which is also compatible with its high mtDNA variability (Engel et al., 2008;Cypriano-Souza et al., 2010. This high genetic diversity is in agreement with other breeding grounds studied in the Southern Hemisphere for both nuclear and mitochondrial markers (e.g. Valsecchi et al., 2002;Garrigue et al., 2004;Pomilla and Rosenbaum, 2006;Olavarría et al., 2007;Rosenbaum et al., 2009;Cypriano-Souza et al., 2017). Nevertheless, in the Southwestern Atlantic population, as well as most other humpback populations, severe reductions of their historical size by commercial whaling are well documented. The lowest number reached for breeding stock A was in the late 1950s, when only around 500 individuals (95% CI = 152 to 3,687) were estimated for this population (Zerbini et al., 2011b). However, our previous study did not detect any significant signal of a genetic bottleneck in this population (Cypriano-Souza et al., 2010), using three standard methods: heterozygosity excess (Cornuet and Luikart, 1996), mode-shift  and M-ratio tests (Garza and Williamson, 2001). That result was corroborated here with an extended data set and an ABC approach, in which by far the best supported scenario was a constant population compared with those in which a population experienced a single size change (expansion or bottleneck) or a bottleneck during the commercial whaling (between 2 and 8 generations ago) followed by an expansion (Figure 2).
As discussed previously, these results are consistent with the hypothesis (Amos, 1996) that the genetic bottleneck caused by commercial whaling was not strong enough to have left a significant signal in the BSA population (Engel et al., 2008;Cypriano-Souza et al., 2010). The magnitude of the genetic bottleneck is related to its duration and the minimum absolute size of the population during the bottleneck (Frankham et al., 2002). In this population, large scale whaling lasted for only about four generations, assuming that the overexploitation of this stock was in the period between 1904 and 1967 (Paiva and Grangeiro, 1970) and a generation time of 18 years, taking into account the range between 12 and 24 years estimated for humpback whales (Chittleborough, 1965;Roman and Palumbi, 2003). In addition, the minimum absolute population size reached was relatively large (N min = 500 individuals in the late 1950s) for the BSA population (Zerbini et al., 2011b). Therefore, fewer generations with a not so small absolute effective population size should have left only weak genetic bottleneck signals and are, therefore, more difficult to detect. However, this estimation based on catch records did not incorporate missing whaling records between 1929 and 1946, a period in which the catch records in the breeding grounds are incomplete, producing likely biased estimates of depletion levels for this population (Morais et al., 2017).
Recently, Phillips et al. (2012) showed a similar result with an ABC analysis for the bowhead whales (Balaena mysticetus), in which a bottleneck scenario was also not supported. In contrast, a recent ABC analysis for the Antarctic fur seal (Arcthocephalus gazella) supported a bottleneck scenario, although they have not detected a bottleneck using standard tests (Hoffman et al., 2011). These different results were expected given the contrasting values for the two key parameters discussed above between the bowhead whale and the Antarctic fur seal: the generation length of the former was estimated~50 years and for the latter~10 years, while the minimum population size during the bottleneck for the former was around 1,000 individuals and for the latter it was as low as~30-60 individuals (Phillips et al., 2012). Interestingly, the inclusion of ancient samples may importantly increase the power to detect recent bottlenecks, as demonstrated by the study with the eastern Pacific gray whales (Eschrichtius robustus) (Alter et al., 2012).
Although the methods used here assumed closed population with discrete generations (no generation overlap), random mating and equal contribution of individuals to the next generation, these assumptions are rarely completely satisfied in natural populations. The humpback whale is a long-lived species with overlapping generations, and some migration between breeding grounds in the Southern Hemisphere has been documented (Rosenbaum et al., 2009). However, a recent study indicates that single-sample estimators of contemporary N e based on linkage disequilibrium were not affected significantly with migration rates up to approximately 5-10% (Waples and England, 2011). Indeed, the migration rates between the Atlantic and Pacific breeding grounds of South America were estimated to be lower than the above limits (Cypriano-Souza et al., 2017). Likewise, low migration rates based on mtDNA data have been estimated between humpbacks from Brazil and from breeding grounds in the Southeastern Atlantic Ocean (H.C. Rosenbaum, personal communication). The humpback whale is an age-structured species, and our estimates based on mixed-age samples reflect the effective size per generation (N e ), but these estimates are approximately equal to population N e when the sample includes as many cohorts as there are in a generation . As we sampled only mature adults from nine different breeding seasons, the number of cohorts in this sample is lower than the generation length. However, our estimates should be fairly robust, albeit perhaps slightly lower, since the species is long-lived and has intermittent breeding .
Our N e estimates with NeEstimator and DIYABC (point estimates 1,078 and 4,170 individuals, respectively) differed significantly. These differences can be explained mostly because these methods are expected to estimate N e on different periods. The estimator from NeEstimator uses the amount of linkage disequilibrium within a population 258 Cypriano-Souza et al.
and is specifically designed to estimate N e from the parental generation of the sample. In contrast, the N e inferred from the ABC approach, considering the constant population scenario, is a long-term (that extends well before the whaling period) population size estimator, that represented the weighted harmonic mean of population sizes over 4N e generations.
The most recent contemporary abundance estimates of the Brazilian humpback whale population, derived from aerial surveys that covered the entire stock range, estimated 6,404 individuals (95% CI = 5,085-8,068) in 2005, 7,689 individuals (95% CI = 6,585-8,931) in 2008, 8,652 individuals (95% C.I. = 7,696-9,682) in 2011, and 12,123 individuals (95% C.I. = 10,811-13,531) in 2015 (Andriolo et al., 2010;Pavanato et al., 2017). However, recent abundance estimates derived from ship surveys estimated 16,410 individuals (95% C.I. = 10,563-25,495) in 2008 (covered the total population range), and 19,429 individuals (95% C.I. = 15,958-23,654) in 2012, which were about 50% higher than that of the aerial surveys (Bortolotto et al., 2016a). The difference between aerial and ship abundance estimates can be explained mainly by two sources of bias (underestimated group sizes and perception bias) of the first methodology, which may bias downwards the aerial abundance estimates. Therefore, assuming that our census size (N c ) estimated by the NeEstimator, between around 2,600 and 6,800, corresponds to the parental generation of our sample (collected between 1999 and 2007), which given the uncertainty in the generation time would roughly correspond to the years between around 1980 and 1990. These values are broadly compatible with the abundance ship surveys extrapolated to the past using the growth rate of 7.4% per annum (Ward et al., 2011).
The study on historical (pre-whaling) abundance of the stock A based on catch records using a Bayesian statistical method estimated the population size to nearly 24,700 individuals (95% CI = 22,804-31,220) before exploitation by modern whaling (Zerbini et al., 2011b). Our point estimate for the long term N c (~15,000) using the ABC approach was smaller than the pre-whaling abundance cited above, although the confidence intervals widely overlap that estimate. Furthermore, it is compatible with the most representative abundance estimate of approximately 16,000 humpback whales derived from shipboard survey in 2008 (Bortolotto et al., 2016a). However, our results present broad confidence intervals usually derived from large uncertainties of several parameters, such as generation time, mutation rate, or the relation between N c and N e . Consequently, more loci, more realistic scenarios (non-instantaneous population growth, gene flow, etc.) and methods are necessary to better estimate the demographic history of this population.
Overall, our results corroborate previous studies that have found high genetic diversity in the humpback whale population wintering off the Brazilian coast and no statisti-cally significant reduction in its genetic diversity caused by modern whaling (Engel et al., 2008;Cypriano-Souza et al., 2010). Compatible with that, and estimated here for the first time, the N e of the Southwestern Atlantic Ocean humpback whale population is relatively large. It has been suggested that the short-term minimum viable population size is N e > 50 and the long-term minimum viable population size is N e > 500 (Franklin and Frankham, 1998). Therefore, the contemporary N e estimates of the Brazilian humpback population obtained here were above of these limits, suggesting good protection for its evolutionary future. However, we suggest that similar population size estimates should be carried out every 12 to 24 years (one generation) in order to obtain comparable N e estimates and to monitor this population. To better accomplish this goal, future studies should attempt to reduce uncertainties of several key parameters, for example increasing the number of loci.