Multi-method approach shows stock structure in Loligo forbesii squid

Abstract

Knowledge of stock structure is a priority for effective assessment of commercially-fished cephalopods. Loligo forbesii squid are thought to migrate inshore for breeding and offshore for feeding and long-range movements are implied from past studies showing genetic homogeneity in the entire neritic population. Only offshore populations (Faroe and Rockall Bank) were considered distinct. The present study applied mitchondrial and microsatellite markers (nine loci) to samples from Rockall Bank, north Scotland, North Sea, various shelf locations in Ireland, English Channel, northern Bay of Biscay, north Spain, and Bay of Cadiz. No statistically significant genetic sub-structure was found, although some non-significant trends involving Rockall were seen using microsatellite markers. Differences in L. forbesii statolith shape were apparent at a subset of locations, with most locations showing pairwise differences and statoliths from north Ireland being highly distinct. This suggests that (i) statolith shape is highly sensitive to local conditions and (ii) L. forbesii forms distinguishable groups (based on shape statistics), maintaining these groups over sufficiently long periods for local conditions to affect the shape of the statolith. Overall evidence suggests that L. forbesii forms separable (ecological) groups over short timescales with a semi-isolated breeding group at Rockall whose distinctiveness varies over time.

ecological stocks, genetics, Loligo, microsatellite, Rockall, statolith shape

Introduction

In the recent years, there has been a global expansion in cephalopod fisheries as conventional finfish stocks are depleted (Caddy and Rodhouse, 1998; Hunsicker et al., 2010; Arkhipkin et al., 2021). This trend has been less pronounced within the North East Atlantic, where the bulk of cephalopod (and especially squid) landings are a result of bycatch during demersal fishing activities (Howard et al., 1987; ICES, 2020). Nevertheless, even here, an increasing number of vessels have, from time-to-time, supplemented their standard fishing with seasonally directed squid fishing (Young et al., 2006; Arkhipkin et al., 2015). The veined squid (Loligo forbesii) is among several commercially valued species of cephalopod in the North East Atlantic, which experiences steady levels of fishing pressure due to its respectable market value and status as a non-quota species. Landings are typically reported together with other long-finned squid species (including L. vulgaris and Alloteuthis spp.), which together have landings of ∼8000–12 000 tonnes per year, with French, Portuguese, Spanish, and the UK fleets accounting for > 90% of the catch (Hastie et al., 2009; ICES, 2020). Loligo forbesii is distributed in neritic waters throughout much of the North East Atlantic, to the south as far as west Africa (Senegal, mainly north of 24ºN) and north to the Faroe Islands, and also inhabits much of the Mediterranean (Guerra et al., 2014; Jereb et al., 2015). Most parts of its distribution overlap with the closely related and commercially valued L. vulgaris, but the relative importance of L. forbesii in landings from Spanish and Portuguese Atlantic coasts declined markedly in the 1990s and 2000s (Chen et al., 2006). Although L. forbesii is caught commercially throughout much of its distribution, including directed artisanal fisheries in Iberia; the Moray Firth and Rockall are amongst the few areas that support targeted trawling fisheries (Pierce, et al., 1994a; Young et al., 2006), with landings in other areas primarily being from bycatch fisheries (Young et al., 2006).

Like most squid, L. forbesii completes its short life-cycle within a maximum of 16 months but rarely takes more than 12 months (Guerra and Rocha, 1994). The population(s) and the fishery landings, thus, tend to show a fairly consistent annual cycle—although at longer scales there is some evidence of shifts, possibly related to differing proportions of winter and summer breeders (Pierce et al., 2005). This species also displays a rapid growth rate with younger individuals capable of increasing in size by up to 8% of their own body weight daily (BW/d) under captive conditions, until maturity, by which time growth rate decreases to 2%–4% BW/d (Forsythe and Hanlon, 1989). Individuals typically die after spawning, i.e. they are semelparous (Lum-Kong et al., 1992; Boyle and Ngoile, 1993), however spawning can take place in several batches over a period of time—a life history strategy termed “intermittent terminal spawning” (Rocha et al., 2001). In northern Europe, Loligo forbesii spawning grounds are thought to exist along the coast of Scotland, the North Sea, west of Ireland and into the western Celtic Sea (Collins et al., 1995; Laptikhovsky et al., 2021). Consistent with it being semelparous, Loligo forbesii, like many squid, has non-overlapping generations, such that recruitment of each successive generation is totally dependent on the survival and successful reproduction of the previous one (Boyle, 1990). However, this comes with a temporally spaced-out spawning and recruitment period (indeed, possibly also winter and summer spawning seasons (Pierce et al., 1994a; Boyle et al., 1995)), which acts as a sort of “bet-hedging” strategy to avoid catastrophic loss of spawning stock or recruits in each cycle (Caddy, 1983). At the same time, this species is potentially susceptible to recruitment overfishing whereby individuals are removed from the stock before they can reproduce, to the detriment of future recruitment (Collins et al., 1997). This is particularly true where knowledge is lacking on how the population may be segregated spatially, with the possibility that more isolated components of the metapopulation (i.e. those which receive few recruits from outside) could be over-fished. An example of this is the Rockall population, which is exploited by a directed fishery (ICES, 2020) but is partially isolated from other stocks (Shaw et al., 1999) and is believed to have undergone localised stock depletions in the past (Young et al., 2004).

Knowledge of stock structure is generally poor for commercially fished cephalopods, making this a priority research area (Jereb et al., 2015; Lishchenko et al., 2021), however, stock structure has been assessed previously in L. forbesii using various methods. Allozymes and other molecular methods showed no genetic structure along the European continental shelf locations (including Rockall and Faroe banks) but the Azores samples were sufficiently different to warrant consideration as a subspecies (Norman et al., 1994; Brierley et al., 1995). More sensitive microsatellite DNA markers concluded that the North East Atlantic population is panmictic with the exception of two areas: Rockall and Faroe Banks (Shaw et al., 1999). This was attributed to the highly migratory nature of this species, enabling gene flow between different areas, except where these were segregated by isolating oceanographic regimes and expanses of deep water, as seen in Rockall/Faroe Banks (Shaw et al., 1999). Because populations can share spawning grounds but can be found in discrete units outside of the spawning period, genetic characteristics do not give a full picture regarding the “stock concept” and important information can be gained by including non-genetic markers (Begg et al., 1999). To illustrate the issue with genetic markers: Shaw et al. (2010) performed microsatellite DNA analysis (eight loci) on separate spawning groups in Loligo reynaudii over an area spanning East to West Cape in South Africa, and found little evidence of genetically separate populations. An explanation was provided by Sauer et al. (2000), who undertook a tagging study of Loligo reynaudii in the same area and provided evidence that spawning grounds are shared by individuals with wide-ranging movements, noting that individuals are capable of moving to spawning sites >200 km away within 18 days. Similarly, the population of Doryteuthis opalescens, has two spawning grounds, one in northern California and one in southern California, but despite this, samples across a large area from British Columbia to Mexico were genetically homogenous (using six microsatellite DNA loci), indicating that geographical barriers to gene flow were absent along ≈ 2500 km of coastline (Reichow and Smith, 2001). Thus, even when ecological differences exist (discrete spawning areas, in this example), populations can constitute a single genetic stock.

