Assessing the role of ontogenetic movement in maintaining population structure in fish using otolith microchemistry

Abstract Identifying the mechanisms maintaining population structure in marine fish species with more than a single dispersing life stage is challenging because of the difficulty in tracking all life stages. Here, a two‐stage otolith microchemistry approach to examining life‐stage movement was adopted, tracking a year‐class from the juvenile to adult stage and inferring larval sources from clustering, in order to consider the mechanisms maintaining population structuring in North Sea cod. Clustering of near‐core chemistry identified four clusters, two of which had either a southern or northern affinity and were similar to juvenile edge chemistry. The other two clusters, common to the central North Sea, had intermediate chemical composition and may have reflected either larval mixing in this region or a lack of geographic heterogeneity in the elemental signature. From the comparison of whole juvenile and the corresponding component of adult otoliths, adults from the southern North Sea mostly recruited from adjacent nursery grounds. In contrast, many adults in the northern North Sea had a juvenile chemistry consistent with the Skagerrak and juveniles from the northern Skagerrak site had a near‐core chemistry consistent with the northern North Sea. Similarities in otolith chemistry were consistent with retention of early life stages at a regional level and also juvenile and adult fidelity. The links between the northern North Sea and Skagerrak indicate natal homing, which when considered in the context of genetic evidence is suggestive of philopatry. The approach used here should be useful in exploring the mechanisms underlying population structuring in other species with multiple dispersive life stages and calcified hard parts.

. Some species may also exhibit a form of philopatry, actively returning to their natal site or region as juveniles or adults (Bonanomi et al., 2016), which may lead to reproductive isolation. Within a generation, individuals may exhibit natal homing to the spawning area they originated from or adopt the migratory behavior of whichever adult group they settle to while originating from a common pool of propagules, that is, the adopted migrant hypothesis (McQuinn, 1997;Petitgas, Secor, McQuinn, Huse, & Lo, 2010). Consequently, an understanding of how population structure is maintained requires knowledge of dispersal among all life stages.
Testing for natal homing is very challenging as it requires the tracking of individuals from fertilization until reproduction (Bradbury & Laurel, 2007). Otolith chemistry has proved useful in examining movements among life stages in regions where there is detectable spatial variation (Campana & Thorrold, 2001).
Unfortunately, as otolith chemistry can vary over time in relation to temperature, salinity, and ontogeny, as well as water mass occupied (Stanley et al., 2015;Walther et al., 2010), microchemistry profiles may not provide a reliable means of spatially tracking an individual. Therefore, connectivity can only be reliably estimated by comparing microchemistry in equivalent parts of the otolith from the same year-class, thus allowing post-dispersed individuals to be assigned to sampled sources (Gillanders, 2002;Munch & Clarke, 2008;Thorisson, Jónsdóttir, Marteinsdottir, & Campana, 2011;Wright, Neat, Gibb, Gibb, & Thordarson, 2006).
Such an approach allows for the estimation of natal homing for those species whose larvae are restricted to estuaries or other confined water masses (Thorrold, Latkoczy, Swart, & Jones, 2001;Walther, Thorrold, & Olney, 2008). However, due to the potential for wide dispersal and high mortality, it is very difficult to track the movements of larvae using otolith chemistry.
This has led some to infer the number of distinct natal sources and the spatial scale over which larval dispersal takes place from a cluster analysis of near-core ablation chemistry (Calò et al., 2016;Gibb, Regnier, Donald, & Wright, 2017).
Atlantic cod, Gadus morhua, is well suited to studying the nature of population structuring because there is already considerable information on their movements (Neuenfeldt et al., 2013) and genetic differentiation (Bradbury et al., 2013), including evidence for natal fidelity (Barth et al., 2017;Bonanomi et al., 2016). In the North Sea, studies of microsatellite DNA have supported a degree of reproductive isolation between the deeper northeastern region between 100 and 200 m and shallower depths, although the neutrality of some key markers has been questioned (Hutchinson, Carvalho, & Rogers, 2001;Nielsen et al., 2009), and single-nucleotide polymorphic markers under selection considerably improve the significance of these differences (Heath et al., 2014;Poulsen, Hemmer-Hansen, Loeschcke, Carvalho, & Nielsen, 2011). At a finer spatial scale, analyses of otolith chemistry indicate that juvenile cod settling off the Scottish coast mainly recruit to their local spawning area . At a similar scale, spawning fidelity has been indicated from tag recapture (Righton, Quayle, Hetherington, & Burt, 2007;Wright, Galley, Gibb, & Neat, 2006) and electronic tagging studies (Neat et al., 2014;Svedäng, André, Jonsson, Elfman, & Limburg, 2010). Physical segregation of cod larvae has been proposed from their distribution relative to frontal areas (Munk, Wright, & Pihl, 2002;Munk et al., 2009) and biophysical model evidence of retention in the northern North Sea (Heath, Kunzlik, Gallego, Holmes, & Wright, 2008), although some larvae are expected to be advected into the Skagerrak (Jonsson, Corell, André, Svedäng, & Moksnes, 2016). Hence, genetic and tagging evidence suggests that the North Sea stock comprises two or more populations: one inhabiting depths >100 m in the northern North Sea and another inhabiting the shallow southern shelf <50 m depth and possibly further coastal groups around the northern UK. Further, both biogeographic barriers and natal homing have been proposed as possible mechanisms shaping population structure in North Sea cod (Heath et al., 2008;Svedäng et al., 2010).
In this study, otolith chemistry was used to examine the role of retention, migrant adoption, and natal homing in shaping population structuring of North Sea cod. Laser ablation inductively coupled plasma mass spectrometry (LA-ICP-MS) was applied to one of the sagittal otoliths from recently settled juveniles to examine the contribution of segregated larval sources (Heath et al., 2008) to nursery grounds using a clustering approach (see Gibb et al., 2017).
To estimate the extent of nursery to spawning ground movement, the other sagittal otolith was used in whole solution (WS-ICP-MS) to define a juvenile baseline for North Sea nursery grounds. Then, adult fish of the same year-class sampled on known spawning grounds during the spawning season were assigned to these nursery grounds based on the chemical signature of the equivalent juvenile region of their otolith, also determined by WS-ICP-MS. The combination of these two approaches carried out on the same individuals gave insights on the roles played by the adopted migrants' hypothesis or natal homing in shaping cod population structure in the North Sea. Finally, otolith evidence was reviewed in the context of reported genetic structuring to consider philopatry.

