Morphological differences yet genetic similarity among seasonal cohorts of Japanese anchovy during the early life stages in the Kii Channel

Abstract

Species in the genus Engraulis show extensive intraspecific as well as interspecific morphological and genetic diversity. Since morphological differences do not necessarily correspond to genetic differences, it is necessary to clarify the relationship between morphological differences and genetic differences for a better understanding of the population structure. Fish morphology at a given standard length differs between cohorts of Japanese anchovy Engraulis japonicus during the early life stages in the Kii Channel, but it is unknown whether the differences are caused by genetic differences or not. The Kii Channel includes the boundary between the Pacific (southern side of the Kii Channel) and the Seto Inland Sea stocks (northern side), but stock separation is based primarily on demographic characteristics. In the present study, genetic analyses were conducted to examine the genetic differences among samples (month and area) based on mitochondrial DNA cytochrome b region (Cyt b), control region (CR) and microsatellite markers. AMOVA revealed that the percentage of genetic variation among samples was low at 0.11% for Cyt b, 0.30% for CR and 0.00% for microsatellite, and no significant genetic variation was observed among samples. Although two clades were identified in the unrooted neighbour-joining tree for Cyt b and CR, both Cyt b and CR sequences were similar between months and between areas. Accordingly, the morphological differences among cohorts can be attributed to phenotypic plasticity. Additionally, there were no genetic differences between samples from the southern side and the northern side of the Kii Channel, suggesting strong genetic connectivity in these areas.

Introduction

Fish species identification is generally based on morphological characteristics and DNA sequence data. It is important to identify the fish to species accurately, since confounding multiple species in a stock assessment would amount to ignoring the different population dynamics of each individual species (Karahan et al., 2014). Fish stocks, the spatial management unit, are commonly defined on the basis of genetic structure, morphological and demographic characteristics, fishing patterns, and/or connectivity patterns (adult movement and larval dispersal) (Cuttitta et al., 2015; Huret et al., 2020). As such, morphological characteristics and genetic structure would provide useful stock assessment information for fisheries management.

Species of the genus Engraulis are widely distributed in marine ecosystems throughout the world (Whitehead et al., 1988), and ecologically and economically important resources (Zhang & Xian, 2015; Rioual et al., 2021; Khan et al., 2022; Michishita et al., 2023). Species in the genus Engraulis show extensive intraspecific as well as interspecific morphological and genetic diversity. In Mersin Bay, eastern Mediterranean Sea, two adult Engraulis forms differentiated by their external morphology and designated as silver anchovy (habitat depth of 9–77 m) and blue anchovy (majority at depth of 69–111 m) showed strong genetic differences (Karahan et al., 2014). Additionally, morphological differences between species in the genus Engraulis from coastal and offshore habitats have been shown to have a heritable genetic basis in the north-eastern Atlantic, the western Mediterranean and the Black Sea (Bonhomme et al., 2021). On the other hand, larval Engraulis encrasicolus showed habitat-specific morphological variations and lack of corresponding genetic variation among larvae from different locations in the Strait of Sicily (Cuttitta et al., 2015). In the northeastern Pacific Ocean, Engraulis mordax is typically divided into three subpopulations based on electrophoretic, morphometric and meristic characteristics of adult fish, but there is no apparent genetic structure (Lecomte et al., 2004; Schwartzkopf et al., 2022). In Japanese anchovy Engraulis japonicus distributed in the western North and central Pacific (Whitehead et al., 1988), two types of adult Japanese anchovy have been differentiated based on the condition factor (fat and normal types), but no genetic difference was observed between these two types in Sagami Bay, Japan (Nakategawa et al., 2009). On the other hand, larval Japanese anchovy have been shown to have spatiotemporal differences in morphology among populations, and the differences suggested to be genetically based on surveys in the southern East China Sea (Chen et al., 2010). As such, since morphological differences do not necessarily correspond to genetic differences, it is necessary to clarify the relationship between morphological differences and genetic differences for a better understanding of the population structure.