Under a more “ecological” stock concept, a stock may be considered a management unit or group of individuals which are characterised by similarities that are not necessarily inherited or related to an isolated breeding group, but include those induced by the environment (Begg and Waldman, 1999; Begg and Brown, 2000). Hence group similarities might be based on shared phenotypes or phenological factors (spawning period, for instance). These factors can be linked and are often biologically meaningful, for example: groups of fish can be identified with phenotypic similarities (otolith widths) that are associated with biological characteristics (different juvenile growth rates, indicated in the otolith) caused by phenological factors (autumn- versus winter-spawning adult stock) (Burke et al., 2008). This approach not only helps define a stock, but provides critical information for stock management since less productive stock components (lower growth rates) cannot be sustainably exploited at the same rate as more productive stock units (Begg et al., 1999). This broader stock definition also includes ecological characteristics such as shared feeding grounds, shared nursery grounds or shared migrations. The latter are also relevant for management since stock exchange between particular geographical areas must be known to manage changes in abundance due to migration (Begg et al., 1999; Power et al., 2005).

Under the ecological stock definition, a wider range of tools to investigate stock structure become available and a multi-disciplinary approach combining several tools can be helpful. For example, Van der Vyver et al. (2016) undertook a combined genetic (microsatellite, four loci) and morphometric study of L. reynaudii in similar locations to Shaw et al.’s (2010) study (see above), where they confirmed no significant population genetic structure. However, it was possible to distinguish phenotypic differences among genetically homogenous groups. Morphometric variation (43 variables on the body, beak, gladius, statolith) was observed in samples from either side of Cape Agulhas and south Angola. This was attributed to environmental influences on growth and body shape, which differed between groups (Van der Vyver et al., 2016). Morphometric analysis of body shape (Silva, 2003) or hard structures (Campana and Casselman, 1993; Libungan et al., 2015a) can frequently indicate differences in phenotype. Since squid are soft-bodied organisms, full body morphometric analyses are difficult in this group, though not impossible (see Pierce et al., 1994b; Pierce et al., 1994c; Van der Vyver et al., 2016; Jones et al., 2019). Studies involving hard parts such as beaks, gladii (i.e. internal chitinous shells) and statoliths have provided promising results. For example, changes in isotope signatures (δ¹³C and δ¹⁵N) from the conus of the gladius to its postracum edge revealed ontogenic migrations in Illex argentinus and Doryteuthis gahi (Rosas-Luis et al., 2017) and in Ommastrephes bartramii (Kato et al., 2016). This approach can also be used to map ontogenic migrations according to changes in the trophic niche (Kato et al., 2016; Rosas-Luis et al., 2017), or can be used to understand inshore–offshore migrations when seawater isotopic signatures differ between habitats (Dance et al., 2014). The analysis of migration patterns from statolith elemental signatures can also be combined with statolith shape analysis to discriminate stocks (Green et al., 2015; Jones et al., 2018a).

Statoliths are paired calcified deposits located in the head, and are the cephalopod's primary balance organ (Young, 1989). The morphology of statoliths is unique to each species and these structures contain valuable ecological and life history information, acting as “black box” life recorders (Arkhipkin, 2005). The morphology of statoliths can vary between cephalopods inhabiting different habitat types (e.g. demersal versus pelagic) due to this organ's involvement in acceleration and gravity reception (Arkhipkin and Bizikov, 2000). As with the equivalent structures in fish (otoliths), care must be taken in interpreting morphologies from different groups, since otoliths have been shown to vary between age groups, sexes (Campana and Casselman, 1993), and genotypes (Vignon and Morat, 2010). Analysis of statolith/otolith shape is nowadays carried out using “geometric morphometrics,” a term coined to describe analysis of point co-ordinates in two or three dimensions (Marcus and Corti, 1996; see Adams et al., 2004 for a review). Semi-automated shape analysis methods have recently been developed, including software to decompose an outline (shape) into a series of representative coefficients, which has become widely used due to its accuracy, robustness and speed (Libungan and Pálsson, 2015). Using such semi-automated techniques, statolith shape variations have been observed in Humboldt squid (Dosidicus gigas) sampled from different areas, which were consistent across sampling years (no statistical stock x year interaction), however, the spatial scale examined was rather large, spanning 50 degrees of latitude (Fang et al., 2018). If this technique could also distinguish stock groupings in neritic species and at smaller spatial scales, this would offer a rapid method of distinguishing ecological stocks since the shape analysis itself can be semi-automated. This approach could inform about areas/species where different ecological groupings can be distinguished, despite a degree of mixing at some stage in the lifecycle to make genetic sub-structuring undetectable.

The objective of the current study was to apply a holistic approach to investigate the stock structure in L. forbesii in the North East Atlantic. The genetic structure of this species was investigated using two genetic markers (COI and microsatellites) at ten locations spanning south Iberia to the north of Scotland and Rockall. In tandem, comparisons of statolith shape—an ecological stock marker, and biological data (length distribution, maturity and sex ratio), were conducted at a sub-set of five locations to the north of the range. Biological data were included to provide a starting point for interpreting groupings indicated by other stock markers. The results of this study are intended to indicate the spatial scale at which separate units of stock may be distinguished in L. forbesii and to interpret these in light of previous research.

Methods

Genetic sample collection

Samples from a total of 425 individuals were collected from 10 locations in the North East Atlantic (see Table 1 and Figure 1). All samples were taken aboard research surveys, with the exception of one sampling location (English Channel) where samples were obtained via port sampling. To obtain reasonable sample sizes (Hale et al., 2012), individuals sampled from multiple research surveys and sampling years/seasonal quarters were grouped to comprise each sampling “location,” assuming that these came from stations that were judged to be sufficiently close together. A full breakdown of spatially referenced sampling stations and survey information is provided in Supplementary Table S1. Tissue sampling took place after sorting of the catch by species on-board each research survey, or in the case of port samples, after returning to the laboratory. Mantle tissue samples were taken after first removing the skin, and tissue was transferred to 2 mL Eppendorf tubes containing 96% ethanol. All tissue samples were stored at −20°C.

Figure 1.

(Left) Sampling stations from research surveys and port sampling at Port en Bessin (see Supplementary Table S1) for L. forbesii microsatellite and COI genetic analysis; (Right) Sampling stations for statolith shape analysis overlain on bathymetry, with sampling year and sampling quarter information included. Red callout box and ellipses display the five sampling locations which spatially overlapped in both genetic and statolith shape analyses (i.e. north Scotland, Rockall, north Ireland, west Ireland, and south Ireland), with sample sizes included.

Table 1.

Sample location details and basic genetic diversity indices including the average no. of individuals genotyped at each locus (N_g), the no. of alleles (N_a), observed (H_o) and expected (H_e) heterozygosities, the fixation index (F), the inbreeding co-efficient (F_IS), and allelic richness (A_r).

Sample location	Code	Year and Quarter	Sample size	Ng	Na	Ho	He	F	FIS	Ar
North Scotland (NS)	1	19 Q4	27	25.9	13.9	0.790	0.868	0.088	0.109	9.5
Rockall (RK)	2	19 Q3	36	35.6	14.9	0.814	0.860	0.050	0.066	9.3
North Ireland (NI)	3	18 Q4/19 Q1,Q4	149	143.4	19.9	0.796	0.878	0.094	0.098	9.7
West Ireland (WI)	4	18 Q4/19 Q1	53	52.4	16.2	0.803	0.866	0.074	0.083	9.5
South Ireland (SI)	5	18 Q1/18 Q4	45	44.8	17.6	0.815	0.872	0.064	0.076	9.9
North Sea (NSea)	6	19 Q1	13	12.8	11.2	0.778	0.859	0.093	0.134	10
English Channel (ECh)	7	19 Q3	29	28.7	13.9	0.809	0.868	0.069	0.086	9.5
North Spain (NSp)	8	19 Q3	25	24.9	13.3	0.804	0.860	0.061	0.086	9.4
Cadiz (Cad)	9	19 Q4	25	23.7	13.1	0.771	0.852	0.091	0.117	9.3
Bay of Biscay (BB)	10	19 Q4	23	20.4	13.3	0.838	0.873	0.041	0.068	10
Mean					14.7	0.802	0.866

