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., Reference Karahan, Borsa, Gucu, Kandemir, Ozkan, Orek, Acan, Koban and Togan2014). 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., Reference Cuttitta, Patti, Maggio, Quinci, Pappalardo, Ferrito, De Pinto, Torri, Falco, Nicosia, Musco, Armeri, Placenti, Tranchida, Mifsud, Bonanno and Mazzola2015; Huret et al., Reference Huret, Lebigre, Iriondo, Montes and Estonba2020). 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., Reference Whitehead, Nelson and Wongratana1988), and ecologically and economically important resources (Zhang & Xian, Reference Zhang and Xian2015; Rioual et al., Reference Rioual, Ofelio, Rosado-Salazar, Dionicio-Acedo, Peck and Aguirre-Velarde2021; Khan et al., Reference Khan, Bal, Battal and Seyhan2022; Michishita et al., Reference Michishita, Gibble, Tubbs, Felton, Gjeltema, Lang and Finkelstein2023). 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., Reference Karahan, Borsa, Gucu, Kandemir, Ozkan, Orek, Acan, Koban and Togan2014). 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., Reference Bonhomme, Meyer, Arbiol, Bănaru, Bahri-Sfar, Fadhlaoui-Zid, Strelkov, Arculeo, Soulier, Quignard and Gagnaire2021). 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., Reference Cuttitta, Patti, Maggio, Quinci, Pappalardo, Ferrito, De Pinto, Torri, Falco, Nicosia, Musco, Armeri, Placenti, Tranchida, Mifsud, Bonanno and Mazzola2015). 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., Reference Lecomte, Grant, Dodson, Rodríguez-Sánchez and Bowen2004; Schwartzkopf et al., Reference Schwartzkopf, Dorval, James, Walker, Snodgrass, Porzio and Erisman2022). In Japanese anchovy Engraulis japonicus distributed in the western North and central Pacific (Whitehead et al., Reference Whitehead, Nelson and Wongratana1988), 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., Reference Nakategawa, Dairiki, Itoi, Sugita and Akimoto2009). 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., Reference Chen, Tzeng and Chiu2010). 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., Reference Yu, Lee, Huang and Chiu2002; Takasuka et al., Reference Takasuka, Oozeki, Kubota, Tsuruta and Funamoto2005; Oozeki et al., Reference Oozeki, Takasuka, Kubota and Barange2007; Zheng et al., Reference Zheng, Zou and Han2015). Therefore, throughout the year multiple cohorts of larvae from different spawning events occur (Wang & Tzeng, Reference Wang and Tzeng1999; Yasue et al., Reference Yasue, Utsumi and Takeda2006, Reference Yasue, Takasuka and Shirakihara2011). In genetic studies on the population structure, Yu et al. (Reference Yu, Lee, Huang and Chiu2002) 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. (Reference Chen, Tzeng and Chiu2010) 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. (Reference Yu, Kong, Guo, Jiang, Zhuang and Jin2005) and Zheng et al. (Reference Zheng, Zou and Han2015) reported that the existence of separate genetic stocks was not detected in the Yellow Sea and East China Sea. Liu et al. (Reference Liu, Gao, Zhuang, Jin, Yokogawa and Zhang2006) reported that no significant population genetic structure occurred in the East China Sea and on the coast of Japan. Funamoto et al. (Reference Funamoto, Kikuchi, Aoki, Watabe and Taniuchi1999) 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, Reference Noichi and Okiyama2014). 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., Reference Kono, Tsukamoto and Zenitani2003; Xie & Watanabe, Reference Xie and Watanabe2007). 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., Reference Yasue, Harada and Takasuka2016). 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, Reference Yoshioka1988). 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., Reference Yasue, Utsumi and Takeda2006). 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., Reference Kinoshita, Yasuda, Watanabe, Watai, Imoto, Kamimura, Kono and Takahashi2022; Kono et al., Reference Kono, Takahashi, Yasuda, Watanabe, Watai, Imoto, Kinoshita and Nishijima2022). 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., Reference Yasue, Utsumi and Gosho2008), 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., Reference Kono, Takahashi, Yasuda, Watanabe, Watai, Imoto, Kinoshita and Nishijima2022). 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., Reference Egan, Bloom, Kuo, Hammer, Tongnunui, Iglésias, Sheaves, Grudpan and Simons2018; Caballero-Huertas et al., Reference Caballero-Huertas, Frigola-Tepe, Coll, Muñoz and Viñas2022). In species of the genus Engraulis, Japanese anchovy has been discriminated from E. encrasicolus through comparison of Cyt b sequence data (Grant et al., Reference Grant, Leslie and Bowen2005; Nakategawa et al., Reference Nakategawa, Dairiki, Itoi, Sugita and Akimoto2009). Cyt b (Yu et al., Reference Yu, Kong, Guo, Jiang, Zhuang and Jin2005; Zheng et al., Reference Zheng, Zou and Han2015), CR (Liu et al., Reference Liu, Gao, Zhuang, Jin, Yokogawa and Zhang2006) and microsatellite (Yu et al., Reference Yu, Lee, Huang and Chiu2002) markers were used to examine population structure of Japanese anchovy. Cyt b (Silva et al., Reference Silva, Lima, Martel and Castilho2014b), CR (Viñas et al., Reference Viñas, Sanz, Peñarrubia, Araguas, García-Marín, Roldán and Pla2014; Cuttitta et al., Reference Cuttitta, Patti, Maggio, Quinci, Pappalardo, Ferrito, De Pinto, Torri, Falco, Nicosia, Musco, Armeri, Placenti, Tranchida, Mifsud, Bonanno and Mazzola2015) and microsatellite (Zarraonaindia et al., Reference Zarraonaindia, Pardo, Iriondo, Manzano and Estonba2009; Ruggeri et al., Reference Ruggeri, Splendiani, Occhipinti, Fioravanti, Santojanni, Leonori, De Felice, Arneri, Procaccini, Catanese, Tičina, Bonanno, Cerioni, Giovannotti, Grant and Barucchi2016; Ouazzani et al., Reference Ouazzani, Benazzou, Charouki, Bonhomme and Chlaida2017) markers as well as the combination of them (Borrell et al., Reference Borrell, Piñera, Sánchez Prado and Blanco2012; Oueslati et al., Reference Oueslati, Fadhlaoui-Zid, Kada, Augé, Quignard and Bonhomme2014; Silva et al., Reference Silva, Horne and Castilho2014a) 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, Reference Yasue and Takasuka2009) 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., Reference Yasue, Harada and Takasuka2016). 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.
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.
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, Reference Noichi and Okiyama2014). 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., Reference Yasue, Harada and Takasuka2016).
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., Reference Funamoto, Kikuchi, Aoki, Watabe and Taniuchi1999; Inoue et al., Reference Inoue, Miya, Tsukamoto and Nishida2000; Silva et al., Reference Silva, Horne and Castilho2014a; 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).
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., Reference Thompson, Higgins and Gibson1994; Kumar et al., Reference Kumar, Stecher, Li, Knyaz and Tamura2018). 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., Reference Inoue, Miya, Tsukamoto and Nishida2001), KF765500.2 (Zhang & Xian, Reference Zhang and Xian2015), and AP017957.1 (Lavoué et al., Reference Lavoué, Bertrand, Wang, Chen, Ho, Motomura, Hata, Sado and Miya2017). BLAST analysis was also conducted to confirm identities with the other species of the genus Engraulis: E. encrasicolus NC_009581.1 (Lavoué et al., Reference Lavoué, Miya, Saitoh, Ishiguro and Nishida2007) and E. mordax NC_041097.1 (Lewis & Lema, Reference Lewis and Lema2019). 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., Reference Rozas, Ferrer-Mata, Sánchez-DelBarrio, Guirao-Rico, Librado, Ramos-Onsins and Sánchez-Gracia2017) to examine genetic diversity. An analysis of molecular variance (AMOVA) considering the six samples as a single group and a pairwise FST analysis were conducted using the Arlequin ver 3.5.2.2 software (Excoffier & Lischer, Reference Excoffier and Lischer2010) to examine genetic differences among samples. The level of statistical significance was set at P < 0.05 after Bonferroni correction for a pairwise FST analysis. Genetic distances were calculated based on the Kimura-2 parameter (K2P) model (Kimura, Reference Kimura1980), and the unrooted neighbour-joining (NJ) tree was constructed using the MEGA X software (Kumar et al., Reference Kumar, Stecher, Li, Knyaz and Tamura2018) 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. (Reference Chiu, Lee, Huang and Yu2002) 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 FST 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 FST analysis. The GenALEx 6.5 software (Peakall & Smouse, Reference Peakall and Smouse2006, Reference Peakall and Smouse2012) 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.
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 2 = 0.842, P < 0.001); February (Sample codes O, U, and W): ln HL = –3.630 + 1.603 ln SL (N = 96, R 2 = 0.936, P < 0.001); April (Sample code X): ln HL = –3.630 + 1.587 ln SL (N = 32, R 2 = 0.973, P < 0.001); August (Sample code Y): ln HL = –3.331 + 1.555 ln SL (N = 31, R 2 = 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.
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.
Genetic statistics for the mitochondrial DNA of Japanese anchovy in the Kii Channel
N, number of individuals; nh, 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 FST analysis revealed that FST 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).
Results of an analysis of molecular variance based on mitochondrial DNA for Japanese anchovy in the Kii Channel
Pairwise FST and P values between samples based on mitochondrial DNA for Japanese anchovy in the Kii Channel
FST 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.
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.
χ 2 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 FST analysis revealed that FST 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).
Genetic statistics for five microsatellite loci of Japanese anchovy in the Kii Channel
N, number of individuals; NA, number of alleles; HO, observed heterozygosity; He, 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 FST analysis revealed that FST 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).
Results of an analysis of molecular variance based on five microsatellite loci for Japanese anchovy in the Kii Channel
* <0.05.
Pairwise FST and P values between samples based on five microsatellite loci for Japanese anchovy in the Kii Channel
FST 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., Reference Moreno, Claramunt and Castro2011). 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., Reference Chen, Tzeng and Chiu2010; Katamachi et al., Reference Katamachi, Ikeda and Uno2015). 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., Reference Cuttitta, Patti, Maggio, Quinci, Pappalardo, Ferrito, De Pinto, Torri, Falco, Nicosia, Musco, Armeri, Placenti, Tranchida, Mifsud, Bonanno and Mazzola2015). 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., Reference Chen, Tzeng and Chiu2010). 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., Reference Yasue, Harada and Takasuka2016) 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, Reference Takahashi and Watanabe2004). 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., Reference Caneco, Silva and Morais2004; Khan et al., Reference Khan, Bal, Battal and Seyhan2022). 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., Reference Borrell, Piñera, Sánchez Prado and Blanco2012) and CO (Yu et al., Reference Yu, Kong, Guo, Jiang, Zhuang and Jin2005) regions, single nucleotide polymorphisms (Zarraonaindia et al., Reference Zarraonaindia, Iriondo, Albaina, Pardo, Manzano, Grant, Irigoien and Estonba2012; Huret et al., Reference Huret, Lebigre, Iriondo, Montes and Estonba2020), nuclear DNA intron sequences (Borsa et al., Reference Borsa, Collet and Durand2004; Bouchenak-Khelladi et al., Reference Bouchenak-Khelladi, Durand, Magoulas and Borsa2008), allozyme analysis (Sanz et al., Reference Sanz, García-Marín, Viñas, Roldán and Pla2008) and restriction fragment length polymorphism analysis of mitochondrial DNA (Magoulas et al., Reference Magoulas, Castilho, Caetano, Marcato and Patarnello2006) 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., Reference Yu, Kong, Guo, Jiang, Zhuang and Jin2005) and 0.4706–0.9935 (Zheng et al., Reference Zheng, Zou and Han2015) in the Yellow Sea and East China Sea, and 0.995–1.00 (Chen et al., Reference Chen, Tzeng and Chiu2010) 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., Reference Liu, Gao, Zhuang, Jin, Yokogawa and Zhang2006). 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., Reference Yu, Lee, Huang and Chiu2002). 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 FST 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., Reference Nakategawa, Dairiki, Itoi, Sugita and Akimoto2009). In captive conditions, the size and nutrient content of eggs and batch intervals were highly constrained by temperature (Yoneda et al., Reference Yoneda, Kitano, Tanaka, Kawamura, Selvaraj, Ohshimo, Matsuyama and Shimizu2014). 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, Reference Takahashi and Watanabe2004). 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., Reference Cuttitta, Patti, Maggio, Quinci, Pappalardo, Ferrito, De Pinto, Torri, Falco, Nicosia, Musco, Armeri, Placenti, Tranchida, Mifsud, Bonanno and Mazzola2015). In the Kii Channel, there is substantial seasonal variability in environmental conditions such as sea temperature and plankton density (Yasue & Takasuka, Reference Yasue and Takasuka2009). 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., Reference Yasue, Harada and Takasuka2016). 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, Reference Yasue and Takasuka2009). 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, Reference Denton and Nicol1965) 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, Reference Kristoffersen and Magoulas2008).
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, Reference Grant and Bowen1998). Although environmental conditions are different between OA and IA (Yoshioka, Reference Yoshioka1988), 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., Reference Funamoto, Kikuchi, Aoki, Watabe and Taniuchi1999; Liu et al., Reference Liu, Gao, Zhuang, Jin, Yokogawa and Zhang2006; Zheng et al., Reference Zheng, Zou and Han2015): 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., Reference Yu, Lee, Huang and Chiu2002; Chen et al., Reference Chen, Tzeng and Chiu2010), but no genetic differences have been observed in other areas, including coastal and offshore areas (Funamoto et al., Reference Funamoto, Kikuchi, Aoki, Watabe and Taniuchi1999; Liu et al., Reference Liu, Gao, Zhuang, Jin, Yokogawa and Zhang2006; Zheng et al., Reference Zheng, Zou and Han2015). 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., Reference Cuttitta, Patti, Maggio, Quinci, Pappalardo, Ferrito, De Pinto, Torri, Falco, Nicosia, Musco, Armeri, Placenti, Tranchida, Mifsud, Bonanno and Mazzola2015; Huret et al., Reference Huret, Lebigre, Iriondo, Montes and Estonba2020), the present study provides useful information on the population structure for defining stocks of Japanese anchovy in Japan.
Acknowledgements
The authors are grateful to Shunsuke Ito (Bioinsight Co., Ltd) for his laboratory assistance. The sea temperature data in OA were supplied by Dr Yasuyuki Baba (Kyoto University). The present study benefited from discussions with Drs Shigeo Harada and Ryu Doiuchi (Wakayama Prefectural Fisheries Experimental Station).
Authors’ contributions
N. Y. designed the study and handled the sample collection. N. Y. and D. N. performed the laboratory work. N. Y., T. Y. and D. N. conducted data analysis. N. Y. and A. T. drafted the initial manuscript. All the authors contributed to elaborating the study processes, interpreting the results, and developing the conclusions.
Financial support
The present study was financially supported by the research and assessment programme for fisheries resources from the Fisheries Agency of Japan and by Challenging Research (Pioneering) (KAKENHI No. 20K20455) from the Japan Society for the Promotion of Science (JSPS).
Competing interest
None.
Ethical Standards
The authors used Japanese anchovy landed by the commercial fishery in the Kii Channel. The present study did not include any animal experimentation or harm any organisms, and thus formal ethical approval is not required.
Data availability
Cyt b and CR sequences are available in DDBJ under the accession numbers LC802150 to LC802632.