Japanese anchovy has an extensive spawning area and a protracted spawning season with larvae present in the water column throughout the year (Yu et al., 2002; Takasuka et al., 2005; Oozeki et al., 2007; Zheng et al., 2015). Therefore, throughout the year multiple cohorts of larvae from different spawning events occur (Wang & Tzeng, 1999; Yasue et al., 2006, 2011). In genetic studies on the population structure, Yu et al. (2002) reported that a significant population differentiation was found between two spawning populations, one off northeastern Taiwan in the Pacific Ocean and the other off southwestern Taiwan in the Taiwan Strait, but not between the temporal (two consecutive years) populations in the northeastern Taiwan in the Pacific Ocean, although the level of geographic differentiation was weak. Chen et al. (2010) revealed a significant genetic structure among three spatiotemporally distinct larval groups in the southern East China Sea (eastern Taiwan and Taiwan Strait). In contrast, Yu et al. (2005) and Zheng et al. (2015) reported that the existence of separate genetic stocks was not detected in the Yellow Sea and East China Sea. Liu et al. (2006) reported that no significant population genetic structure occurred in the East China Sea and on the coast of Japan. Funamoto et al. (1999) reported that no separate genetic stocks were found for the coast of Japan and in the Kuroshio-Oyashio transition region. As such, although several studies have examined the population structure of Japanese anchovy genetically, it is necessary to investigate different cohorts to clarify the population structure in more detail.

Two species of Engraulidae occur in the Kii Channel, Japan: Japanese anchovy and Encrasicholina punctifer, with Japanese anchovy being dominant. During the early life stages, the morphology of Japanese anchovy is similar to that of E. punctifer, but differences in melanophores on the ventral midline of the caudal peduncle enable rapid species identification based on morphology (Noichi, 2014). Larval Japanese anchovy occur throughout the year in the Kii Channel, and they are caught from this region by the commercial ‘shirasu’ fishery. The shirasu fishery efficiently catches larvae and juveniles of small pelagic fish such as Japanese anchovy, Sardinops melanostictus and Etrumeus micropus in the shallow waters for food resources (called shirasu in Japanese), and this fishery is economically important in Japan (Kono et al., 2003; Xie & Watanabe, 2007). Fish morphology at a given standard length (developmental stage) differs among cohorts during the early life stages in relation to sea temperature (Yasue et al., 2016). However, it is unknown whether this difference is caused by genetic differences or phenotypic plasticity. Since fish morphology can show marked differences among cohorts, it is often pointed out by fishermen that a certain portion of fish may represent a separate species such as E. punctifer. The Kii Channel has two environmentally distinct areas (Yoshioka, 1988). The southern side (outer area, OA) has oceanic conditions typical of an open sea area, such as deep waters, minimal effect of river discharge, and frequent intrusion of the Kuroshio branch currents. In contrast, the northern side (inner area, IA) has environmental conditions typical of a coastal sea area, such as shallow waters, large river discharge and mixing with the Seto Inland Sea water. In OA, Japanese anchovy spawn mainly from February to September with peak spawning in March, while in IA Japanese anchovy spawn mainly from April to September with peak spawning in June (Yasue et al., 2006). If the spawning season is fixed for each cohort, genetic differences may occur among cohorts. For stock assessment in waters around Japan, Japanese anchovy is assumed to constitute three separate stocks: Pacific, Seto Inland Sea, and Tsushima Warm Current stocks. This stock separation is based primarily on demographic characteristics. The Kii Channel includes the boundary between the Pacific and the Seto Inland Sea stocks (Kinoshita et al., 2022; Kono et al., 2022). Japanese anchovy distributed in OA has been included as part of the Pacific stock, and that in IA as part of the Seto Inland Sea stock. Larvae hatched in OA could disperse/migrate to IA as they grow in spring (Yasue et al., 2008), although no details of dispersal/migration are known for Japanese anchovy. Larvae that are distributed in IA in spring and during summer to autumn have been clearly distinguished as different cohorts (Kono et al., 2022). If Japanese anchovy migrate to the areas where they were born during the spawning season, genetic differences may exist among cohorts in OA and IA. Here emerged the need to clarify the issue of morphological and genetic differences in Japanese anchovy in the Kii Channel, where the population structure is controversial in terms of cohorts, dispersal/migration, and oceanographic characteristics. In other words, Japanese anchovy in the Kii Channel is an ideal model to understand the relationship between morphological characteristics and genetic structure in small pelagic fish.

Mitochondrial DNA cytochrome b region (Cyt b), control region (CR) and microsatellite loci have provided powerful markers for phylogenetic and population structure studies of fishes (e.g. Egan et al., 2018; Caballero-Huertas et al., 2022). In species of the genus Engraulis, Japanese anchovy has been discriminated from E. encrasicolus through comparison of Cyt b sequence data (Grant et al., 2005; Nakategawa et al., 2009). Cyt b (Yu et al., 2005; Zheng et al., 2015), CR (Liu et al., 2006) and microsatellite (Yu et al., 2002) markers were used to examine population structure of Japanese anchovy. Cyt b (Silva et al., 2014b), CR (Viñas et al., 2014; Cuttitta et al., 2015) and microsatellite (Zarraonaindia et al., 2009; Ruggeri et al., 2016; Ouazzani et al., 2017) markers as well as the combination of them (Borrell et al., 2012; Oueslati et al., 2014; Silva et al., 2014a) have also been used to examine the population structure of other species of the genus Engraulis. These markers are effective for identifying species and detailing the population structure.