DNA extraction and COI Analysis

DNA was extracted using the Invitrogen™ PureLink™ Genomic DNA Mini Kit, following the Mammalian Tissue and Mouse/Rat Tail Lysate protocol. DNA was eluted in 100 μL Genomic Elution Buffer and stored at –20°C.

DNA barcoding, which targets the 650 bp fragment of the COI gene, was performed using forward LCO1490: 5′-GGTCAACAAATCATAAAGATATTGG-3′ and reverse primers HC02198: 5′-TAAACTTCAGGGTGACCAAAAAATCA-3′ of Folmer et al. (1994). Primers were purchased from Biolegio (Nijmegen, Netherlands) and diluted using dH₂O to a concentration of 100 μM. Each PCR contained 12.5 μL of Thermo Scientific™ DreamTaq Green PCR Master Mix (2X), 0.5 μM of each primer, 2.5 μL of DNA and 9 μL H₂O, resulting in final reaction volume of 25 μL. A negative control to ensure cross-contamination did not occur was incorporated to the procedure by adding 2.5 μL of dH₂O. The PCR conditions were 94°C for 2 min, 35 cycles of 94°C for 40s,  50°C for 40s, and 72°C for 90s,  followed by 72°C for 10 min (Allcock et al., 2007). PCR products were visualised on a 1% (w/v) agarose gel using SYBR safe DNA Gel stain Invitrogen™ and bands of 650 bp were obtained. PCR products were cleaned using Invitrogen™ PureLink™ PCR Purification Kit according to the manufacturer's instructions. Purified PCR products were then standardised to 12 ng/μL using a Biochrom SimpliNano NanoDrop Spectrophotometer in accordance with the DNA sequencing facility specifications. Samples were prepared for sequencing by adding 5 μL of each purified PCR product to 5 μM forward primer LCO1490 resulting in a 10 μL reaction volume. Samples were sent to Eurofins Genomics (Germany) for DNA Sequencing on an ABI 3730XL DNA Analyser. All sequences were viewed and trimmed using Sequence Scanner v2 before input into MEGA X software (Molecular Evolutionary Genetics Analysis across computing platforms, Kumar et al., 2016) where they were aligned using MUltiple Sequence Comparison by Log- Expectation (MUSCLE) and inspected for mitochondrial diversity. Using 558 bp of COI sequence, a median-joining network was constructed using the median algorithm of Bandelt et al. (1999), implemented in the programme POPART v.1.7.1 (Leigh and Byrant, 2015). DnaSP v.10.01 (Librado and Rozas, 2009) was used to estimate nucleotide diversity across the COI dataset.

Microsatellite analysis

A total of nine microsatellite loci was used to assess the nuclear genetic diversity: Lfor 2, Lfor 3, Lfor 4, Lfor 5, Lfor 8 (Shaw, 1997), Lfor 1 (Shaw et al., 1999), and Lfor 11, Lfor 12, and Lfor 16 (Emery et al., 1999). The reverse primer of each set was “pigtailed” using the following bases “5”-GTTTCTT’ to enhance non-templated nucleotide addition according to Brownstein et al. (1996). All primers were purchased in a lyophilised state from Applied Biosystems (Life Technologies) and eluted to 100 μM. The forward primer of each set also had a fluorescent label attached (see below).

The following primer sets were amplified together in three multiplexes: Mix A: Lfor 5 (0.25 μM, 6-FAM™), Lfor 1 (0.5 μM, PET™), Lfor 12 (0.5 μM, 6-FAM™), Lfor 16 (0.5 μM NED™); Mix B: Lfor 4 (0.5 μM, VIC®), Lfor 8 (0.5 μM, 6-FAM™) and Mix C: Lfor 2 (0.33 μM, 6-FAM™), Lfor 3 (0.33 μM, VIC®) Lfor 11 (1 μM,6-FAM™). Each multiplex contained 5 μL Thermoscientific DreamTaq Hot Start Green PCR Master Mix (2X), the relevant quantities of each forward and reverse primer as listed above, 1 μL of DNA and H₂0, resulting in a final reaction volume of 10 μL. Samples were analysed in duplicate from the PCR stage to ensure repeatability. A touchdown PCR was used, with the annealing temperature dropping by 0.5°C for each cycle. This consisted of 94°C for 3 min, 16 cycles of 94°C for 30 s, 61°C for 30 s and 72°C for 30 s followed by 24 cycles of 94°C for 30 s, 53°C for 30 s and 72°C for 30 s and a final extension of 72°C for 10 min. The PCR products were diluted (1:4) using nuclease-free water prior to fragment analysis, and 1 μL was then added to 15 μL of Hi Di Formamide with 0.2 μL of GeneScan™ 500 LIZ™ dye size standard.

Fragment analysis was performed on an ABI Genetic Analyser 3500 and alleles were scored using GENEMAPPER Software v5 (Applied Biosystems). Any samples which gave inconsistent scores between replicates or failed entirely were repeated from the PCR step.

Statolith shape & biological data sampling

All statolith samples came from the same research surveys as genetics tissue samples, apart from the IAMS 2020, when only statoliths were sampled (no genetics samples were taken from that cruise since the same area had already been sampled–Figure 1; Supplementary Table S1). However, the same individual squid were not necessarily sampled for both statoliths and genetics at each sampling station. Statoliths were collected during demersal fishery surveys (2018–2020) at a subset of the locations that were sampled for genetics analysis (Figure 1, 9). RKS19, WSIBTS19, and IGFS18 are all components of the International Bottom Trawl Surveys (IBTS) and used the standard Grande Overture Verticale (GOV 36/47) net with mesh sizes ranging from 100 mm at the trawl opening to a 20 mm liner in a 25 mm cod end (Stokes et al., 2015). IAMS20 used a commercially derived otter trawl with mesh sizes of 200 mm at the trawl wings reducing to 100 mm in the cod-end (Reid et al., 2007). Fishing times during RKS19, WSIBTS19 and IGFS18 were standardised to roughly 30 min at 4 knots and only occurred during daylight hours. Fishing times during IAMS20 were roughly 1 h at 4 knots and occurred over 24 h.