| Otolith preparation
Juvenile otoliths were decontaminated as per Gibb, Gibb, and Wright (2007). One sagittal otolith from all 431 juveniles collected was stored for WS-ICP-MS. The second sagittal otolith from 15 randomly selected juveniles from each site (except the German coast) was prepared for LA-ICP-MS (n = 135). Samples were randomized with respect to capture site to prevent sample batch bias before being mounted in epoxy resin (Araldite ® CY212, Agar Scientific) in groups of 5-6. Otoliths were sectioned transversely through the core using an Isomet ® 1000 precision linear saw (Buehler, Coventry, UK). The juvenile component of the otolith from spawning adults (age 2 or 3) was prepared for WS-ICP-MS.
Otoliths were cleaned, embedded in epoxy resin, and sectioned as per the juveniles, before microsampling of the juvenile region using a New Wave MicroMill™, following the procedure of Charlier et al. (2006). The otolith was drilled to a depth of 500 μm over two passes (each pass to a depth of 250 μm) to remove most of the juvenile component, which was easily identifiable under reflected light. As a secondary measure to ensure that only material from the 0-group part of the otolith was included, a threshold core length, width and depth limit were predefined electronically, set by averaging the cross-sectional area of otoliths from a select group (n = 4) of 0-group juveniles representing the lower size range (<10 cm TL).  (Heath et al., 2014;Poulsen et al., 2011) are also shown (squares), highlighting the segregation around depths >100 m in the northern North Sea. with a New Wave Research UP213 laser ablation system, using argon gas as the carrier. Pre-ablation runs were undertaken on the standard and each otolith to remove any extraneous impurities.

