Molecular Ecology (2015)
doi: 10.1111/mec.13092
Heterozygosity–fitness correlations in a declining seabird population L V A R O B A R R O S * and P A L O M A M O R A N † ALBERTO VELANDO,* A *Departamento de Ecoloxıa e Bioloxıa Animal, Universidade de Vigo, Campus As Lagoas, 36310 Vigo, Spain, †Departamento de Bioquımica, Xenetica e Inmunoloxıa, Universidade de Vigo, Campus As Lagoas, 36310 Vigo, Spain
Abstract Loss of genetic diversity is thought to lead to increased risk of extinction in endangered populations due to decreasing fitness of homozygous individuals. Here, we evaluated the presence of inbreeding depression in a long-lived seabird, the European shag (Phalacrocorax aristotelis), after a severe decline in population size by nearly 70%. During three reproductive seasons, 85 breeders were captured and genotyped at seven microsatellite loci. Nest sites were monitored during the breeding season to estimate reproductive success as the number of chicks surviving to full-size-grown per nest. Captured birds were tagged with a ring with an individual code, and resighting data were collected during 7-year period. We found a strong effect of multilocus heterozygosity on female reproductive performance, and a significant, although weaker, effect on breeder survival. However, our matrix population model suggests that this relatively small effect of genetic diversity on breeder survival may have a profound effect on fitness. This highlights the importance of integrating life history consequences in HFC studies. Importantly, heterozygosity was correlated across loci, suggesting that genomewide effects, rather than single loci, are responsible for the observed HFCs. Overall, the HFCs are a worrying symptom of genetic erosion in this declining population. Many long-lived species are prone to extinction, and future studies should evaluate the magnitude of fitness impact of genetic deterioration on key population parameters, such as survival of breeders. Keywords: birds, bottleneck, fitness, genetic diversity, inbreeding, population dynamics Received 5 October 2014; revision received 20 January 2015; accepted 21 January 2015
Introduction Loss of genetic diversity is a major conservation concern in small, isolated (i.e. endangered) populations. Genetic drift and inbreeding often lead to lower individual fitness and, hence, increased risk of extinction (Soule 1987; H€ oglund 2009). Inbreeding increases the expression of recessive deleterious mutations and the loss of heterozygote advantage (Charlesworth & Charlesworth 1987). In natural populations, estimates of individual inbreeding levels have been quantified by measuring multilocus heterozygosity (MHL) at neutral makers (Coltman & Slate 2003; Aparicio et al. 2006; Chapman et al. 2009). A number of previous studies Correspondence: Alberto Velando, Fax: +34986812556; E-mail: [email protected] © 2015 John Wiley & Sons Ltd
have reported a relationship between MHL and fitnessrelated traits (so-called heterozygosity–fitness correlations; HFCs), although the relationship reported in these studies is usually relatively weak (reviewed in Chapman et al. 2009). Nevertheless, the significance of observed HFCs as a proxy for inbreeding depression has been the subject of intense debate in the last decade (e.g. Balloux et al. 2004; Miller & Coltman 2014). Heterozygosity–fitness correlations may be brought about by mechanisms other than genomewide heterozygosity (Hansson & Westerberg 2002). HFCs may arise if some of the scored loci are under selection and heterozygous genotypes at these loci have a fitness advantage (‘direct effect’; David 1998; Hansson & Westerberg 2002); however, this effect is not expected in noncoding DNA, such as microsatellites (Jarne & Lagoda 1996). Although neutral markers may not have direct effects
. B A R R O S and P . M O R A N 2 A. VELANDO, A on fitness, they may be associated by linkage disequilibrium to loci under selection (Lynch & Walsh 1998). Hence, an apparent heterozygote advantage may result from specific loci influencing fitness (‘local effect’; David 1998; Hansson & Westerberg 2002), rather than from effects of overall homozygosity in the genome (‘general effect’ Hansson & Westerberg 2002; Keller & Waller 2002). For genomewide effects, a correlation in heterozygosity across loci within individuals (i.e. identity disequilibrium, ID) generated by inbreeding is expected (Szulkin et al. 2010; Kardos et al. 2014). Thus, the magnitude of identity disequilibrium may be used as an indication of the general effects. A recent metaanalysis showed that the average effect size of HFCs was associated with the genomewide identity disequilibrium observed in the population (Miller & Coltman 2014). The demographic history of a population (such as the occurrence of genetic drift, population admixture or bottlenecks) is expected to influence genomewide heterozygosity (Szulkin et al. 2010). Most HFC studies have been conducted in large outbred populations, but their significance in small bottlenecked populations has been questioned (Grueber et al. 2008). Neutral markers in populations recently affected by a large reduction in size are likely to reflect genomewide heterozygosity (Bierne et al. 2000; Szulkin et al. 2010). Inbreeding depression also increases under stressful conditions (Jimenez et al. 1994; Fox & Reed 2011), but it may be difficult to detect by HFCs because genetic diversity is probably also reduced in bottlenecked populations (Frankham et al. 2002; Grueber et al. 2008). Additionally, populations recently suffering severe reductions may display a higher degree of linkage disequilibrium than large, stable populations (Hartl & Clark 1997; Wang et al. 1998). The presence of linkage disequilibrium increases the probability of local effects (Brouwer et al. 2007), making HFCs difficult to interpret (Grueber et al. 2008). Thus, the relevance of HFCs in threatened populations must still be evaluated (but see Forcada & Hoffman 2014). Here, we examined the occurrence of inbreeding depression in a long-lived seabird, the European shag (Phalacrocorax aristotelis), after a large decrease in the size of the population (Rias Baixas, Galicia, NW Spain) following an oil spill (see Fig. S2, Supporting information). A previous European scale study of the population genetic structure of the species identified a genetic cluster at the southern limit (including the Rıas Baixas population; Barlow et al. 2011), reflecting limited levels of contemporary gene flow due geographic constraints. Preliminary analysis of dispersal patterns also suggested that breeding shags in the Rıas Baixas population are currently isolated from other populations (Velando & Alvarez 2004). This coastal seabird is a
year-round local resident in the study area (Velando & Munilla 2008) and is highly philopatric, with most of the birds recruiting within 5 km of their natal site (Velando & Freire 2002; Barros et al. 2013). This population of European shags is therefore at risk of inbreeding due to its population structure, as with other inshore seabirds that display similar life characteristics (see Friesen et al. 2007). Indeed, a previous study suggested that inbreeding occurs in this population, as revealed by homozygosity in microsatellite data (Barlow et al. 2011). On 13 November 2002, the Prestige oil tanker was wrecked and about 63 000 tonnes of heavy oil were released into the marine environment. The European shag population in the Rıas Baixas was strongly affected by the initial spillage (Velando et al. 2005a,b; Martınez-Abraın et al. 2006), with long-term consequences at population level (Velando & Munilla 2008; Barros et al. 2014). In addition to direct mortality, longterm consequences have been attributed to stressful conditions after the oil spill, including delayed effects in response to sublethal oil exposure and reduced food availability (Velando et al. 2005a, 2010; Velando & Munilla 2008). Prior to the Prestige oil spill, the population was declining at an annual rate of 5% (Velando & Freire 2002), and it was classified as ‘endangered’ in the Red Book of Birds of Spain (Velando & Alvarez 2004). Five years after the Prestige oil spill, the breeding population of the shag population in the Rıas Baixas was about 70% lower than prespill counts (Alvarez & Velando 2007, see Fig. S2, Supporting information). In this study, we investigate the effect of multilocus heterozygosity on reproductive output and adult survival (using capture–mark–recapture analyses). Here, we used genetic and population data to explore the relationships between heterozygosity and two fitnessrelated traits: reproductive success and adult survival. We estimated a locus-weighted measure of homozygosity (HL; Aparicio et al. 2006) based on seven microsatellite markers. Under a general effect, we predicted a negative relationship between HL and reproductive success and/or survival. Finally, we investigated identity disequilibrium in the population by determining whether heterozygosity is correlated across loci, as expected by genomewide effects of inbreeding (Szulkin et al. 2010; Kardos et al. 2014).
Methods Study area and sampling This study was carried out in the breeding colonies of the European shag in the Rıas Baixas (Galicia, Spain) between 2004 and 2013. This area holds the largest population (80% of the total) of European shags in the © 2015 John Wiley & Sons Ltd
HETEROZYGOSITY AND FITNESS IN A SEABIRD 3 Atlantic Iberian Peninsula (Velando & Alvarez 2004). Breeding colonies are found in the Illas Atl anticas National Park, including three groups of islands (Illas Cıes, Illa de Ons and Sagres; see Fig. 1 in Barros et al. 2013). During three breeding seasons, 85 incubating adult birds were captured by hand (20 individuals in 2004, 54 in 2007 and 11 in 2012, see Table S1, Supporting information), and blood samples were obtained from each bird. The birds were tagged with a numbered metal ring and a coloured plastic ring with an individual two-digit combination to facilitate identification at a distance. Blood samples were stored at low temperature (4ªC) until they were centrifuged, on the same day. The blood cells separated by centrifugation were stored frozen at 80 °C.
DNA extraction and microsatellite analysis Total genomic DNA from red blood cells was extracted by a chaotropic NaI-based method (Gedik & Collins 2005). Sex was identified from blood cell DNA by PCR amplification of two CHD genes (CHD1W and CHD1Z), using a pair of primers (2550F and 2718R; Fridolfsson & Ellegren 1999). All samples were genotyped at three dinucleotide (Phaari02, Phaari05 and Phaari06) and four tetranucleotide (Phaari08, Phaari11, Phaari12 and Phaari16) microsatellite loci (see Table S2, Supporting information, described in Molecular Ecology Resources Primer Development Consortium 2010; Barlow et al. 2011). The 50 end of each forward primer was fluorescently labelled with either 6-FAM, HEX or NED, and genotypes were resolved on an automatic ABI 3100 capillary DNA sequence. Genotypes were scored using GENEMAPPER (Applied Biosystems). The data set was checked for genotyping errors, and null-allele frequencies per population were estimated using MICROCHECKER 2.2.3 (van Oosterhout et al. 2004).
Heterozygosity variation Observed (Ho) and expected (He) heterozygosity at each locus (Table 1) were calculated using GENALEX 6.4 (Peakall & Smouse 2006). In a Markov chain Monte Carlo (MCMC) approach (1000 dememorizations, 100 batches, 1000 iterations), GENEPOP 4.0.10 (Raymond & Rousset 1995; Rousset 2008) was used to test for deviations from Hardy–Weinberg equilibrium (HWE) per locus in each sampling year (Table 1, Fisher’s method). We also evaluated the genetic signature of the major reduction in size experienced by the Rıas Baixas shag population (see Supplementary material). As rare alleles, which contribute marginally to heterozygosity, are more likely to be lost following a bottleneck, a tran© 2015 John Wiley & Sons Ltd
sient excess in heterozygosity can be used as a genetic signature of recent bottlenecks (Luikart & Cornuet 1998; Peery et al. 2012). Individual multilocus heterozygosity was estimated by calculating the homozygosity weighted by locus (HL), which was determined from the multilocus genotypes of 7 microsatellite loci using CERNICALIN (Aparicio et al. 2006). HL improves heterozygosity estimates by weighting the contribution of each locus in relation to the overall homozygosity value (Rijks et al. 2008).
Reproduction Nest sites of 66 captured breeding shags (38 males and 29 females) were monitored three to five times during the breeding season (March to June) in the year of capture (2004, 2007, 2012). Nest sites were marked with epoxy resin. In all cases, nests were monitored at least once during the incubation period and twice during the chick-rearing period. This procedure enabled estimation of reproductive success as the number of chicks surviving to full-size-grown per nest (age >35 days, Velando et al. 2000; Barros et al. 2014). Reproductive success was analysed by a generalized linear model for count data with Poisson errors and a natural logarithm link. The model included sampling year, island (Illas Cıes, Illa de Ons and Sagres) and HL. We also explored the quadratic effect of HL on reproductive success, but this nonlinear effect was not significant. To test sex-specific effects of HL on reproductive success, the sex of the captured adult and its interaction with HL were included in the model. Note that extra-pair paternity is frequent in this species (Graves et al. 1993), and therefore, reproductive success cannot be confidently estimated in males without paternity assignment.
Capture–recapture and survival estimation Since our resighting scheme started in 2007, survival was estimated on adult birds captured in 2007 and 2012 (see Table S1, Supporting information). During the period 2007–2013, we collected resighting data by intensive field monitoring of marked shags in the Parque Nacional das Illas Atl anticas (all colonies were visited more than three times per year in 2007, 2008, 2009, 2012 and 2013, and once per year in 2010, 2011). In our sampling localities, European shags with colour rings are easily detected and the probability of resighting is relatively high (Velando 2000; Velando & Freire 2002; Noguera et al. 2012; Barros et al. 2013). Adult shags exhibit strong fidelity to the colony where they bred for first time, and adult dispersal rates are therefore very low (Aebischer et al. 1995). This study is part of a large capture–recap-
. B A R R O S and P . M O R A N 4 A. VELANDO, A Table 1 Characteristics of the seven microsatellites used in this study 2004 (n = 20)
2007 (n = 54)
2012 (n = 11)
Locus
Size range
Na
Ho
He
Ho
He
Ho
He
Phaari02 Phaari05 Phaari06 Phaari08 Phaari11 Phaari12 Phaari16
97–109 156–158 207–231 153–189 120–192 233–265 186–226
3 6 3 10 14 8 14
0.100 0.800 0.450 0.900 0.800 0.950 0.800
0.095 0.585 0.399 0.809 0.884 0.816 0.853
0.074 0.796 0.000 0.852 0.704 0.759 0.667
0.071 0.730 0.071 0.757 0.880 0.790 0.873
0.182 1.000 0.091 0.909 0.636 0.818 0.909
0.165 0.640 0.087 0.789 0.769 0.702 0.831
Na, number of alleles; Ho, observed Heterozygosity; He, unbiased expected heterozygosity. In bold, significant deviations in expected genotype frequencies from Hardy–Weinberg expectation, at 0.05 significance level.
ture survey with more than 2500 resightings across the Atlantic coast of Iberian Peninsula (see Barros et al. 2013). Study birds were resighted 67 times and none at breeding sites other than the natal colony. We used the MARK program (White & Burnham 1999; White 2000) and Cormack–Jolly–Seber models (Lebreton et al. 1992) to estimate the probability of apparent local survival and for the recapture probability. We first selected an ‘umbrella model’, that is the most parsimonious model, from a starting model that included the relevant environmental and biological data in the study system. As expected, survival was not affected by sampling island and recapture probability was similar in both sexes (Table 2). Thus, we modelled apparent survival probability (Ф) as a function of sex and year and recapture probability (p) as function of year and colony (Ф(sex + t), p(island + t)). The inflation factor (^c) in the starting model, estimated by performing 1000 bootstrap simulations, indicated that this model showed a negligible overdispersion (^c = 1.19), and the bootstrap goodness of fit test confirmed that the starting model adequately fitted the data (P = 0.164). After model selection for the umbrella model, we included HL. Quadratic term of HL did not improve the final model. We used the Akaike information criterion (Burnham & Anderson 2002) for model selection and ranking of competing models.
Fitness and selection To compare the HFCs of each vital rate, fitness and selection were estimated in a matrix-modelling framework (Lande 1982). The intrinsic population growth rate (k) was estimated from the dominant eigenvalue of the Leslie matrix (Caswell 2001), incorporating all the vital rates at all the age classes of the life cycle. The sensitivity index of k may be interpreted as directional selection (Lande 1982; Caswell 2001). As reproductive success cannot be confidently estimated in males,
Table 2 Summary of the model selection statistics for the umbrella model (see Methods) and models incorporating homozygosity weighted by locus (HL) on survival probabilities. Models were ranked according to their AICc value, with the best supported model in bold Model
AICc
DAIC
(a) umbrella model 264.60 0.000 Ф(sex), p(t) 265.63 1.028 Ф(.), p(t) 266.19 1.590 Ф(sex + t), p(t) 267.54 2.939 Ф(t), p(t) 268.69 4.086 Ф(island), p(t) 270.10 5.50 Ф(sex + t), p(island + t) 279.56 14.95 Ф(sex + t), p(sex + t) 281.73 17.12 Ф(island + t), p(t) 282.52 17.92 Ф(sex + t), p(.) 286.54 21.94 Ф(sex + t), p(island) 287.48 22.88 Ф(sex + t), p(sex) (b) models incorporating genetic terms 259.83 0.000 Ф(sex + HL), p(t) 261.26 1.427 Ф(sex*HL), p(t) 264.60 4.769 Ф(sex), p(t)
AICc Weight
K
0.393 0.235 0.178 0.090 0.051 0.025 0.000 0.000 0.000 0.000 0.000
8 7 10 9 9 12 18 16 5 7 8
0.632 0.310 0.058
9 10 8
intrinsic population growth rates were only estimated in females. A stage (age) matrix model (Fig. S3, Supporting information) was built according to shag population dynamics in the study population (Velando & Freire 2002; Velando et al. 2005a; Martınez-Abraın et al. 2006). Thus, the matrix model included four age classes (juveniles, first-year olds, second-year olds and adults). A balanced chick sex ratio was assumed in the model (Velando et al. 2002). In this model, adult female survival (Sa, Fig. S3, Supporting information) and reproduction (f, Fig. S3, Supporting information) and their variance were dependent on HL, according to our estimation based on genotyped females (see results). The remaining demographic parameters were established in accordance with © 2015 John Wiley & Sons Ltd
H E T E R O Z Y G O S I T Y A N D F I T N E S S I N A S E A B I R D 11 Keller LF, Waller DM (2002) Inbreeding effects in wild populations. Trends in Ecology and Evolution, 17, 230–241. Lande R (1982) A quantitative genetic theory of life-history evolution. Ecology, 63, 607–615. Lebreton JD, Burnhan KP, Clobert J, Anderson DR (1992) Modelling survival and testing biological hypotheses using marked animals: a unified approach with case studies. Ecological Monographs, 62, 67–118. Lieutenant-Gosselin M, Bernatchez L (2006) Local heterozygosity-fitness correlations with global positive effects on fitness in threespine stickleback. Evolution, 60, 1658–1668. Luikart G, Cornuet JM (1998) Empirical evaluation of a test for identifying recently bottlenecked populations from allele frequency data. Conservation Biology, 12, 228–233. Lynch M, Walsh B (1998) Genetics and Analysis of Quantitative Traits. Sinauer Associates Inc, Sunderland, Massachusetts. Martınez-Abraın A, Velando A, Genovart M et al. (2006) Sexspecific mortality of European shags during an oil spill: demographic implications for the recovery of colonies. Marine Ecology Progress Series, 318, 271–276. Miller JM, Coltman DW (2014) Assessment of identity disequilibrium and its relation to empirical heterozygosity fitness correlations: a meta-analysis. Molecular Ecology, 23, 1899– 1909. Molecular Ecology Resources Primer Development Consortium (2010) Permanent genetic resources added to molecular ecology resources database 1(April), pp. 2010–31, May 2010. Molecular Ecology Resources, 10, 1098–1105. Noguera JC, Kim S-Y, Velando A (2012) Pre-fledgling oxidative damage predicts recruitment in a long-lived bird. Biology Letters, 8, 61–63. Olano-Marin J, Mueller JC, Kempenaers B (2011) Heterozygosity and survival in blue tits (Cyanistes caeruleus): contrasting effects of presumably functional and neutral loci. Molecular Ecology, 20, 4028–4041. van Oosterhout C, Hutchinson WF, Wills DP, Shipley P (2004) MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Molecular Ecology Notes, 4, 535–538. Ortego J, Calabuig G, Cordero PJ, Aparicio JM (2007) Egg production and individual genetic diversity in lesser kestrels. Molecular Ecology, 16, 2383–2392. Peakall R, Smouse PE (2006) GENALEX 6: genetic analysis in Excel. Population genetic software for teaching and research. Molecular Ecology Notes, 6, 288–295. Peery MZ, Kirby R, Reid BN et al. (2012) Reliability of genetic bottleneck tests for detecting recent population declines. Molecular Ecology, 21, 3403–3418. Pemberton J (2004) Measuring inbreeding depression in the wild: the old ways are the best. Trends in Ecology & Evolution, 19, 613–615. Raymond M, Rousset F (1995) GENEPOP Version 1.2: population genetics software for exact tests and ecumenicism. Journal of Heredity, 86, 248–249. Rijks JM, Hoffman JI, Kuiken T, Osterhaus ADME, Amos W (2008) Heterozygosity and lungworm burden in harbour seals (Phoca vitulina). Heredity, 100, 587–593. Robinson SP, Simmons LW, Kennington WJ (2013) Estimating relatedness and inbreeding using molecular markers and pedigrees: the effect of demographic history. Molecular Ecology, 22, 5779–5792.
© 2015 John Wiley & Sons Ltd
Rousset F (2008) Genepop’007: a complete reimplementation of the Genepop software for Windows and Linux. Molecular Ecology Resources, 8, 103–106. Samollow PB, Soule ME (1983) A case of stress related heterozygote superiority in nature. Evolution, 37, 646–649. Seddon N, Amos W, Mulder RA, Tobias JA (2004) Male heterozygosity predicts territory size, song structure and reproductive success in a cooperatively breeding bird. Proceedings of the Royal Society B, 271, 1823–1829. Slate J, David P, Dodds KG et al. (2004) Understanding the relationship between the inbreeding coefficient and multilocus heterozygosity: theoretical expectations and empirical data. Heredity, 93, 255–265. Soule ME (1987) Viable Populations for Conservation. Cambridge University Press, Cambridge. Szulkin M, Bierne N, David P (2010) Heterozygosity-fitness correlations: a time for reappraisal. Evolution, 64, 1202–1217. Tiira K, Laurila A, Enberg K et al. (2006) Do dominants have higher heterozygosity? Social status and genetic variation in brown trout, Salmo trutta. Behavioral Ecology & Sociobiology, 59, 657–665. V€ alim€ aki K, Hinten G, Hanski I (2007) Inbreeding and competitive ability in the common shrew (Sorex araneus). Behavioral Ecology & Sociobiology, 61, 997–1005. Velando A (2000) The importance of the hatching date on dominance of young shags. Animal Behaviour, 60, 181–185. Velando A, Alvarez D (2004) Cormoran Mo~ nudo, Phalacrocorax aristotelis aristotelis. In: Libro Rojo de las Aves de Espa~ na (eds Madro~ no A, Gonzalez C, Atienza JC), pp. 60–62. Ministerio de Medio Ambiente-SEO/Birdlife, Madrid. Velando A, Freire J (2002) Population modelling of European shag at their southern limit: conservation implications. Biological Conservation, 107, 59–69. Velando A, Munilla I (2008) Plan de Conservacion del Cormoran Mo~ nudo en el Parque Nacional de las Islas Atlanticas. Universidade de Vigo, Vigo. Velando A, Ortega-Ruano JE, Freire J (1999) Chick mortality in European shag Stictocarbo aristotelis related to food limitations during adverse weather events. Ardea, 87, 51–59. Velando A, Graves J, Freire J (2000) Sex-specific growth in the European Shag Stictocarbo aristotelis, a seabird with size dimorphism. Ardea, 88, 127–136. Velando A, Graves J, Ortega-Ruano JE (2002) Sex ratio in relation to timing of breeding, and laying sequence in a dimorphic seabird. Ibis, 144, 9–16. Velando A, Munilla I, Leyenda PM (2005a) Short-term indirect effects of the ˋPrestigeˊ oil spill on European shags: changes in availability of prey. Marine Ecology Progress Series, 302, 263–27. (2005b) Velando A, Alvarez D, Mouri~ no J, Arcos F, Barros A Population trends and reproductive success of the European shag Phalacrocorax aristotelis on the Iberian Peninsula following the Prestige oil spill. Journal of Ornithology, 146, 116–120. Velando A, Munilla I, L opez-Alonso M, Freire J, Perez C (2010) EROD activity and stable isotopes in seabirds to disentangle marine food web contamination after the Prestige oil spill. Environmental Pollution, 158, 1275–1280. Wang JL, Caballero A, Hill WG (1998) The effect of linkage disequilibrium and deviation from Hardy-Weinberg proportions on the changes in genetic variance with bottlenecking. Heredity, 81, 174–186.
. B A R R O S and P . M O R A N 6 A. VELANDO, A d.f. = 1, P = 0.80) or years (HL*year interaction, Wald v2 = 0.34, d.f. = 1, P = 0.56).
Incorporation of HL in the model explained 2.88% of the residual deviance.
HL and adult survival
HL and fitness
The umbrella model (Table 2a) indicated that recapture probability differed between years (DAIC = 16.33 from a model with time-constant recapture, Fig. S3, Supporting information). Recapture probability was not affected by sex or sampling island (Table 2a). The inclusion of survival differences among sampling islands did not improve the model (Table 2a). The model of survival probability with the lowest AICc included differences between sexes (Table 2a), but this model differed by
doi: 10.1111/mec.13092
Heterozygosity–fitness correlations in a declining seabird population L V A R O B A R R O S * and P A L O M A M O R A N † ALBERTO VELANDO,* A *Departamento de Ecoloxıa e Bioloxıa Animal, Universidade de Vigo, Campus As Lagoas, 36310 Vigo, Spain, †Departamento de Bioquımica, Xenetica e Inmunoloxıa, Universidade de Vigo, Campus As Lagoas, 36310 Vigo, Spain
Abstract Loss of genetic diversity is thought to lead to increased risk of extinction in endangered populations due to decreasing fitness of homozygous individuals. Here, we evaluated the presence of inbreeding depression in a long-lived seabird, the European shag (Phalacrocorax aristotelis), after a severe decline in population size by nearly 70%. During three reproductive seasons, 85 breeders were captured and genotyped at seven microsatellite loci. Nest sites were monitored during the breeding season to estimate reproductive success as the number of chicks surviving to full-size-grown per nest. Captured birds were tagged with a ring with an individual code, and resighting data were collected during 7-year period. We found a strong effect of multilocus heterozygosity on female reproductive performance, and a significant, although weaker, effect on breeder survival. However, our matrix population model suggests that this relatively small effect of genetic diversity on breeder survival may have a profound effect on fitness. This highlights the importance of integrating life history consequences in HFC studies. Importantly, heterozygosity was correlated across loci, suggesting that genomewide effects, rather than single loci, are responsible for the observed HFCs. Overall, the HFCs are a worrying symptom of genetic erosion in this declining population. Many long-lived species are prone to extinction, and future studies should evaluate the magnitude of fitness impact of genetic deterioration on key population parameters, such as survival of breeders. Keywords: birds, bottleneck, fitness, genetic diversity, inbreeding, population dynamics Received 5 October 2014; revision received 20 January 2015; accepted 21 January 2015
Introduction Loss of genetic diversity is a major conservation concern in small, isolated (i.e. endangered) populations. Genetic drift and inbreeding often lead to lower individual fitness and, hence, increased risk of extinction (Soule 1987; H€ oglund 2009). Inbreeding increases the expression of recessive deleterious mutations and the loss of heterozygote advantage (Charlesworth & Charlesworth 1987). In natural populations, estimates of individual inbreeding levels have been quantified by measuring multilocus heterozygosity (MHL) at neutral makers (Coltman & Slate 2003; Aparicio et al. 2006; Chapman et al. 2009). A number of previous studies Correspondence: Alberto Velando, Fax: +34986812556; E-mail: [email protected] © 2015 John Wiley & Sons Ltd
have reported a relationship between MHL and fitnessrelated traits (so-called heterozygosity–fitness correlations; HFCs), although the relationship reported in these studies is usually relatively weak (reviewed in Chapman et al. 2009). Nevertheless, the significance of observed HFCs as a proxy for inbreeding depression has been the subject of intense debate in the last decade (e.g. Balloux et al. 2004; Miller & Coltman 2014). Heterozygosity–fitness correlations may be brought about by mechanisms other than genomewide heterozygosity (Hansson & Westerberg 2002). HFCs may arise if some of the scored loci are under selection and heterozygous genotypes at these loci have a fitness advantage (‘direct effect’; David 1998; Hansson & Westerberg 2002); however, this effect is not expected in noncoding DNA, such as microsatellites (Jarne & Lagoda 1996). Although neutral markers may not have direct effects
. B A R R O S and P . M O R A N 2 A. VELANDO, A on fitness, they may be associated by linkage disequilibrium to loci under selection (Lynch & Walsh 1998). Hence, an apparent heterozygote advantage may result from specific loci influencing fitness (‘local effect’; David 1998; Hansson & Westerberg 2002), rather than from effects of overall homozygosity in the genome (‘general effect’ Hansson & Westerberg 2002; Keller & Waller 2002). For genomewide effects, a correlation in heterozygosity across loci within individuals (i.e. identity disequilibrium, ID) generated by inbreeding is expected (Szulkin et al. 2010; Kardos et al. 2014). Thus, the magnitude of identity disequilibrium may be used as an indication of the general effects. A recent metaanalysis showed that the average effect size of HFCs was associated with the genomewide identity disequilibrium observed in the population (Miller & Coltman 2014). The demographic history of a population (such as the occurrence of genetic drift, population admixture or bottlenecks) is expected to influence genomewide heterozygosity (Szulkin et al. 2010). Most HFC studies have been conducted in large outbred populations, but their significance in small bottlenecked populations has been questioned (Grueber et al. 2008). Neutral markers in populations recently affected by a large reduction in size are likely to reflect genomewide heterozygosity (Bierne et al. 2000; Szulkin et al. 2010). Inbreeding depression also increases under stressful conditions (Jimenez et al. 1994; Fox & Reed 2011), but it may be difficult to detect by HFCs because genetic diversity is probably also reduced in bottlenecked populations (Frankham et al. 2002; Grueber et al. 2008). Additionally, populations recently suffering severe reductions may display a higher degree of linkage disequilibrium than large, stable populations (Hartl & Clark 1997; Wang et al. 1998). The presence of linkage disequilibrium increases the probability of local effects (Brouwer et al. 2007), making HFCs difficult to interpret (Grueber et al. 2008). Thus, the relevance of HFCs in threatened populations must still be evaluated (but see Forcada & Hoffman 2014). Here, we examined the occurrence of inbreeding depression in a long-lived seabird, the European shag (Phalacrocorax aristotelis), after a large decrease in the size of the population (Rias Baixas, Galicia, NW Spain) following an oil spill (see Fig. S2, Supporting information). A previous European scale study of the population genetic structure of the species identified a genetic cluster at the southern limit (including the Rıas Baixas population; Barlow et al. 2011), reflecting limited levels of contemporary gene flow due geographic constraints. Preliminary analysis of dispersal patterns also suggested that breeding shags in the Rıas Baixas population are currently isolated from other populations (Velando & Alvarez 2004). This coastal seabird is a
year-round local resident in the study area (Velando & Munilla 2008) and is highly philopatric, with most of the birds recruiting within 5 km of their natal site (Velando & Freire 2002; Barros et al. 2013). This population of European shags is therefore at risk of inbreeding due to its population structure, as with other inshore seabirds that display similar life characteristics (see Friesen et al. 2007). Indeed, a previous study suggested that inbreeding occurs in this population, as revealed by homozygosity in microsatellite data (Barlow et al. 2011). On 13 November 2002, the Prestige oil tanker was wrecked and about 63 000 tonnes of heavy oil were released into the marine environment. The European shag population in the Rıas Baixas was strongly affected by the initial spillage (Velando et al. 2005a,b; Martınez-Abraın et al. 2006), with long-term consequences at population level (Velando & Munilla 2008; Barros et al. 2014). In addition to direct mortality, longterm consequences have been attributed to stressful conditions after the oil spill, including delayed effects in response to sublethal oil exposure and reduced food availability (Velando et al. 2005a, 2010; Velando & Munilla 2008). Prior to the Prestige oil spill, the population was declining at an annual rate of 5% (Velando & Freire 2002), and it was classified as ‘endangered’ in the Red Book of Birds of Spain (Velando & Alvarez 2004). Five years after the Prestige oil spill, the breeding population of the shag population in the Rıas Baixas was about 70% lower than prespill counts (Alvarez & Velando 2007, see Fig. S2, Supporting information). In this study, we investigate the effect of multilocus heterozygosity on reproductive output and adult survival (using capture–mark–recapture analyses). Here, we used genetic and population data to explore the relationships between heterozygosity and two fitnessrelated traits: reproductive success and adult survival. We estimated a locus-weighted measure of homozygosity (HL; Aparicio et al. 2006) based on seven microsatellite markers. Under a general effect, we predicted a negative relationship between HL and reproductive success and/or survival. Finally, we investigated identity disequilibrium in the population by determining whether heterozygosity is correlated across loci, as expected by genomewide effects of inbreeding (Szulkin et al. 2010; Kardos et al. 2014).
Methods Study area and sampling This study was carried out in the breeding colonies of the European shag in the Rıas Baixas (Galicia, Spain) between 2004 and 2013. This area holds the largest population (80% of the total) of European shags in the © 2015 John Wiley & Sons Ltd
HETEROZYGOSITY AND FITNESS IN A SEABIRD 3 Atlantic Iberian Peninsula (Velando & Alvarez 2004). Breeding colonies are found in the Illas Atl anticas National Park, including three groups of islands (Illas Cıes, Illa de Ons and Sagres; see Fig. 1 in Barros et al. 2013). During three breeding seasons, 85 incubating adult birds were captured by hand (20 individuals in 2004, 54 in 2007 and 11 in 2012, see Table S1, Supporting information), and blood samples were obtained from each bird. The birds were tagged with a numbered metal ring and a coloured plastic ring with an individual two-digit combination to facilitate identification at a distance. Blood samples were stored at low temperature (4ªC) until they were centrifuged, on the same day. The blood cells separated by centrifugation were stored frozen at 80 °C.
DNA extraction and microsatellite analysis Total genomic DNA from red blood cells was extracted by a chaotropic NaI-based method (Gedik & Collins 2005). Sex was identified from blood cell DNA by PCR amplification of two CHD genes (CHD1W and CHD1Z), using a pair of primers (2550F and 2718R; Fridolfsson & Ellegren 1999). All samples were genotyped at three dinucleotide (Phaari02, Phaari05 and Phaari06) and four tetranucleotide (Phaari08, Phaari11, Phaari12 and Phaari16) microsatellite loci (see Table S2, Supporting information, described in Molecular Ecology Resources Primer Development Consortium 2010; Barlow et al. 2011). The 50 end of each forward primer was fluorescently labelled with either 6-FAM, HEX or NED, and genotypes were resolved on an automatic ABI 3100 capillary DNA sequence. Genotypes were scored using GENEMAPPER (Applied Biosystems). The data set was checked for genotyping errors, and null-allele frequencies per population were estimated using MICROCHECKER 2.2.3 (van Oosterhout et al. 2004).
Heterozygosity variation Observed (Ho) and expected (He) heterozygosity at each locus (Table 1) were calculated using GENALEX 6.4 (Peakall & Smouse 2006). In a Markov chain Monte Carlo (MCMC) approach (1000 dememorizations, 100 batches, 1000 iterations), GENEPOP 4.0.10 (Raymond & Rousset 1995; Rousset 2008) was used to test for deviations from Hardy–Weinberg equilibrium (HWE) per locus in each sampling year (Table 1, Fisher’s method). We also evaluated the genetic signature of the major reduction in size experienced by the Rıas Baixas shag population (see Supplementary material). As rare alleles, which contribute marginally to heterozygosity, are more likely to be lost following a bottleneck, a tran© 2015 John Wiley & Sons Ltd
sient excess in heterozygosity can be used as a genetic signature of recent bottlenecks (Luikart & Cornuet 1998; Peery et al. 2012). Individual multilocus heterozygosity was estimated by calculating the homozygosity weighted by locus (HL), which was determined from the multilocus genotypes of 7 microsatellite loci using CERNICALIN (Aparicio et al. 2006). HL improves heterozygosity estimates by weighting the contribution of each locus in relation to the overall homozygosity value (Rijks et al. 2008).
Reproduction Nest sites of 66 captured breeding shags (38 males and 29 females) were monitored three to five times during the breeding season (March to June) in the year of capture (2004, 2007, 2012). Nest sites were marked with epoxy resin. In all cases, nests were monitored at least once during the incubation period and twice during the chick-rearing period. This procedure enabled estimation of reproductive success as the number of chicks surviving to full-size-grown per nest (age >35 days, Velando et al. 2000; Barros et al. 2014). Reproductive success was analysed by a generalized linear model for count data with Poisson errors and a natural logarithm link. The model included sampling year, island (Illas Cıes, Illa de Ons and Sagres) and HL. We also explored the quadratic effect of HL on reproductive success, but this nonlinear effect was not significant. To test sex-specific effects of HL on reproductive success, the sex of the captured adult and its interaction with HL were included in the model. Note that extra-pair paternity is frequent in this species (Graves et al. 1993), and therefore, reproductive success cannot be confidently estimated in males without paternity assignment.
Capture–recapture and survival estimation Since our resighting scheme started in 2007, survival was estimated on adult birds captured in 2007 and 2012 (see Table S1, Supporting information). During the period 2007–2013, we collected resighting data by intensive field monitoring of marked shags in the Parque Nacional das Illas Atl anticas (all colonies were visited more than three times per year in 2007, 2008, 2009, 2012 and 2013, and once per year in 2010, 2011). In our sampling localities, European shags with colour rings are easily detected and the probability of resighting is relatively high (Velando 2000; Velando & Freire 2002; Noguera et al. 2012; Barros et al. 2013). Adult shags exhibit strong fidelity to the colony where they bred for first time, and adult dispersal rates are therefore very low (Aebischer et al. 1995). This study is part of a large capture–recap-
. B A R R O S and P . M O R A N 4 A. VELANDO, A Table 1 Characteristics of the seven microsatellites used in this study 2004 (n = 20)
2007 (n = 54)
2012 (n = 11)
Locus
Size range
Na
Ho
He
Ho
He
Ho
He
Phaari02 Phaari05 Phaari06 Phaari08 Phaari11 Phaari12 Phaari16
97–109 156–158 207–231 153–189 120–192 233–265 186–226
3 6 3 10 14 8 14
0.100 0.800 0.450 0.900 0.800 0.950 0.800
0.095 0.585 0.399 0.809 0.884 0.816 0.853
0.074 0.796 0.000 0.852 0.704 0.759 0.667
0.071 0.730 0.071 0.757 0.880 0.790 0.873
0.182 1.000 0.091 0.909 0.636 0.818 0.909
0.165 0.640 0.087 0.789 0.769 0.702 0.831
Na, number of alleles; Ho, observed Heterozygosity; He, unbiased expected heterozygosity. In bold, significant deviations in expected genotype frequencies from Hardy–Weinberg expectation, at 0.05 significance level.
ture survey with more than 2500 resightings across the Atlantic coast of Iberian Peninsula (see Barros et al. 2013). Study birds were resighted 67 times and none at breeding sites other than the natal colony. We used the MARK program (White & Burnham 1999; White 2000) and Cormack–Jolly–Seber models (Lebreton et al. 1992) to estimate the probability of apparent local survival and for the recapture probability. We first selected an ‘umbrella model’, that is the most parsimonious model, from a starting model that included the relevant environmental and biological data in the study system. As expected, survival was not affected by sampling island and recapture probability was similar in both sexes (Table 2). Thus, we modelled apparent survival probability (Ф) as a function of sex and year and recapture probability (p) as function of year and colony (Ф(sex + t), p(island + t)). The inflation factor (^c) in the starting model, estimated by performing 1000 bootstrap simulations, indicated that this model showed a negligible overdispersion (^c = 1.19), and the bootstrap goodness of fit test confirmed that the starting model adequately fitted the data (P = 0.164). After model selection for the umbrella model, we included HL. Quadratic term of HL did not improve the final model. We used the Akaike information criterion (Burnham & Anderson 2002) for model selection and ranking of competing models.
Fitness and selection To compare the HFCs of each vital rate, fitness and selection were estimated in a matrix-modelling framework (Lande 1982). The intrinsic population growth rate (k) was estimated from the dominant eigenvalue of the Leslie matrix (Caswell 2001), incorporating all the vital rates at all the age classes of the life cycle. The sensitivity index of k may be interpreted as directional selection (Lande 1982; Caswell 2001). As reproductive success cannot be confidently estimated in males,
Table 2 Summary of the model selection statistics for the umbrella model (see Methods) and models incorporating homozygosity weighted by locus (HL) on survival probabilities. Models were ranked according to their AICc value, with the best supported model in bold Model
AICc
DAIC
(a) umbrella model 264.60 0.000 Ф(sex), p(t) 265.63 1.028 Ф(.), p(t) 266.19 1.590 Ф(sex + t), p(t) 267.54 2.939 Ф(t), p(t) 268.69 4.086 Ф(island), p(t) 270.10 5.50 Ф(sex + t), p(island + t) 279.56 14.95 Ф(sex + t), p(sex + t) 281.73 17.12 Ф(island + t), p(t) 282.52 17.92 Ф(sex + t), p(.) 286.54 21.94 Ф(sex + t), p(island) 287.48 22.88 Ф(sex + t), p(sex) (b) models incorporating genetic terms 259.83 0.000 Ф(sex + HL), p(t) 261.26 1.427 Ф(sex*HL), p(t) 264.60 4.769 Ф(sex), p(t)
AICc Weight
K
0.393 0.235 0.178 0.090 0.051 0.025 0.000 0.000 0.000 0.000 0.000
8 7 10 9 9 12 18 16 5 7 8
0.632 0.310 0.058
9 10 8
intrinsic population growth rates were only estimated in females. A stage (age) matrix model (Fig. S3, Supporting information) was built according to shag population dynamics in the study population (Velando & Freire 2002; Velando et al. 2005a; Martınez-Abraın et al. 2006). Thus, the matrix model included four age classes (juveniles, first-year olds, second-year olds and adults). A balanced chick sex ratio was assumed in the model (Velando et al. 2002). In this model, adult female survival (Sa, Fig. S3, Supporting information) and reproduction (f, Fig. S3, Supporting information) and their variance were dependent on HL, according to our estimation based on genotyped females (see results). The remaining demographic parameters were established in accordance with © 2015 John Wiley & Sons Ltd
H E T E R O Z Y G O S I T Y A N D F I T N E S S I N A S E A B I R D 11 Keller LF, Waller DM (2002) Inbreeding effects in wild populations. Trends in Ecology and Evolution, 17, 230–241. Lande R (1982) A quantitative genetic theory of life-history evolution. Ecology, 63, 607–615. Lebreton JD, Burnhan KP, Clobert J, Anderson DR (1992) Modelling survival and testing biological hypotheses using marked animals: a unified approach with case studies. Ecological Monographs, 62, 67–118. Lieutenant-Gosselin M, Bernatchez L (2006) Local heterozygosity-fitness correlations with global positive effects on fitness in threespine stickleback. Evolution, 60, 1658–1668. Luikart G, Cornuet JM (1998) Empirical evaluation of a test for identifying recently bottlenecked populations from allele frequency data. Conservation Biology, 12, 228–233. Lynch M, Walsh B (1998) Genetics and Analysis of Quantitative Traits. Sinauer Associates Inc, Sunderland, Massachusetts. Martınez-Abraın A, Velando A, Genovart M et al. (2006) Sexspecific mortality of European shags during an oil spill: demographic implications for the recovery of colonies. Marine Ecology Progress Series, 318, 271–276. Miller JM, Coltman DW (2014) Assessment of identity disequilibrium and its relation to empirical heterozygosity fitness correlations: a meta-analysis. Molecular Ecology, 23, 1899– 1909. Molecular Ecology Resources Primer Development Consortium (2010) Permanent genetic resources added to molecular ecology resources database 1(April), pp. 2010–31, May 2010. Molecular Ecology Resources, 10, 1098–1105. Noguera JC, Kim S-Y, Velando A (2012) Pre-fledgling oxidative damage predicts recruitment in a long-lived bird. Biology Letters, 8, 61–63. Olano-Marin J, Mueller JC, Kempenaers B (2011) Heterozygosity and survival in blue tits (Cyanistes caeruleus): contrasting effects of presumably functional and neutral loci. Molecular Ecology, 20, 4028–4041. van Oosterhout C, Hutchinson WF, Wills DP, Shipley P (2004) MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Molecular Ecology Notes, 4, 535–538. Ortego J, Calabuig G, Cordero PJ, Aparicio JM (2007) Egg production and individual genetic diversity in lesser kestrels. Molecular Ecology, 16, 2383–2392. Peakall R, Smouse PE (2006) GENALEX 6: genetic analysis in Excel. Population genetic software for teaching and research. Molecular Ecology Notes, 6, 288–295. Peery MZ, Kirby R, Reid BN et al. (2012) Reliability of genetic bottleneck tests for detecting recent population declines. Molecular Ecology, 21, 3403–3418. Pemberton J (2004) Measuring inbreeding depression in the wild: the old ways are the best. Trends in Ecology & Evolution, 19, 613–615. Raymond M, Rousset F (1995) GENEPOP Version 1.2: population genetics software for exact tests and ecumenicism. Journal of Heredity, 86, 248–249. Rijks JM, Hoffman JI, Kuiken T, Osterhaus ADME, Amos W (2008) Heterozygosity and lungworm burden in harbour seals (Phoca vitulina). Heredity, 100, 587–593. Robinson SP, Simmons LW, Kennington WJ (2013) Estimating relatedness and inbreeding using molecular markers and pedigrees: the effect of demographic history. Molecular Ecology, 22, 5779–5792.
© 2015 John Wiley & Sons Ltd
Rousset F (2008) Genepop’007: a complete reimplementation of the Genepop software for Windows and Linux. Molecular Ecology Resources, 8, 103–106. Samollow PB, Soule ME (1983) A case of stress related heterozygote superiority in nature. Evolution, 37, 646–649. Seddon N, Amos W, Mulder RA, Tobias JA (2004) Male heterozygosity predicts territory size, song structure and reproductive success in a cooperatively breeding bird. Proceedings of the Royal Society B, 271, 1823–1829. Slate J, David P, Dodds KG et al. (2004) Understanding the relationship between the inbreeding coefficient and multilocus heterozygosity: theoretical expectations and empirical data. Heredity, 93, 255–265. Soule ME (1987) Viable Populations for Conservation. Cambridge University Press, Cambridge. Szulkin M, Bierne N, David P (2010) Heterozygosity-fitness correlations: a time for reappraisal. Evolution, 64, 1202–1217. Tiira K, Laurila A, Enberg K et al. (2006) Do dominants have higher heterozygosity? Social status and genetic variation in brown trout, Salmo trutta. Behavioral Ecology & Sociobiology, 59, 657–665. V€ alim€ aki K, Hinten G, Hanski I (2007) Inbreeding and competitive ability in the common shrew (Sorex araneus). Behavioral Ecology & Sociobiology, 61, 997–1005. Velando A (2000) The importance of the hatching date on dominance of young shags. Animal Behaviour, 60, 181–185. Velando A, Alvarez D (2004) Cormoran Mo~ nudo, Phalacrocorax aristotelis aristotelis. In: Libro Rojo de las Aves de Espa~ na (eds Madro~ no A, Gonzalez C, Atienza JC), pp. 60–62. Ministerio de Medio Ambiente-SEO/Birdlife, Madrid. Velando A, Freire J (2002) Population modelling of European shag at their southern limit: conservation implications. Biological Conservation, 107, 59–69. Velando A, Munilla I (2008) Plan de Conservacion del Cormoran Mo~ nudo en el Parque Nacional de las Islas Atlanticas. Universidade de Vigo, Vigo. Velando A, Ortega-Ruano JE, Freire J (1999) Chick mortality in European shag Stictocarbo aristotelis related to food limitations during adverse weather events. Ardea, 87, 51–59. Velando A, Graves J, Freire J (2000) Sex-specific growth in the European Shag Stictocarbo aristotelis, a seabird with size dimorphism. Ardea, 88, 127–136. Velando A, Graves J, Ortega-Ruano JE (2002) Sex ratio in relation to timing of breeding, and laying sequence in a dimorphic seabird. Ibis, 144, 9–16. Velando A, Munilla I, Leyenda PM (2005a) Short-term indirect effects of the ˋPrestigeˊ oil spill on European shags: changes in availability of prey. Marine Ecology Progress Series, 302, 263–27. (2005b) Velando A, Alvarez D, Mouri~ no J, Arcos F, Barros A Population trends and reproductive success of the European shag Phalacrocorax aristotelis on the Iberian Peninsula following the Prestige oil spill. Journal of Ornithology, 146, 116–120. Velando A, Munilla I, L opez-Alonso M, Freire J, Perez C (2010) EROD activity and stable isotopes in seabirds to disentangle marine food web contamination after the Prestige oil spill. Environmental Pollution, 158, 1275–1280. Wang JL, Caballero A, Hill WG (1998) The effect of linkage disequilibrium and deviation from Hardy-Weinberg proportions on the changes in genetic variance with bottlenecking. Heredity, 81, 174–186.
. B A R R O S and P . M O R A N 6 A. VELANDO, A d.f. = 1, P = 0.80) or years (HL*year interaction, Wald v2 = 0.34, d.f. = 1, P = 0.56).
Incorporation of HL in the model explained 2.88% of the residual deviance.
HL and adult survival
HL and fitness
The umbrella model (Table 2a) indicated that recapture probability differed between years (DAIC = 16.33 from a model with time-constant recapture, Fig. S3, Supporting information). Recapture probability was not affected by sex or sampling island (Table 2a). The inclusion of survival differences among sampling islands did not improve the model (Table 2a). The model of survival probability with the lowest AICc included differences between sexes (Table 2a), but this model differed by