Catches were sorted on-board and the L. forbesii portion of the catch was put aside for biological sampling (mantle length, weight, sex, and maturity). When feasible, biological sampling was carried out for all individuals, but for stations with abundant squid, the subsampling procedure outlined in ICES (2017) was used, i.e. a representative subsample comprising >10 times the number of length classes in the sample (based on 1 cm increments in mantle length) (Gerritsen and McGrath, 2007). Maturity was determined on-board research vessels using the WKMSCEPH common Teuthida maturity scale (Table 2; ICES, 2010) and individuals reaching minimum maturity stage 2a were frozen (−20°C) for subsequent statolith analysis. Selecting a minimum maturity level ensured reliable identification of species and avoided potential artefacts caused by ontogenic effects on statolith shape and size (Fang et al., 2018). A total of 417 individuals (males and females) was selected for statolith analysis across five sampling areas (Rockall, n = 93, north Scotland, n = 96, north Ireland, n = 125, west Ireland, n = 56, south Ireland, n = 47). All statoliths were extracted using the methods of Clarke (1978) using fine-nosed forceps. Extreme care was taken at this stage as the wing portion of the statolith is exceptionally fragile and prone to chipping. Both left and right statoliths were identified and stored (dry) separately. In the majority of cases, only the right statolith was selected for subsequent analysis, however, if the right statolith was damaged or otherwise unusable, the mirrored image of the left statolith, obtained using editing software, was used. Mirror images were justified on the basis that studies which compared left and right statoliths have found that these did not differ in shape (Barcellos and Gasalla, 2015; Green et al., 2015; Fang et al., 2018; Guo et al., 2021). If both statoliths displayed damage, the sample was discarded. From the 417 sampled individuals, 10 were rejected from the analysis due to visible damage on both left and right statolith. The remaining 407 viable statoliths (Rockall, n = 90; north Scotland, n = 96; north Ireland, n = 122; west Ireland, n = 56; south Ireland, n = 43) were imaged with a stereomicroscope (Olympus SZX16, 400x) using reflected-oblique light against a black background. Exposure time and white balance were calibrated in order to achieve the strongest contrast between the statolith and background. The statolith from the right side of each individual was positioned with the lateral dome facing left and anterior rostral lobe facing right, with the wing facing upwards (Figure 2). A 1 mm scale was included for each magnification, to create a calibration measurement of the number of pixels along the 1 mm scale that was later used for calculating the area, length, perimeter, and width of each statolith, obtained using “ImageJ” (v1.52A) software (Schneider et al., 2012).

Figure 2.

Orientation of statolith for imaging showing labelled sub-areas of this structure. The statolith from the right side of each individual was positioned with the lateral dome facing left and anterior rostral lobe facing right with the wing facing upwards.

Table 2.

ICES WKMSCEPH maturity scale (ICES, 2010) for both males and females used in data collection on scientific surveys.

Male	Female	Maturity Stage
No visible maturation	No visible maturation	0
Small translucent testis. Spermatophoric complex visible	Nidamental glands translucent	1a
	Ovary is semi- transparent
Testis enlarged. Spermatophoric complex more developed	Nidamental glands opaque, Ovary enlarged	2a
	Some Oocytes present
Testis at full size. Needham's sack with partially developed spermatophores	Nidamental gland large	2b
	Ovary with some non-hydrated oocytes present
Testis at full size. Needham's sack full of spermatophores	Nidamental gland large	3a
	Ovary packed with large, hydrated oocytes
Very few/no spermatophores in Needham's sack	Very few/no eggs in oviduct Nidamental glands degenerating	3b

Data analysis

Microsatellite statistical analysis

Micro-Checker v2.2.3 (Van Oosterhout et al., 2004) was used to detect the presence of null alleles, stuttering errors or allelic dropout. The full dataset was analysed initially, followed by subdivision of the dataset into the 10 sampling locations and repetition of the analysis. Genetic diversity was assessed using GenAlEx v6.503 (Peakall and Smouse, 2006, 2012) software. Number of alleles (N_a), Observed heterozygosity (H_o), Expected heterozygosity (H_e) and Fixation index (F) were calculated on both full and subdivided datasets and displayed in Table 1. Allelic frequency and principal co-ordinates analysis were also performed using GenAlEx v6.503 software. Genotypic linkage disequilibrium was tested for all locus pairs in each location and for each locus pair across all locations using Fisher's Method (1000 demorisations and 5000 iterations) in GENEPOP v4.2. The inbreeding co-efficient F_IS was calculated in FSTAT v2.9.4 (Goudet, 1995) by randomising the alleles between all individuals in a sample and comparing this to the observed data to determine whether there were any deviations from Hardy–Weinberg Equilibrium, using 10 000 permutations. Population structure was assessed using STRUCTURE v2.3.4 software (Pritchard et al., 2000), which uses Bayesian clustering methods on allele frequencies to find the most likely number of populations (K) present within the dataset. This was performed on both the full and subdivided (by sampling location) dataset. The range of K was set between 1 and 5, with each model consisting of 5 runs. Each run contained a burn-in period of 200 000 with 500 000 MCMC reps. An admixture model was used with the correlated allele frequencies parameter selected and the data were run both with and without the LOCPRIOR function which is used to infer prior location information of the samples. The programme STRUCTURE Harvester (Earl and Von Holdt, 2012) was used to implement the Evanno method (Evanno et al., 2005), in which the L(K) and Δ(K) values were used to identify the most likely K value. Evanno et al. (2005) showed that the use of Δ(K) does not allow the assessment of K = 1 as a potential solution so that assessment of the L(K) graph is the most appropriate analysis when K = 1. POPHELPER R package (Francis, 2017) was used to graphically illustrate the results obtained using STRUCTURE software. To further examine the presence of genetic structure, a principal coordinates analysis (PCoA) in GenAIex v.6.5b (Peakall and Smouse, 2006) was used as a multivariate approach to supplement the STRUCTURE analysis.

Genetic differentiation was evaluated with the software FSTAT v2.9.4 (Goudet, 1995) using randomisation of genotypes. A pairwise test of F_ST values was generated under 1000 permutations and Bonferroni correction for multiple tests was applied. A Mantel test was performed in GenAlEx v6.503 software (Peakall and Smouse, 2006, 2012) to test for isolation by distance in all samples for which geographic distance was available (i.e. excluding the English Channel location, as these samples were taken by means of port sampling and accurate GPS co-ordinates of the catch location were not available). Geographical distances were calculated based on the distances between sampling stations. GPS co-ordinates from a total of 57 sampling stations were entered into GenAlEx to produce a tri matrix table of logged geographical distances between each station. Linearised F_ST values were then plotted against the logged geographical distance.

Statolith shape and biological data analysis

Shape analysis was completed using the ShapeR package (Libungan and Pálsson, 2015) in R (R Core Team, 2019). The first stage of this process was to extract the outlines for all statoliths and write them into the ShapeR package. This operation converts images into grey-scale and picks out the outline using a pre-set threshold pixel value. From this, shape coefficients for each statolith are produced using the discrete wavelet-transform method, which decomposes images into a number of “wavelets” that have an amplitude that begins at zero, increases, and then decreases back to zero and is irregular in shape and compactly supported (Graps, 1995). Due to their irregular shape, wavelets are ideal for analysing complex shapes with many discontinuities or sharp edges, such as statoliths (Fang et al., 2018) and are better suited to this task than Fourier-transform shape analysis which does not approximate sharp edges as accurately (Graps, 1995; Libungan et al., 2015b). The ShapeR package uses this technique to produce a total of 64 wavelet coefficients (five wavelet levels) which describe the shape for each statolith (Libungan and Pálsson, 2015; Libungan et al., 2015b). To take into account any allometric relationships with mantle length, wavelet coefficients which show interaction (p < 0.05) with length, are omitted automatically and the remaining coefficients are standardised for length. Variation in statolith shape was then visually represented by plotting mean statolith shape for all five groups (i.e. sampling locations) using length-standardised wavelet coefficients. The matrix of wavelet standardised coefficients was analysed amongst groups using canonical analysis of principal coordinates (CAP) based on Euclidean dissimilarity indices (Libungan, 2015). Each individual statolith was ordinated using CAP, allowing distances between groups of statoliths to be examined. A cluster plot was used to visually represent the grouping of CAP results in two dimensions. Following this, the partitioning of variation among groups based on their canonical score was tested with an ANOVA-like permutation test (permANOVA) using the Vegan package for R (Oksanen et al., 2013). The permANOVA was then repeated for each pair of sampling locations in order to investigate pairwise associations in statolith shape. To protect against the influence of type one errors, the p-value accepted for significance in this pairwise analysis was adjusted using a Bonferroni correction.