In the present study, species identification was performed using Cyt b sequence data for fish identified as Japanese anchovy based on morphology to genetically examine whether fish occurring throughout the year are Japanese anchovy in the Kii Channel. Then, the population structure of Japanese anchovy was analysed based on Cyt b and CR sequence data and microsatellite data to test whether the difference in morphology among cohorts is caused by genetic differences and to understand the genetic connectivity between OA and IA.

Materials and methods

Field survey and sea temperature data

Japanese anchovy were collected by commercial fisheries in OA and IA during 15 April 2019 to 4 August 2022 (Table 1; Figure 1). The shirasu pair trawl (Yasue & Takasuka, 2009) and the lift-net fisheries usually target shoals of fish with standard length (SL) of ca. 13–30 (larvae and juveniles) and 50–150 (juveniles and adults), respectively. Sampling for species identification was conducted in OA and IA throughout the year (sampling events: 25, N = 3–32 for each sample), including adult fish (sample codes H, I, and R) for comparison with fish during the early life stages. The relationships between SL and head length (HL) and between SL and the degree of guanine deposition on the peritoneal surface differ among cohorts during the early life stages, and inter-cohort differences are clear between the high- and low-sea temperature periods (Yasue et al., 2016). Therefore, sampling was conducted in September (sample codes J, T, and V) and February (sample codes O, U, and W) in OA and in April (sample code X) and August (sample code Y) in IA for morphological and genetic analyses (sampling events: 8, N = 31–32 for each sample). Fish collected in the four months are regarded as different seasonal cohorts in the Kii Channel (March–May: spring, June–August: summer, September–November: autumn, December–February: winter). The samples (sample codes J, T, V, O, U, W, X, and Y) were also used to compare genetic data between OA and IA. Fish samples were preserved in chilled sea water after sampling.

Table 1. Sample profiles of Japanese anchovy for morphological and genetic analyses in the Kii Channel

^a The same fish were used in the morphological analysis.

^b One fish was not sequenced for mitochondrial DNA control region.

Cyt b, mitochondrial DNA cytochrome b region; CR, mitochondrial DNA control region; ST, sea temperature at 10 m depth.

Figure 1. Study areas in the Kii Channel. Shaded areas, solid circles, and the bold line represent sampling areas for Japanese anchovy, sea temperature stations, and the boundary between the southern side (outer area, OA) and the northern side (inner area, IA), respectively.

Sea temperature (ST, °C) data of OA were obtained from records of the ST survey measured by self-registering thermometer at a depth of 10 m attached to the observation tower (Shirahama Oceanographic Observatory, Disaster Prevention Research Institute, Kyoto University, Wakayama, Japan). ST data of IA were obtained from records of the ST survey of Wakayama Prefectural Fisheries Experimental Station measured by self-registering thermometer located at a depth of 10 m attached to a fixed commercial set net.

Morphological analysis

Broad size ranges of Japanese anchovy during the early life stages were selected from the daily catch. Species identification was made on the basis of morphology: melanophores on the ventral midline of the caudal peduncle (Noichi, 2014). SL and HL of each fish were measured to the nearest 0.1 mm using digital calipers (under a stereo microscope in the case of ⩽30 mm SL). The degree of guanine deposition was observed using a stereo microscope at a magnification of ⩽63 × . Fish were assigned to two developmental stages, as follows. GD-0 stage: no guanine deposition observed on the peritoneal surface (lateral part of the abdominal cavity near the pelvic fin base); GD-1 stage: guanine deposition observed on the peritoneal surface (Yasue et al., 2016).

The relationships between SL and HL and between SL and the degree of guanine deposition were examined for each area (OA and IA). An analysis of covariance (ANCOVA) was conducted to compare the difference of ln-transformed HL on ln-transformed SL between September (sample codes J, T, and V) and February (sample codes O, U, and W) in OA. ANCOVA was also conducted to compare the difference of ln-transformed HL on ln-transformed SL between April (sample code X) and August (sample code Y) in IA.

Mitochondrial DNA analysis