| Otolith
Material was ablated from two pits (55 μm diameter) on each otolith, one adjacent to, but not covering, the primordium (reflecting natal environment, at an average of 11 increments from the core) and one at the edge (reflecting local environment at capture). At the beginning, after each otolith and at the end of each experimental run of five to six otoliths, an argon gas blank and a glass standard reference material (NIST 612; National Institute of Standards and Technology) were taken for calibration and instrument drift correction. Blank subtracted count data were gathered for each sample ablation in Elan ® 3.4 software and converted to element concentrations (μg/g) in the otoliths by manual calculations using the internal standardization equation described by Longerich, Jackson, and Günther (1996), with concentrations of each element being determined relative to the NIST 612 standard. All isotopes were calculated as ratios relative to 44 Ca, compensating for any variation in ablation yield between samples and standards.

| Juvenile baseline and corresponding component of adult otoliths (WS-ICP-MS)
The preparation of the whole juvenile otoliths and the material milled from the juvenile component of adult otoliths was designed to account for differences in sample size and weight. Juvenile whole otoliths were weighed to the nearest 0.001 mg and then digested in 1 ml 50% ultrahigh purity nitric acid (Seastar). Material from the juvenile component of the adult samples was also digested in a comparable acid, while maintaining the material ratio similar to that of the whole juvenile otoliths. To quantify potential contamination during the cleaning and preparation process, a series of blank acid-washed vials and samples of resin were prepared in the same manner. All samples were randomly ordered with respect to location to prevent sample batch bias although due to the stage-related difference in material analyzed, juvenile otoliths were run as different batches from the 2and 3-group otoliths. All element data were normalized as a ratio of calcium.

| Statistical analyses
As the data did not fully meet multivariate parametric assumptions of normality and homogeneity of variance/covariance, methods that relaxed or did not require these assumptions were used.
Differences in elemental ratios at the ablated otolith edge or in the whole otolith among juvenile sites were tested using Kruskal-Wallis and a Dunn multiple comparison test using the dunn.test The potential for geographic variation in ablated edge chemistry to discriminate juveniles between sampling locations was assessed using random forest (RF) classification (Breiman, 2001), estimated with the randomForest package (v4.6.10) in R 3.3.1. The built-in cross-validation method of the RF approach (Out-of-Bag classification error; Breiman, 2001) was used to assess the accuracy with which sampling locations could be discriminated from each other (Gibb et al., 2017). Variable importance was assessed using the "mean decrease in accuracy", a measure that reflects the loss in classification accuracy when a given variable is left out of the analysis (averaged across all trees; Breiman, 2001).

| Sources of juveniles
Clustering analysis of the near-core ablation data was performed to gain insights into the number of sources of juveniles and the movements of larvae between the locations sampled in this study. A RF clustering method was used in an unsupervised way, hereafter called RF clustering (Régnier et al., 2017). RF clustering is a two-step process. In the first step, unsupervised RF was applied to the standardized dataset. Capture location (relevant for settled individuals) was ignored in this process; instead, a synthetic dataset was created through random sampling from the product of empirical marginal distributions of the different elements. Then, RF was used to separate observed from synthetic data (used as the class factor) and to produce a similarity matrix, defined by the frequency at which two individuals end up in the same terminal node of the trees (Breiman, 2001). In the second step, the dissimilarity matrix (dissimilarity = √(1 − similarity)) was used as an input for partitioning around medoid clustering and the appropriate number of clusters was determined using the Dunn index (see Régnier et al., 2017).

| Adult assignment
The whole juvenile otolith chemistry, derived from WS-ICP-MS, was used as a baseline for assignment of the corresponding juvenile component of adult chemistry. The contribution of each baseline area to the adult samples was estimated from an assignment analysis in a Bayesian framework (or MCMC; Munch & Clarke, 2008). Here, the assumptions of strict multivariate normality and homogeneity of covariance matrices are relaxed by accounting for uncertainty in the baseline.
Individuals were considered to be assigned to one of the 10 sites if the assignment probability was greater than an arbitrary threshold of 0.5.

| Spatial variation in elemental signature of juveniles
The ablated edge chemistry of the juvenile otolith showed signifi- Buchan, close to the 100-m contour, differed from the more southerly Tay station which had an affinity with southern Dogger sites.
There were no misassignments between the northern and southern sites (Table 1). Variable importance measured as the mean decrease in classification accuracy revealed that Mn was the most relevant element (24.5% mean decrease in classification accuracy), followed by Ba (16.1%), Sr (15%), and Mg (7.9%).