Data from biological sampling (dorsal mantle length (DML), sex and maturity) were also compared between locations for the same individuals which had been analysed using statoliths (n = 407). Sex and maturity ratios were compared between all locations using a chi-square test. Prior to analysis, maturity information was converted to a binary (0, 1) scale with individuals of maturity 2a being designated 0 (“maturing”) and those of 2b and 3a being designated 1 (“mature”). DML data were analysed separately for each sex, comparing between locations (male: n = 191, female: n = 216) after first being tested for normality and homoscedasticity using Kolmogorov–Smirnov and Levene's tests. Where transformations failed to normalise the data, non-parametric tests (Kruskal–Wallis) were used to analyse between-location variation in DML of both males and females, followed by pairwise comparisons using a Dunn–Bonferroni post-hoc test.

To help interpret L. forbesii statolith shape comparisons between locations, a small sample of Loligo vulgaris was also analysed to provide context, i.e. to place the shape variability across L. forbesii samples into perspective when compared with a closely related species. For this, a sample of L. vulgaris (n = 20) was obtained during the IGFS18 in the Celtic Sea (ICES Division 7g) from 7th to 10th December 2018. This was analysed following the same statolith processing and statistical analysis steps as before and was compared with L. forbesii data across all locations.

Results

Mitochondrial DNA (COI) Analysis

A total of eight haplotypes were found in this study. The most common haplotype, H1, occurred throughout the sampled area, with the remaining seven haplotypes occurring in relatively few samples (Supplementary Figure S1). Across the 558 bp of sequence, 547 sites were monomorphic, 11 polymorphic (including a total of 11 mutations), seven sites included singletons and four sites were parsimony-informative. Additional metrics included an overall nucleotide diversity, Pi = 0.00621 ± 0.00142 and gene or haplotype diversity, Hd = 1.000 ± 0.063.

Microsatellite analysis

Most samples were successfully genotyped in duplicate with success rates ranging from 88% in Lfor 8 to 100% in both Lfor 1 and Lfor 3. There was no evidence that stuttering or allele dropout resulted in scoring errors. There was, however, some evidence of null alleles (based on the combined probability for homozygote frequencies at all allele size classes, p < 0.001) in the following three loci: Lfor 8, Lfor 12, and Lfor 16. Analysis of the subdivided dataset by sampling location, showed that none of the loci consistently produced significant levels of null alleles across all locations. Statistical analysis was completed with and without the worst affected loci for null alleles, Lfor 8 and Lfor 12, which were found to have a negligible impact on the overall dataset. Hence the results for the full dataset including all nine loci are presented below, unless otherwise indicated.

Genetic diversity indices are displayed in Table 1 for nine microsatellite loci at 10 sampling locations (n = 425 individuals). The sampling locations exhibited high mean values for the N_a index (numbers of alleles) ranging from 11.2 (North Sea) to 19.9 (north Ireland), averaging 14.7 overall. Allelic richness (A_r) was also high, ranging from 9.3 to 10. Observed and expected heterozygosities per sampling location were very similar overall, respectively ranging from 0.771–0.838 (average 0.802) to 0.852–0.878 (average 0.866), indicating high diversity within all locations.

None of the inbreeding co-efficient (F_IS) values obtained from the 10 sampling locations deviated from Hardy–Weinberg Equilibrium at p < 0.0006 following Bonferroni correction for multiple comparisons (Table 1). Significant linkage disequilibrium was detected at only one pair of loci (north Ireland Lfor 5 x Lfor 2) at p < 0.0001 following Bonferroni correction. No significant linkage disequilibrium was detected for any locus pair across all populations at p < 0.0001 after Bonferroni correction.

Genetic structure

Pairwise F_STvalues between sampling locations were very low, both with (−0.009 to 0.005) and without (−0.007 to 0.008) the worst affected loci for null alleles (Lfor 8 and Lfor 12) (Table 3). This indicates that the presence of null alleles in loci Lfor 8 and Lfor 12 did not impact on the overall result that genetic structure is lacking between sampling locations. When all nine loci were included, the pairwise F_ST resulted in one significant pairwise test (Rockall versus north Ireland) at p < 0.001 after Bonferroni correction for multiple tests was applied (data not shown), but when Lfor 8 and Lfor 12 were removed from the data, none of the sampling locations were significantly differentiated from one another once Bonferroni correction had been applied (Table 3). Despite not meeting the threshold for statistical significance, seven comparisons out of 45 were “significant” at p < 0.05 prior to Bonferroni correction, six of which involved Rockall (Table 3). These six locations were: north Scotland, north Ireland, west Ireland, English Channel, north Spain and Cadiz, with pairwise F_ST values versus Rockall ranging from 0.003 to 0.008 (Table 3). The remaining comparison showing p < 0.05 prior to Bonferroni correction was west Ireland versus Cadiz (pairwise F_ST of 0.004–Table 3).

Table 3.

Results of pairwise F_ST values between sampling locations using seven microsatellite loci (i.e. omitting loci Lfor 8 and Lfor 12, which were affected by null alleles). Underlined values indicate statistical significance at p < 0.05 which became non-significant after Bonferroni correction was applied (cut-off for statistical significance at p < 0.001). Sample location abbreviations are listed in full in Table 1.

Sample location	NS	RK	NI	WI	SI	NSea	ECh	NSp	Cad	BB
NS	–
RK	0.003	–
NI	−0.001	0.005	–
WI	−0.004	0.004	0.001	–
SI	−0.001	0.003	−0.001	-0.002	–
NSea	−0.005	−0.002	−0.004	0.000	−0.006	–
ECh	−0.003	0.007	0.001	0.000	−0.001	−0.002	–
NSp	−0.002	0.003	−0.001	0.000	−0.003	−0.005	-0.003	–
Cad	−0.002	0.008	0.000	0.004	0.004	0.002	−0.003	−0.006	–
BB	−0.004	0.003	−0.004	−0.003	−0.004	−0.003	−0.005	−0.005	−0.007	–

The STRUCTURE results showed that there was no evidence of genetic partitioning in the dataset and the population displays high levels of admixture—see bar plots for K = 2, K = 3, K = 4, and K = 5 (Figure 3a). STRUCTURE analysis indicated a population of K = 1, supported by the L(K) graph, which showed the most likely value to be K = 1, based on posterior probabilities (Figure 3b). The PCoA (Supplementary Figure S2) was also used to assess whether any genetic patterns were present and this showed no evidence of clustering by sampling locations. Mantel's test revealed significant Isolation By Distance (IBD), with a weak negative trend and a very low fit (R² = 0.0059, p < 0.05, Supplementary Figure S3).

Figure 3.

(a) STRUCTURE assignment of individuals across all subdivided sample sites using LOCPRIOR function for the following values of K (K = 2), (K = 3), (K = 4), and (K = 5) displayed via the POPHELPER package. The numbers on the x-axis refer to the sample locations listed in Table 1. (b) Mean likelihood L(K) ± standard deviation per K value produced by STRUCTURE Harvester and visualised in POPHELPER indicating that 1 is the true value of K.