After the morphological analysis, a piece of muscle including skin was extracted from each fish. The piece of muscle was preserved in 99.5% ethanol until the subsequent DNA analysis. Total DNA was isolated using DNeasy blood & tissue kit (QiaGen). Cyt b (1141 bp) and CR (1024 bp) were amplified by the polymerase chain reaction (PCR) using KOD FX (Toyobo) and primer pairs (Funamoto et al., 1999; Inoue et al., 2000; Silva et al., 2014a; Table 2). Amplification reactions were performed in a total volume of 20 μl containing 4.4 μl of template DNA, 1.2 μl of each primer (10 μM), 10 μl of 2× PCR buffer for KOD FX, 4 μl of 2 mM dNTPs and 0.4 μl of KOD FX. The thermal cycle consisted of an initial denaturation at 94°C for 2 min, followed by 35 cycles of 98°C for 10 s, 52°C for 30 s and 68°C for 1 min. The PCR products were purified AMPure XP beads (Beckman Coulter), and cycle sequence reactions were performed using BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems) and sequence primers. Cycle sequencing reactions were performed in a total volume of 10.35 μl containing 7.65 μl of template DNA, 0.5 μl of each primer (10 μM), 0.4 μl of BigDye Terminator 3.1 Ready Reaction Mix and 1.8 μl of 5× Sequencing Buffer. The thermal cycle consisted of an initial denaturation at 96°C for 1 min, followed by 30 cycles of 96°C for 10 s, 50°C for 5 s and 60°C for 4 min. Then, PCR products were sequenced using a 3730xL DNA analyser (Applied Biosystems).

Table 2. Primer nucleotide sequences for mitochondrial DNA cytochrome b region, control region and microsatellite markers in Japanese anchovy

#, fluorescent-dye.