| Sources of juveniles
Although juveniles at a site could have originated from multiple sources, there was no significant difference in Ba, Mn, and Mg between the ablated near-core and edge chemistry in eastern Skagerrak juveniles or for Ba and Mn in Buchan juveniles (Wilcoxon's matched-pair test; p > 0.1; Table 2). Ba levels were also consistently high between the near core and edge of northern sites (Shetland and Fisher Bank) and conversely consistently low in the more southern Danish coast juveniles. Near-core Ba levels of north Skagerrak juveniles were also significantly higher than those in the edge (Wilcoxon's matched-pair test; p < 0.001) and similar to the high levels found on the otolith edge of the juveniles of Shetland and Fisher Bank.
The RF clustering of near-core ablation chemistry distinguished between four and seven natal clusters based on early larval elemental values. Given the four regions of similar edge chemistry, the four larval cluster solutions were considered further. The mean and standard deviation of each element value per cluster are given in Table 3.

Bank having high Ba and low Mn, and the southern sites Southern
Bight, German Bight, and Danish coast having high Mn.
Differences in whole-solution juvenile chemistry among sites were reflected in the RF assignment. Classification accuracy by station ranged from 0.25 to 0.87, although individuals were classified to the northern North Sea or southern North Sea to an accuracy of >0.83, as most misassignments were to nearby grounds (Table 4).
Classification accuracy was only 0.59 for the Scottish coastal ground of Buchan due to some similarity in elemental composition with east Skagerrak. Similarly, east Skagerrak also had 0.13 assigned to Buchan, as well as 0.08 and 0.13 assigned to southern and northern sites, respectively. However, classification error was far lower in the north Skagerrak site, with none assigned to Buchan and only 0.06 assigned to both southern and northern sites. Based on the classification error, there is a low likelihood of misassigning southern-tonorthern juvenile sites, but less certainty with assignments to east Skagerrak and Buchan.