Biological data

No significant between-location variation in sex ratio was found (χ² = 7.220, 4 df, n = 407, p > 0.05). However, maturity ratios (maturing/mature) varied significantly between locations with Rockall containing a higher frequency of “maturing” individuals (stage 2a) and all other locations containing more “mature” individuals (stages 2b and 3a) (χ² = 72.618, 4 df, n = 407, p < 0.05, Figure 4). The K-S test statistic indicated that neither male (K-S D = 0.133, p < 0.05) nor female (K-S D = 0.115, p < 0.05) DML data were normally distributed and a Levene's test revealed that the assumption of homogeneity of variance was met only in males (F = 1.364, n = 191, 4 df, p > 0.05) but not in females (F = 4.994, n = 216, 4 df, p < 0.05), despite attempts at transformation, hence a Kruskal–Wallis test was used on untransformed DML data of both sexes. This showed significant differences in DML between sampling locations in males (H = 9.231, 4 df, n = 191, p < 0.05, Figure 4) but not in females (H = 36.669, 4 df, n = 216, p > 0.05). Pairwise comparisons using a Dunn–Bonferroni post-hoc test conducted on male DML data (Table 4 with significance cut-off adjusted for multiple comparisons) showed significant differences in DML between south Ireland (median DML 170 mm) and north Ireland (median DML 220 mm), south Ireland and west Ireland (median DML 275 mm), south Ireland and Rockall (median DML 290 mm), and between north Scotland (median DML 190 mm) and north Ireland, west Ireland, and Rockall. Thus, the largest median size was seen in males from Rockall, followed by west Ireland, north Ireland, north Scotland, and south Ireland, however, the data distributions were highly overlapping and multi-modal (Figure 4).

Figure 4.

Histograms of dorsal mantle length in male Loligo forbesii. The magnitude of the x- and y-axis is the same in each location to aid comparison. Arrows indicate median length for each location.

Table 4.

Dunn–Bonferroni pairwise comparisons between locations based on Kruskal–Wallis comparison of male DML in Loligo forbesii. SI = south Ireland, NS = north Scotland, NI = north Ireland, WI = west Ireland, RK = Rockall. Values in bold/underline are statistically significant.

Comparison	Median DML (mm)	Range (mm)	Test Statistic	Std. Test Statistic	Adj. Sig.
S Ireland = N Scotland	SI: 170	SI: 410	15.339	1.083	1.000
	NS: 190	NS: 260
S Ireland < N Ireland	SI: 170	SI: 410	−51.376	−3.810	0.001
	NI: 220	NI: 300
S Ireland < W Ireland	SI: 170	SI: 410	−64.708	−3.599	0.003
	WI: 275	WI: 350
S Ireland < Rockall	SI: 170	SI: 410	65.906	4.655	0.000
	RK: 290	RK: 220
N Scotland < N Ireland	NS: 190	NS: 260	−36.037	−3.332	0.009
	NI: 220	NI: 300
N Scotland < W Ireland	NS: 190	NS: 260	−49.369	−3.071	0.021
	WI: 275	WI: 350
N Scotland < Rockall	NS: 190	NS: 260	−50.567	−4.343	0.000
	RK: 290	RK: 220
N Ireland = W Ireland	NI: 220	NI: 300	13.332	0.861	1.000
	WI: 275	WI: 350
N Ireland = Rockall	NI: 220	NI: 300	14.530	1.343	1.000
	RK: 290	RK: 220
W Ireland = Rockall	WI: 275	WI: 350	1.198	0.075	1.000
	RK: 290	RK: 220

Statolith shape, CAP ordination, and mean plots

The statolith outline was reconstructed with > 98.5% accuracy using five wavelet levels (Supplementary Figure S4a). Visual representation of mean statolith shape under discrete wavelet reconstruction showed observable differences between all locations along the wing, rostral angle, lateral dome and dorsal dome regions of the statolith (Figure 5). An intraclass correlation (ICC) plot displaying the proportion of variance among groups (i.e. locations) in wavelet coefficients confirmed that most of the variation among groups can be attributed to angles 140–180°, 190–230° (lateral/dorsal dome), and 310–360° (wing/rostrum) (Supplementary Figure S4b). North Ireland was strongly distinguished from other locations in overall mean statolith shape, especially in the dorsal and lateral dome and wing section (Figure 5). Highest perimeter-area ratios, which correspond to shape complexity in the statolith (measured in two dimensions) were observed in south Ireland, whereas these ratios were lower in north Scotland, north Ireland, and Rockall (Supplementary Table S2). Figure 6displays canonical principle coordinates (CAP) scores based on Euclidean dissimilarity indices, including mean scores per location ± standard error. The first discriminating axis (CAP1) explained 80.5% of variation among locations while the second axis (CAP2) explained 15.8% of variation. Canonical values for both Rockall and north Scotland oriented to the upper right and centre of the plot while south and west Ireland values were situated in the centre and lower right of the plot. North Ireland, on the centre left, was found to be most distant from the other locations. Squid size did not affect this ordination since there was no correlation between DML and CAP1 score (Supplementary Figure S5). With the exception of Rockall versus north Scotland, comparisons of scores between pairs of locations were all significantly different from one another by permANOVA (F = 27.423, 4 df, adjusted p-value cut-off = 0.005, n = 407, Table 5).

Figure 5.

Mean statolith shape of L. forbesii (males and females) in north Ireland (black), north Scotland (red), Rockall (green), south Ireland (blue) and west Ireland (cyan) under discrete wavelet reconstruction. The numbers show angle in degrees (°) based on polar coordinates where the centroid of the statolith is the center point of the polar coordinates.

Figure 6.

Canonical scores on discriminating axes CAP1 and CAP2 for each L. forbesii location: north Ireland (I, black), north Scotland (N, red), Rockall (R, green), south Ireland (S, blue) and west Ireland (W, cyan). Large black letters represent the mean canonical value (± standard error) for each location and smaller coloured letters represent individual squid (SI = south Ireland, NS = north Scotland, NI = north Ireland, WI = west Ireland, RK = Rockall).

Table 5.

PermANOVA results of CAP ordination of wavelet coefficients between all 10 pairs of sampling locations. Significance between each pair of areas was analysed using a Bonferroni–adjusted p-value cut-off of 0.005. All comparisons were statistically significant except where indicated (NS).

Comparison	df	F	P
S Ireland—N Scotland	1	7.34	< 0.001
S Ireland—N Ireland	1	29.04	< 0.001
S Ireland—W Ireland	1	2.94	0.004
S Ireland—Rockall	1	9.82	< 0.001
N Scotland—N Ireland	1	57.79	< 0.001
N Scotland—W Ireland	1	11.37	< 0.001
N Scotland—Rockall	1	2.11	0.040 (NS)
N Ireland—W Ireland	1	40.58	< 0.001
N Ireland—Rockall	1	60.97	< 0.001
W Ireland—Rockall	1	14.08	< 0.001