DNA sequences were aligned using the CLUSTAL W program implemented in the MEGA X software (Thompson et al., 1994; Kumar et al., 2018). Basic Local Alignment Search Tool (BLAST) analysis (https://blast.ncbi.nlm.nih.gov/Blast.cgi; accessed online 13 August 2024) was conducted to confirm identities with Cyt b sequences of Japanese anchovy NC_003097.1 (Inoue et al., 2001), KF765500.2 (Zhang & Xian, 2015), and AP017957.1 (Lavoué et al., 2017). BLAST analysis was also conducted to confirm identities with the other species of the genus Engraulis: E. encrasicolus NC_009581.1 (Lavoué et al., 2007) and E. mordax NC_041097.1 (Lewis & Lema, 2019). For Cyt b and CR sequence data, number of haplotypes, haplotype diversity and nucleotide diversity were calculated for each sample (sample codes J, O, V, W, X, and Y) using the DnaSP version 6 software (Rozas et al., 2017) to examine genetic diversity. An analysis of molecular variance (AMOVA) considering the six samples as a single group and a pairwise F_ST analysis were conducted using the Arlequin ver 3.5.2.2 software (Excoffier & Lischer, 2010) to examine genetic differences among samples. The level of statistical significance was set at P < 0.05 after Bonferroni correction for a pairwise F_ST analysis. Genetic distances were calculated based on the Kimura-2 parameter (K2P) model (Kimura, 1980), and the unrooted neighbour-joining (NJ) tree was constructed using the MEGA X software (Kumar et al., 2018) to examine phylogenetic relationships among fish. The confidence in tree nodes was assessed with 1000 bootstraps.

Microsatellite DNA analysis

For the total DNA isolated from each fish, six microsatellite loci developed by Chiu et al. (2002) were amplified by PCR using Type-it Microsatellite PCR (QiaGen) and fluorescent dye-labelled primers (Table 2). The six loci were amplified in four independent PCR: 3-plex reaction for EJ2, EJ35 and EJ41.1; independent reactions for EJ9, EJ27.1 and EJ27.2. In the 3-plex reaction for EJ2, EJ35 and EJ41.1, amplification reactions were performed in a total volume of 10 μl containing 4 μl of template DNA, 1 μl primer (EJ2: 2 μM; EJ35: 1.6 μM; EJ41.1: 1.2 μM), and 5 μl of 2x Type-it Multiplex PCR Master Mix. The thermal cycle consisted of an initial denaturation at 95°C for 5 min, followed by 31 cycles of 95°C for 30 s, 60°C for 90 s and 72°C for 30 s. In the independent reactions for EJ9, EJ27.1 and EJ27.2, amplification reactions were performed in a total volume of 10 μl containing 4 μl of template DNA, 1 μl of each primer (2 μM) and 5 μl of 2x Type-it Multiplex PCR Master Mix. The thermal cycle consisted of an initial denaturation at 95°C for 5 min, followed by 35 cycles of 95°C for 30 s, 54°C for 90 s and 72°C for 30 s. Then, PCR products were sequenced using a 3730xL DNA analyser. Allele types were determined using the GeneMapper version 4.1, 5.0 or 6.0 software (Applied Biosystems). However, EJ35 was excluded from the following data analysis, because it was not possible to determine the allele types.

Number of alleles, observed heterozygosity and expected heterozygosity were calculated for each locus of each sample (sample codes T, U, V, and W) to examine genetic diversity. Hardy-Weinberg equilibrium (HWE) was tested for each locus of each sample. AMOVA considering the four samples as a single group and a pairwise F_ST analysis were conducted to examine the genetic differences among samples. The level of statistical significance was set at P < 0.05 after Bonferroni correction for a pairwise F_ST analysis. The GenALEx 6.5 software (Peakall & Smouse, 2006, 2012) was used to conduct these analyses.

Results

Fish morphology

A significant difference in ln-transformed HL was revealed between September (sample codes J, T, and V) and February (sample codes O, U, and W) in OA (ANCOVA, df = 1, 189, F = 618.39, P < 0.001; Figure 2). HL was larger and therefore the trunk length (SL minus HL) was smaller at the same SL in September. ANCOVA also revealed a significant difference in ln-transformed HL between April (sample code X) and August (sample code Y) in IA (ANCOVA, df = 1, 60, F = 388.57, P < 0.001). HL was larger and therefore the trunk length was smaller at the same SL in August.

Figure 2. Relationships between ln-transformed standard length (SL) and ln-transformed head length (HL) for Japanese anchovy in the southern side (outer area, OA) and the northern side (inner area, IA) of the Kii Channel. September (Sample codes J, T, and V): ln HL = –2.983 + 1.474 ln SL (N = 96, R ² = 0.842, P < 0.001); February (Sample codes O, U, and W): ln HL = –3.630 + 1.603 ln SL (N = 96, R ² = 0.936, P < 0.001); April (Sample code X): ln HL = –3.630 + 1.587 ln SL (N = 32, R ² = 0.973, P < 0.001); August (Sample code Y): ln HL = –3.331 + 1.555 ln SL (N = 31, R ² = 0.946, P < 0.001).

In OA, the maximum SL of a GD-0 stage fish was 26.9 mm, and GD-1 stage fish occurred from 25.2 mm in September (sample codes J, T, and V; Figure 3). The maximum SL of a GD-0 stage fish was 36.5 mm, and GD-1 stage fish occurred from 32.5 mm in February (sample codes O, U, and W). Guanine deposition started at a smaller SL in September compared with February. In IA, the maximum SL of a GD-0 stage fish was 28.2 mm, and GD-1 stage fish occurred from 26.9 mm in August (sample code Y). The maximum SL of a GD-0 stage fish was 37.0 mm, and GD-1 stage fish occurred from 36.6 mm in April (sample code X). Guanine deposition started at a smaller SL in August compared with April. As such, fish morphology at standard length differed between months in both OA and IA.

Figure 3. Relationships between standard length and guanine deposition for Japanese anchovy in the southern side (outer area, OA) and the northern side (inner area, IA) of the Kii Channel. Guanine deposition 0: no guanine deposition observed on the peritoneal surface (GD-0 stage); Guanine deposition 1: guanine deposition observed on the peritoneal surface (GD-1 stage). September: Sample codes J, T, and V; February: Sample codes O, U, and W; April: Sample code X; August: Sample code Y.

Mitochondrial DNA cytochrome b region

All fish were successfully sequenced for Cyt b (N = 242). As a result of BLAST analysis, the Cyt b sequences exhibited high identity with those of Japanese anchovy registered in GenBank. BLAST analysis revealed 97.81–99.74% identity with Cyt b sequences of Japanese anchovy NC_003097.1, 97.55–99.47% with KF765500.2 and 97.81–99.47% with AP017957.1. On the other hand, the Cyt b sequences exhibited 95.44–96.76% identity with those of E. encrasicolus NC_009581.1 and 84.05–85.19% with E. mordax NC_041097.1.

Haplotype diversity ranged from 0.998 to 1.000, and nucleotide diversity ranged from 0.008 to 0.011 for Cyt b (Table 3). Genetic diversity of the six samples (sample codes J, O, V, W, X, and Y) was high.

Table 3. Genetic statistics for the mitochondrial DNA of Japanese anchovy in the Kii Channel

N, number of individuals; n_h, number of haplotypes; h, haplotype diversity; π: nucleotide diversity.

AMOVA revealed that the percentage of genetic variation among samples was 0.11%, and no significant genetic variation was observed among samples (P = 0.326; Table 4). Pairwise F_ST analysis revealed that F_ST values were low ranging from −0.006 to 0.014, and no significant genetic variation was observed between samples (P ⩾ 0.05 after Bonferroni correction; Table 5).

Table 4. Results of an analysis of molecular variance based on mitochondrial DNA for Japanese anchovy in the Kii Channel

Table 5. Pairwise F_ST and P values between samples based on mitochondrial DNA for Japanese anchovy in the Kii Channel

F_ST values are above the diagonal, and P values are below the diagonal.

K2P genetic distances between individuals ranged from 0 to 2.78%. Two clades (clades 1 and 2) were identified in the unrooted NJ tree, although bootstrap values were not high (Figure 4). The frequency of the clades 1 and 2 were 49.2% (N = 119) and 50.8% (N = 123), respectively (Table 6). In the monthly comparison, there were no significant differences in the frequency of clades 1 and 2 in both September (sample codes J and V; P = 0.264) and February (sample codes O and W; P = 0.527) of OA. There were also no significant differences in the frequency of clades 1 and 2 in both April (sample code X; P = 0.925) and August (sample codes Y; P = 0.930) of IA. In comparison between areas, there were no significant differences in the frequency of clades 1 and 2 in both OA (sample codes J, O, V and W; P = 0.731) and IA (sample codes X and Y; P = 0.996). As such, Cyt b sequences were similar between months and between areas.

Figure 4. Unrooted neighbour-joining tree showing the relationships of sequences of mitochondrial DNA cytochrome b region and control region for Japanese anchovy in the Kii Channel. Bootstrap values above 90% are shown in the tree.

Table 6. χ ² test on the differences of ratios between clade 1 and clade 2

Mitochondrial DNA control region

241 fish were sequenced for CR, although one fish (sample code X) was not sequenced. Haplotype diversity was 1.000, and nucleotide diversity ranged from 0.010 to 0.011 for CR (Table 3). Genetic diversity of the six samples (sample codes J, O, V, W, X, and Y) was high.

AMOVA revealed that the percentage of genetic variation among samples was 0.30%, and no significant genetic variation was observed among samples (P = 0.201; Table 4). Pairwise F_ST analysis revealed that F_ST values were low ranging from −0.004 to 0.015, and no significant genetic variation was observed between samples (P ⩾ 0.05; Table 5).

K2P genetic distances between individuals ranged from 0 to 2.59%. Two clades (clades 1 and 2) were identified in the unrooted NJ tree, although bootstrap values were not high (Figure 4). The frequency of the clades 1 and 2 were 48.1% (N = 116) and 51.9% (N = 125), respectively (Table 6). As in case of Cyt b, there were no significant differences in the frequency of clades 1 and 2 for months and areas (P ⩾ 0.05). CR sequences were similar between months and between areas.

Microsatellite DNA

Number of alleles ranged from 16 (EJ41.1 in sample code T) to 41 (EJ9 in sample code T) (Table 7). Observed heterozygosity ranged from 0.548 (EJ27.2 in sample code U) to 0.938 (EJ9 in sample code W). Significant departures from HWE over loci EJ2, EJ27.1, and EJ27.2 were observed in the four samples (sample code T, U, V, and W; P < 0.05). On the other hand, locus EJ9 was at HWE in the four samples (P ⩾ 0.05). Significant departures from HWE were observed in the sample code U (P < 0.05), but the others were at HWE for locus EJ41.1 (P ⩾ 0.05).

Table 7. Genetic statistics for five microsatellite loci of Japanese anchovy in the Kii Channel

N, number of individuals; N_A, number of alleles; H_O, observed heterozygosity; H_e, expected heterozygosity; P, P value for the test of Hardy-Weinberg equilibrium; *<0.05.

AMOVA revealed that the percentage of genetic variation among samples was 0.00%, and no significant genetic variation was observed among samples (P = 0.817; Table 8). Pairwise F_ST analysis revealed that F_ST values were low ranging from −0.004 to 0.015, and no significant genetic variation was observed between samples (P ⩾ 0.05 after Bonferroni correction; Table 9).

Table 8. Results of an analysis of molecular variance based on five microsatellite loci for Japanese anchovy in the Kii Channel

* <0.05.

Table 9. Pairwise F_ST and P values between samples based on five microsatellite loci for Japanese anchovy in the Kii Channel

F_ST values are above the diagonal, and P values are below the diagonal.

Discussion

The present study challenges a controversial issue of morphological and genetic differences in Japanese anchovy in the Kii Channel, through a combination of morphological analysis and genetic analysis. This approach was enabled by sampling fish throughout the year. The fish samples of the present study are mainly categorized in the transition period from larvae to juveniles (Moreno et al., 2011). The use of samples during the early life stages could be a bias source for genetic population structure studies, because samples might include fish of the same brood (Chen et al., 2010; Katamachi et al., 2015). However, the commercial shirasu pair trawl targets shoals of fish, and the catches are large (>20 kg) enough to include multiple shoals of fish. Additionally, wide size ranges of fish were selected from the daily catch. Therefore, the fish samples were assumed to be randomly collected from sea areas for genetic analysis.

Seasonal differences were detected in the relationships of head length and guanine deposition to standard length in OA and IA. In the Strait of Sicily, the most useful characteristics to discriminate the different larval groups were identified as the length of the mouth and the body diameter for larval E. encrasicolus (ca. 5.0–5.8 mm SL) (Cuttitta et al., 2015). Body height has been suggested to use for group discrimination in addition to body size for larval Japanese anchovy (19.7–36.9 mm SL) in the southern East China Sea (Chen et al., 2010). HL and guanine deposition were adopted as morphological characteristics in the present study, because head size, degree of guanine deposition and ossification at a given standard length differ among cohorts in the Kii Channel (Yasue et al., 2016) and also because this method does not need the double staining process to prepare skeletal specimens for each fish. Environmental conditions such as sea temperature and food availability affect the degree of guanine deposition of Japanese anchovy (Takahashi & Watanabe, 2004). Although inter-cohort variability in morphological characteristics is influenced by changes in environmental conditions, the present study focused on the morphological differences observed at the same SL and tested whether the differences among seasonal cohorts are caused by genetic differences or not. The data plots (values) were continuously distributed in the relationships between ln-transformed SL and ln-transformed HL (Figure 2), but the plots showed low overlapping between September and February in OA and between April and August in IA. Despite the overlaps among samples, morphologically distinct groups were found in E. encrasicolus (Caneco et al., 2004; Khan et al., 2022). A possible difference in morphological characteristics between areas (OA and IA) was not examined, since samples were collected from different months for morphological and genetic analyses in the present study. In any case, the difference can be assumed to be small compared to the difference in morphological characteristics among months.

In the genetic analysis, Cyt b, CR and microsatellite (EJ2, EJ9, EJ27.1, EJ27.2, EJ41.1) markers were used to determine the population structure of Japanese anchovy. 13 out of 20 cases deviated from HWE in the microsatellite markers, indicating heterozygote deficiency. One possible reason for the observed heterozygote deficiency is the existence of null alleles. Mitochondrial DNA 16S (Borrell et al., 2012) and CO (Yu et al., 2005) regions, single nucleotide polymorphisms (Zarraonaindia et al., 2012; Huret et al., 2020), nuclear DNA intron sequences (Borsa et al., 2004; Bouchenak-Khelladi et al., 2008), allozyme analysis (Sanz et al., 2008) and restriction fragment length polymorphism analysis of mitochondrial DNA (Magoulas et al., 2006) have been also used for population structure studies of species of the genus Engraulis. Since similar results were obtained for Cyt b, CR and microsatellite markers in the present study, these markers would be valid to determine the population structure of Japanese anchovy using these markers.

The Cyt b sequences exhibited high identity with those of Japanese anchovy registered in GenBank, showing ⩾97.55% identity with Cyt b sequences of NC_003097.1, KF765500.2, and AP017957.1. Additionally, the Cyt b sequences exhibited relatively low identity with those of the other species of the genus Engraulis, showing ⩽96.76% identity. Accordingly, all the 242 fish were genetically identified as Japanese anchovy. Japanese anchovy was distinguished from E. encrasicolus and E. mordax based on the Cyt b sequences. Thus, species identification based on the morphological characteristics was consistent with that based on the genetic markers. It was genetically confirmed that Japanese anchovy during the early life stages occurs throughout the year in the Kii Channel.

In Japanese anchovy, haplotype diversity of Cyt b was 0.950–0.981 (Yu et al., 2005) and 0.4706–0.9935 (Zheng et al., 2015) in the Yellow Sea and East China Sea, and 0.995–1.00 (Chen et al., 2010) in the southern East China Sea. Haplotype diversity of CR was 0.99–1.00 in the northeastern Taiwan in Pacific Ocean and southwestern Taiwan in the Taiwan Strait (Liu et al., 2006). Observed heterozygosity of microsatellite loci (EJ2, EJ41.1, EJ9, EJ27.1, EJ27.2) was 0.317 (EJ2)–0.783 (EJ41.1) near northeastern Taiwan in the Pacific Ocean and near southwestern Taiwan in the Taiwan Strait (Yu et al., 2002). The values of haplotype diversity (0.998–1.000) and observed heterozygosity (0.548–0.938) in the Kii Channel were comparable to those in these areas. AMOVA revealed that the percentage of genetic variation among samples was low at 0.11% for Cyt b, 0.30% for CR and 0.00% for microsatellite, and no significant genetic variation was observed among samples. Pairwise F_ST analysis revealed no significant genetic variation between samples. Two clades were identified in the unrooted NJ tree for Cyt b and CR, which suggests the possibility that there is genetic structure. However, both Cyt b and CR sequences were similar between months and between areas. It should be noted that although the unrooted NJ tree was constructed, the confidence of tree nodes in phylogenetic relationships among fish was not high. Therefore, we concluded that there were no genetic differences among cohorts in the Kii Channel.

Although fish morphology at a given standard length differed among seasonal cohorts, there were no genetic differences in Japanese anchovy in the Kii Channel. Their morphological differences among seasonal cohorts would be most reasonably attributable to phenotypic plasticity, as follows. Japanese anchovy exhibited substantial variability in biological characteristics that can be attributed to phenotypic plasticity in previous studies. Since two types of adult Japanese anchovy differentiated by the condition factor showed no genetic difference, their morphological differences could be due to the differences in their migration pattern and food environments (Nakategawa et al., 2009). In captive conditions, the size and nutrient content of eggs and batch intervals were highly constrained by temperature (Yoneda et al., 2014). Additionally, the degree of guanine deposition tended to be accelerated and growth rates were higher with increasing temperature and food availability during the early life stages (Takahashi & Watanabe, 2004). In other species of the genus Engraulis, for example, E. encrasicolus larvae showed morphological variations and lack of corresponding genetic variation in the Strait of Sicily, and the observed differences in morphology were linked to environmental conditions such as food availability (Cuttitta et al., 2015). In the Kii Channel, there is substantial seasonal variability in environmental conditions such as sea temperature and plankton density (Yasue & Takasuka, 2009). Seasonal variability in sea temperature was considerable, ranging from 14.3 to 28.3°C in the present study, although food availability was not quantitatively estimated. Inter-cohort differences of anchovy morphology were clear between the high- and low-sea temperature periods (Yasue et al., 2016). Differences in sea temperature among seasons are considered to be a major factor causing the morphological differences of Japanese anchovy during the early life stages. In fact, seasonal variability in larval growth rate was explained by sea temperature variability (Yasue & Takasuka, 2009). From the above, the observed morphological differences among seasonal cohorts can be attributed to the environmental variability, that is, phenotypic plasticity. The differences in fish morphology may reflect a survival strategy during the early life stages of Japanese anchovy. For example, the reflection of light from the surface (guanine) of fish can serve as camouflage (Denton & Nicol, 1965) potentially reducing risks of detection by predators and offering an advantage to fish at warmer sea temperatures. In defining fish stocks, morphological variation is influenced by environmental conditions, and therefore may not be very useful in comparing groups that experience different environmental conditions (Kristoffersen & Magoulas, 2008).

There were no genetic differences between OA and IA, suggesting strong genetic connectivity in these areas. In general, lower levels of differentiation between marine fish populations are attributed to higher dispersal potential during planktonic egg, larval, or adult stages in the life history, coupled with an absence of physical barriers to movement between ocean basins or adjacent continental margins (Grant & Bowen, 1998). Although environmental conditions are different between OA and IA (Yoshioka, 1988), such differences would not be a substantive barrier to gene flow of Japanese anchovy. The population structure of adult Japanese anchovy was not examined, but the present results were consistent with previous studies on adult Japanese anchovy (Funamoto et al., 1999; Liu et al., 2006; Zheng et al., 2015): the existence of separate genetic stocks was not detected in waters around Japan. To summarize the knowledge about the population structure of Japanese anchovy, population differentiation was reported in waters around Taiwan (Yu et al., 2002; Chen et al., 2010), but no genetic differences have been observed in other areas, including coastal and offshore areas (Funamoto et al., 1999; Liu et al., 2006; Zheng et al., 2015). We concluded that there were no genetic differences between the Pacific and Seto Inland Sea stocks of Japanese anchovy in the Kii Channel.

In conclusion, the combination of morphological and genetic analyses showed that the distinct morphological differences among the seasonal cohorts of Japanese anchovy are due to high phenotypic plasticity rather than genetic differences and that their cohorts are genetically connected around the border of the two stocks in the Kii Channel. Since fish stocks are commonly defined on the basis of various factors (Cuttitta et al., 2015; Huret et al., 2020), the present study provides useful information on the population structure for defining stocks of Japanese anchovy in Japan.