| D ISCUSS I ON
Through a combination of clustering of near-core ablation chemistry to consider natal sources and assignment of the juvenile component of adult otoliths to baseline nursery sites, the present study provides a framework for examining how population structure is maintained in species with more than a single dispersive life stage. The approach is dependent on spatial variation in otolith elemental composition (Elsdon et al., 2008) and the assumption that differences in elemental values near the core reflect multiple larval origins (Calò et al., 2016). The elements Mn and Ba gave the most discriminant accuracy and exhibited a latitudinal gradient that provided a natural tag for northern and southern sites.
Smaller-scale differences in Sr and Mg allowed a further degree Sea (Gibb et al., 2007;. Although the observed differences in elemental values cannot be explained, variation in water elemental concentrations has been reported over similar spatial scales (e.g., Balls, Cofino, Schmidt, Topping, & Wilson, 1993), and there are oxygen, salinity, and temperature gradients (Berx & Hughes, 2009;Queste, Fernand, Jickells, & Heywood, 2013) that may account for the observed spatial variation in some of the considered elements (Bath et al., 2000). For example, the regional importance of Atlantic inflow in the north and freshwater from the Baltic and major rivers in the south persists all year around (Munk et al., 2002(Munk et al., , 2009Winther & Johannessen, 2006 (Balls et al., 1993), which is consistent with the higher concentrations of this element in both water (Dehairs, Baeyens, & Van Gansbeke, 1989) and otoliths from the south and east of the study area. Ba concentration tends to increase with depth (IPCS 1990), which is consistent with the highest concentrations of otolith Ba from the northern sites that have a strong Atlantic influence. Importantly, the gradients in otolith elemental values enabled distinction of the northern waters that adults from the Viking genetic unit inhabit (Heath et al., 2014), as well as exploration of the finer-scale structuring suggested by adult movements (Neat et al., 2014;Righton et al., 2007;Wright, Galley, et al., 2006).
As natal clusters were not obtained from sampled larvae, we do not know whether similarities in elemental values between edge and near-core ablation chemistry reflect a similar spatial distribution of larvae and juveniles. Spawning grounds tend to occur further offshore than the areas where cod settle as juveniles (Gibb et al., 2007;Gonzalez-Irusta & Wright, 2016), and because of freshwater input into shallow coastal areas, the chemical landscape of settling cod may be expected to be more varied than that experienced by hatching larvae. This may be one explanation for the narrower range in near-core chemistry compared to edge chemistry.
However, different element signatures between life stages may develop due to ontogenetic variation in element incorporation (De Pontual, Lagardère, Amara, Bohn, & Ogor, 2003). In addition to possible ontogenetic effects, elements affected by temperature should not be directly comparable between the near-core and edge chemistry because they were formed around 11-19°C at the edge and around 7°C near the core, based on the temperature during spawning (Gonzalez-Irusta & Wright, 2016). Stanley et al. (2015) found that otoliths of juvenile cod had around a third higher Mn and Mg in individuals held at 12°C compared to 5°C. Consequently, the lower near-core values of these elements in our study may be partially explained by the difference in temperature. However, the threefold variation in median Mn and Mg values among sites suggests that spatial differences would be far more important than temperature in explaining the observed core and edge differences. Variation in salinity across the study area (33.21-34.5) was also too low to expect that this variable had any significant effect when considered in the context of the study by Stanley et al. (2015).
Clustering of near-core chemistry indicated more than one contributing larval source, consistent with genetic evidence (Heath et al., 2014;Nielsen et al., 2009;Poulsen et al., 2011). The evidence for segregated larval sources between north and south is inconsistent with a common propagule pool of larvae. The scale TA B L E 4 Cross-validated classification probabilities of whole juvenile chemistry for 10 juvenile sites calculated on 10 independent samples of 80 individuals. Different shaded values refer to classification probability for sites in the same region, and the sum classification probability for these regions is also given of similarities between near-core clusters and juvenile edge chemistry, as well as between the juvenile component of adults and whole juvenile chemistry, was largely consistent with the spatial extent of the deep northern and shallow southern population units proposed by Heath et al. (2014) and evident from tagging (Neat et al., 2014;. While this appears suggestive of natal fidelity within the boundaries of proposed population units, inferred movements appeared more complex in the northern Skagerrak, as many juveniles had near-core chemistry that differed significantly from the edge, and many adults in the North Sea had a juvenile Skagerrak signature. The segregation of near-core cluster and edge chemistry between juveniles north and south of the 50-m-deep contour was consistent with predictions from a biophysical model (Heath et al., 2008). In the southern North Sea, spawning is concentrated in coastal waters (Gonzalez-Irusta & Wright, 2016) and egg and larval dispersal tends to be confined by adjacent recurrent salinity fronts (Munk et al., 2002(Munk et al., , 2009. The similarity in juvenile ablated edge and near-core chemistry from southern sites therefore appears consistent with the degree of hydrographic retention in this region. Following settlement, juvenile cod tend to exhibit a degree of site fidelity until 1-year-old (Hawkins, Soofiani, & Smith, 1985;Riley & Parnell, 1984). Tagging studies also indicate that the southern extent of the Viking unit and northern extent of the Dogger unit do not overlap (Neat et al., 2014;Righton et al., 2007;Wright, Galley, et al., 2006). This regional fidelity could explain why most adults in the south were assigned to nearby juvenile sites. Consequently, the Dogger unit in the southern and central North Sea appears to be maintained by a combination of hydrographic retention of larvae and regional fidelity of the settled juvenile and adult cod.
F I G U R E 5 Chart showing assignment of adults to juvenile sites. Pie chart color represents juvenile site assignment, and the size is scaled to the square root of number sampled (maximum n = 46). The depth contours are also shown The similarity between juvenile edge and near-core chemistry seen in Buchan off the Scottish east coast, together with past evidence for a local nursery origin of spawning cod , appears consistent with an even finer spatial scale of natal fidelity than that seen in the southern North Sea. High densities of juvenile cod do occur in the coastal waters of this region (Gibb et al., 2007), which are in close proximity to spawning grounds (Gonzalez-Irusta & Wright, 2016). Unfortunately, too few adults were caught close to the Scottish east coast to verify fidelity during this study.
However, given the similarity in juvenile natal and edge values, it is unlikely that past evidence of the links between adults and nearby juvenile sites  can be explained by the adopted migrant hypothesis (McQuinn, 1997;Petitgas et al., 2010), as such juveniles would be expected to have had a different natal signature than just the region they settled to. Moreover, differences in maturation schedules between this region and the southern North Sea, based on field investigations (Wright et al., 2011) and common environment experiments (Harrald, Wright, & Neat, 2010), are suggestive of local adaptation within a reproductively isolated population.
Patterns of overlap in the chemical signatures between the juvenile component of northern North Sea adults and juveniles from the northern Skagerrak site appear consistent with natal homing. As with juveniles in the northern North Sea, a high proportion of northern Skagerrak juveniles had high Ba near the core, contrary to the corresponding edge chemistry. Advection of larvae into the Skagerrak is quite likely as southwesterly winds dominate during the spawning season of cod, which tends to take eggs and larvae south and east (Winther & Johannessen, 2006) and often into the Skagerrak (Jonsson et al., 2016). The high assignment of northern North Sea spawners to the Skagerrak juvenile signature suggests that the latter area acts as a nursery and could explain why densities of this life stage are so low in the former region. Previous genetic, otolith microchemistry and archival tagging studies in the Skagerrak have demonstrated that there is a mix of resident and nonresident juvenile cod, with a shift in the latter toward the North Sea by around age 2 to 3 (Knutsen et al., 2004;Svedäng, Righton, & Jonsson, 2007;Svedäng et al., 2010). The present study helps explain these earlier findings, as most sampled juveniles in the northern Skagerrak appeared to belong to the deep northern North Sea "Viking" unit. Once recruited to this Viking region, electronic tagging records suggest that cod remain in depths >100 m all year around and experience a much smaller temperature range than the Dogger population unit (Neat et al., 2014).
Juvenile cod from the Northern Skagerrak site have also been found to largely originate from the North Sea (Barth et al., 2017). While future genetic investigation is needed to confirm the links between northern and southern North Sea spawning grounds, this apparent return migration is suggestive of natal philopatry.
Knowledge of life-stage connectivity is clearly important in understanding the changes in population dynamics. In the case of cod, past studies have treated the North Sea as a single unit and correlated the center of distribution with environmental conditions (e.g., Rindorf & Lewy, 2006;Rutterford et al., 2015). However, this study suggests that such a treatment could fail to account for important differences in substock dynamics.
Importantly, this study demonstrates that juveniles move both northward and southward within a year-class, dependent on the nursery area they settle to and their region of origin. Hence, reported northern shifts in North Sea stock distribution are more likely to reflect higher population growth in the Viking than in the Dogger unit. Indeed, accounting for such structuring, Holmes, Millar, Fryer, and Wright (2014) found evidence that the Dogger spawning stock biomass had undergone a much greater decline than that seen in the Viking unit. The links between spawning units and juveniles evident from this study now make it possible to consider the relationships between recruitment and spawning stock at a population level.
The challenge of characterizing population structuring in fully marine species is evident from the paucity of published studies.
Indeed, the question of natal homing in cod has long been debated (Barth et al., 2017;Svedäng et al., 2007Svedäng et al., , 2010 with the only convincing genetic evidence being recently reported (Barth et al., 2017;Bonanomi et al., 2016). While genetic approaches allow for more generalized conclusions about multigeneration philopatry, they can be biased by rare mixing events and the evolutionary timescales over which genetic drift occurs. Genetic differences may also provide little insight into ecologically relevant connectivity (i.e. levels sufficient to influence demographic rates), because only a few individuals need to be exchanged per generation to maintain genetic homogeneity between subpopulations (Waples & Gaggiotti, 2006;Wright, 1931).
Natural chemical tags, such as those used in this study, provide a contemporary indication of movements between life stages and could be used in combination with genetic tools to improve the understanding of the role of natal philopatry in structuring marine fish populations.

ACK N OWLED G EM ENTS
We gratefully acknowledge the assistance of Craig Robinson and Pauline Learmouth in the ICP-MS analysis; staff and crew onboard the NS-IBTS surveys; and Scottish commercial vessels Phoenix and Ocean Bounty for sample collection. This study was funded by Scottish Government projects MF0760 and SP004. Two reviewers provided helpful comments on an earlier version of the manuscript.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
PW conceived the ideas and led the writing of the manuscript; PW and TR designed the methodology and analyzed the data. JA, FG, and SD analyzed samples. All authors contributed critically to the drafts and gave the final approval for publication.
Data are available from Marine Scotland DATA http://data.marine.gov.scot.