The shapes of L. forbesii and L. vulgaris statoliths were also analysed using CAP ordination scores. This time, CAP1 explained 14.9% of the disparity between species/locations while CAP2 explained 76.6% of this variation. Differentiation between L. vulgaris and L. forbesii shows that canonical scores for L. forbesii from all locations oriented much the same as in Figure 6 and L.vulgaris scores clustered between the L. forbesii scores at Rockall/north Scotland and south/west Ireland but distant from the L. forbesii from north Ireland (Supplementary Figure S6). A permANOVA showed that the shape coefficients in the two species were statistically different from one another (F = 23.408, 5 df, p < 0.05).

Discussion

The genetic markers examined in the present study showed no significant evidence of sub-structure in Loligo forbesii from 10 locations across the NE Atlantic but some non-significant evidence of structure was observed at Rockall. By contrast, good spatial differentiation was seen in the statolith shape markers, highlighting the fact that ecologically separated groups can be identified in this species. New insight from both genetic and ecological markers allows us to re-interpret the stock structure of L. forbesii, both in locations previously shown to be isolated i.e. Rockall, as well as in locations previously shown to be homogeneous i.e. continental shelf areas, especially Malin Shelf in the north of Ireland.

Discussing the genetic marker results first, very low levels of diversity were observed in the mitochondrial DNA (mtDNA) COI gene in L. forbesii with only 8/425 sequences showing haplotypic variation. The majority of the samples contained the same haplotype, H1, which occurred throughout the dataset and formed the centre of a network with low levels of sequence diversity. Other genes should be examined in future studies, given the low variability seen in the COI gene, which is useful for DNA barcoding studies (Gebhardt and Knebelsberger, 2015; Maggioni et al., 2020), but relatively uninformative for studies such as this, where a faster evolving section of the mtDNA such as the control region may provide additional data. Interestingly, while other loliginids also showed low mtDNA variation (Aoki et al., 2008; Ibáñez et al., 2012), analysis of 55 Loligo vulgaris samples from Galicia in north west Spain using this gene resulted in 35 haplotypes, indicating high COI diversity in this species at that location (Garcia-Mayoral et al., 2020).

High genetic diversity was observed in the microsatellite markers at nine loci in the present study (mean H_o = 0.802, H_e = 0.866 in L. forbesii at all locations) and each microsatellite locus was highly polymorphic with 11.2–19.9 alleles per locus, depending on the location. These values are similar to previous observations in L. forbesii by Shaw et al. (1999) and they also fall within the range of values reported in Todaropsis eblanae (H_{o =} 0.83–0.93, H_e = 0.85–0.93, alleles per locus 9.6–27, Dillane et al., 2005), but are a bit higher than Sepioteuthis australis (H_o = 0.519–0.585, H_e = 0.735–0.777, alleles per locus 10.71–12, Smith et al., 2015), and Dosidicus gigas (H_o = 0.657–0.759 H_e = 0.815–0.833, alleles per locus not provided, Sanchez et al., 2020). Microsatellites revealed no statistically significant population sub-structure between sample locations in the present study, indicating high gene flow in L. forbesii across the sampled range. Nevertheless, six of the seven comparisons that were significantly different prior to Bonferroni correction involved Rockall versus various shelf locations. Given this non-random pattern, the fact that the Bonferroni correction is known to be conservative with increased Type II error (Narum, 2006), and previous results showing subtle genetic structure at Rockall, this pattern is interesting and should be considered.

Shaw et al. (1999) examined seven microsatellite markers in L. forbesii, revealing significant pairwise F_ST differences between samples from Rockall and Faroe Banks, and between Rockall and two shelf locations—west of Scotland (but only in one of two years sampled) and English Channel. Other differences in that study, i.e. between Rockall–north west Spain and Rockall–west of Scotland (in the second of the sampled years), were no longer significant after Bonferroni correction was applied (Shaw et al., 1999). Shaw et al. (1999) suggested the isolation arose due to deep water and oceanography which could present a potential barrier to isolate offshore populations like Rockall from the shelf population. Our results contrast with this: we found no statistically significant evidence of genetic sub-structure between similar sampling locations, including an additional six locations which were not sampled before; the North Sea, the northern Bay of Biscay, north, west and south Ireland and south Spain. It should be noted that there was not complete overlap between the markers used in both cases: of the nine microsatellite markers used in the present study, only five markers overlapped directly with Shaw et al. (1999) study. So although we had two additional markers in total, we also had four different markers in comparison to Shaw et al. (1999). That said, the loci which were most variable in the past regarding Rockall (Lfor 1 and Lfor3; Shaw et al., 1999) were included in the present study. Null alleles (caused by a mismatch between the primer and the target sequence resulting in alleles which fail to amplify) are known to inflate F_ST values (Chapuis and Estoup, 2007). Previous studies have shown that null alleles are common in mollusc populations with large effective population sizes (Li et al., 2003; Kaukinen et al., 2004; Astanei et al., 2005; Ramos et al., 2018). In the present study, prior to removal of loci affected by null alleles, there was one comparison involving Rockall versus northern Ireland which was significant at p < 0.001 after Bonferroni correction. But once loci which were shown to be most affected by null alleles (Lfor 8 and Lfor 12) had been removed, and once a Bonferroni correction had been applied, no significant differences remained between Rockall and the other sampling locations. Further analysis to detect genetic structure using STRUCTURE analysis and Principal Coordinate Analysis (PCoA) also showed no differentiation. Sampling of squid is opportunistic and is dependent on national groundfish survey dates but this should not have impacted on results. Samples were taken in 2019 at all locations, apart from south Ireland (2018), and all locations contained genetic samples taken in Q4 apart from Rockall, north Spain and English Channel (sampled in Q3) and North Sea (Q1) (details of sampling quarters and years are given in Supplementary Table S1).

Unlike genetic markers, differences in L. forbesii statolith shape were apparent between four out of five locations analysed by statolith shape markers. The north Ireland statoliths were highly distinct from the rest, but other locations also had pairwise differences (with the exception of Rockall versus north Scotland). What this suggests is that a) statolith shape is highly sensitive to local conditions and b) L. forbesii forms distinguishable groups (based on statolith shape statistics), maintaining these for long enough for local conditions to affect the shape of the statolith. Thus, statolith shape is an “ecological marker” in L. forbesii which can indicate a spatial structure. Rockall, north Scotland and north Ireland were all sampled for statoliths in the same year (2019, Q3/Q4) so temporal sampling variation cannot explain differences between north Ireland and these others (Rockall in Q3 was sampled only 6–8 weeks apart from north Scotland and north Ireland in Q4 in 2019). Furthermore, the remaining locations, west and south Ireland, were sampled in Q4 (2018) plus Q1 of 2020, i.e. the same “biological year” as the other locations because L. forbesii recruits in Summer/Autumn in the study area and has an approximately annual lifecycle (Collins et al., 1995; Pierce et al., 2005). Influences on statolith shape are poorly understood but appear to be environmentally induced and linked to growth of the individual. The statocyst structure permits the animal to sense acceleration and movement and to adjust its eye and head movements accordingly (Arkhipkin and Bizikov, 2000). Because demersal squid like Loligo must avoid high-velocity impact against the seafloor, they need extra control at lower accelerations, which may be provided by an oar-like rostrum in this group, that helps during “rolling plane” movements of the animal. Arkhipkin and Bizikov (2000) suggest that habitat and swimming behaviour are more important determinants of statolith shape than evolutionary relatedness. Indeed, in the present study, the L. vulgaris shape coefficients fell closer to those of L. forbesii at some locations, than did L. forbesii from north Ireland to its conspecifics in other locations (Supplementary Figure S6). Clear shape differences between sampling locations were observed in many areas of the statolith in the present study: the lateral and dorsal dome, rostral angle and wing portion.

Why statolith shape was most distinct in north Ireland from the rest is very difficult to say and possible explanations are highly speculative without experimental evidence of swimming behaviour and its effect on statolith shape. Bathymetry was not markedly different across the sampled locations, mostly spanning 50–200 m depth (Figure 1). The Malin Shelf, where “north Ireland” samples were obtained, falls under the influence of both Shelf Edge Current (SEC) and Irish Coastal Current (ICC), both of which are northward-flowing, with mean current speeds varying seasonally in the range 10–20 cm s⁻¹ and 20–25 cm s^–1, respectively (Lynch et al., 2004; Jones et al., 2018b). Hydrography in this area is complex (Jones et al., 2018b) and may elevate recruitment in zooplankton and fish (Porter et al., 2018).

Our statolith shape analysis showed differences between all locations except Rockall versus north of Scotland, close to the area sampled by Shaw et al. (1999), and somewhat at odds with the latter study, which distinguished Rockall from the west of Scotland shelf (albeit in only one of two years sampled). As we saw above, our genetic analysis also did not make Rockall distinct, although almost all the non-significant genetic trends involved Rockall. Our proposed explanation is that there appears to be a semi-isolated breeding group at Rockall. If Rockall gets periodically depleted of native stock e.g. due to high levels of directed fishing pressure (Pierce and Boyle, 2003; Young et al., 2004), or if there are movements of residents away from this area and it gets repopulated by relatively few local individuals, it could lead to low effective population size and genetic drift, assuming a degree of isolation. The number of inwards migrants, and how recent these are, might explain why sometimes Rockall is distinct (i.e. in one of two sampled years–Shaw et al., 1999), whereas at other times it is not distinct [Brierley et al., 1995; Shaw et al., 1999 (second year studied); present study]. The fact that Rockall is a small and relatively restricted area is relevant and so are the wide year-to-year fluctuations in Scottish L. forbesii landings from Rockall compared to the mainland (ICES, 2020). Large population size fluctuations tend to increase genetic drift (Frankham, 1995). The free movement of individuals between Rockall and north/west Scotland may be aided by the Rockall slope current running northwards on the eastern side of the Rockall Channel, and in the opposite direction, by southward running currents on its western side (Huthnance, 1986; Smilenova et al., 2020). Statolith shape similarities between Rockall and north Scotland may be explained by similar habitats/swimming conditions, rather than migrations and further data are required to settle this. Many population genetic studies on loliginids have shown apparently contradictory findings but these may instead be a feature of a high degree of lifecycle flexibility. This includes examples in Doryteuthis pealeii (Buresch et al., 2006; Shaw et al., 2010), L. reynaudii (Shaw et al., 2010; Van der Vyver et al., 2016) and D. opalescens (Reichow and Smith, 2001; Cheng et al., 2021). The complex picture of loliginid population genetic structure may be resolved with intensive within-season and within-spawning location sampling, incorporating advanced molecular analysis such as SNPs (single nucleotide polymorphisms), e.g. Cheng et al. (2021). Alternatively, tools besides genetic ones may be useful, or even necessary, to determine stock structure due to the flexible life cycles seen in squid.

A final piece of evidence is provided by biological data, where Rockall stands out, somewhat. Sex ratios were similar across all locations and, while the Rockall samples were statistically less mature, this was to be expected since Rockall was sampled earlier, in Q3 rather than Q4 when most other maturity data were taken. Male sizes (DML) were nevertheless on the large side at Rockall, in fact, Rockall males were significantly larger (median DML 290 mm) than males in north Scotland (median DML 190 mm), despite being less mature. Given what is known about the variability of growth, size at maturity, and the resulting polymodal size distribution in this species (e.g. Forsythe and Hanlon, 1989; Boyle et al., 1995; Forsythe, 2004), this could be indicative of better growth conditions or different growth–temperature relationships at Rockall (e.g. different seasonality in growth). Larger size at Rockall only applied to males, however, and females sizes were not different between locations. Male size distributions are complex in L. forbesii, which is also sexually dimorphic with males typically growing faster and longer than females (Rocha and Guerra, 1999). Overall, despite the statistically significant differences observed, further data are needed to provide clarity on whether there are growth (or maturity) differences at Rockall.

Information from isotopic or elemental analysis of the statolith might help to refine migratory patterns in L. forbesii in future, e.g. Green et al. (2015) working with arrow squid Nototodarus gouldi was able to place individual squid into sampling “locations” with an accuracy of up to ∼75% using statolith shape coefficients, but elemental sampling permitted additional insights such as the fact that squid with a statolith shape that was typical for a particular location were not usually “born” there, but had originated in other locations. Such observations also perfectly illustrate why genetic structure cannot be detected in many squid species and emphasise that strong ecological distinctions may be detectable using statoliths as “black boxes” (Arkhipkin, 2005).

Conclusions/management implications

We present evidence that L. forbesii squid form ecological stocks. Statolith shape statistics indicated that squid sampled in four out of five regional locations formed distinguishable groups, which indicates that they are reasonably resident over the period of time that statolith shape differences develop. Only squid in Rockall/north Scotland were statistically indistinct on the basis of statolith shape. Genetic comparisons across a geographically wider set of sampling locations did not present any statistically significant signal, although several comparisons involving Rockall came close to the threshold for statistical signifiance. We found no evidence of biological (size, maturity) differences between the ecological stocks, apart from a few differences in males, particularly at Rockall, however, this needs further study given the polymodal size structure and scope for micro-cohorts in these squid (e.g. Collins et al., 1999). Rockall, which supports a targeted fishery and appears to have larger males than other areas, may be have a semi-isolated breeding population, that is periodically colonised from shelf sources—this would explain why it is sometimes genetically distinct (see Shaw et al., 1999) and sometimes not (Brierley et al., 1995; Shaw et al., 1999; present study). Were increased management intervention to be considered in this non-quota species, we should be aware that separate ecological stocks are distinguishable at the spatial scales sampled in this study, but that a proportion of individuals in these stocks are nevertheless mobile and move between sampling locations at some or all life stages. A second consideration for management is that Rockall appears to be semi-isolated, may be dependent on colonisation from outside sources, and recolonisation may be periodic. This is important, given directed fisheries at Rockall and localised stock depletions at this location in the past (Young et al., 2004). Future research into seasonal stock definitions in Loligo forbesii, as well as presence/absence of spawning locally within Rockall, would be useful.

Data availability

The data underlying this article will be shared on reasonable request to the corresponding author.

Author contributions

AMP, ALA conceived the idea for the paper; ES, LB, EA, AL, DO, MP, CPR, JPR, IS, JV contributed samples; ES, LB collected the data; ES, LB, DOM, MP, AMP performed statistical analysis; LB, ES, AMP led the writing and all authors contributed to the writing.

Competing interests

The authors declare no competing interests.

Funding

This research comes from the Cephs & Chefs project (https://www.cephsandchefs.com/) and was funded by the European Regional Development Fund through the Interreg Atlantic Area Programme grant number EAPA_282/2016

ACKNOWLEDGEMENTS

Thank you to Mary Gannon, National University of Ireland Galway, for help with project management. The authors wish to acknowledge Carlos Veloy, Dr Kai Wieland from DTU, David Stokes and Dr Eoin Kelly from the Marine Institute, and the staff and crew of the RV Celtic Explorer for their assistance with sampling.

Author notes

co-first author