Pseudodiaptomus marinus Sato, 1913 in the Black Sea: morphology, genetic analysis, and variability in seasonal and interannual abundance

Sample collection and zooplankton processing

Sevastopol Bay is located in the northern part of the Black Sea at the southwestern tip of the Crimean Peninsula (Fig. 1). It is a semi-enclosed estuarine-type bay with restricted water exchange, is about 7 km long and 1 km wide at its widest point, and has an average depth of 12 m. According to long-term hydrological observations, Sevastopol Bay’s average sea surface temperature (SST) varies between 6 and 8 °C during winter and between 23 and 26 °C during summer (Garmashov, 2019). During the studied period, SST values stayed within this range of average values. The salinity values did not significantly change and were within the range of 17–18 ppt. The bay is a port area greatly affected by maritime traffic and anthropogenic factors, including NIS invasions. Regular investigations of the bay’s zooplankton were conducted in 1976, 1980, 1989, 1990, 1995, and 1996. They resumed in 2002 and have continued to this day.

Zooplankton samples were collected between 2016 and 2018 during routine surveys at three stations located within or adjacent to the Sevastopol Bay area (Fig. 1). In the bay, samples from Stations 2 and 3 (St. 2 and 3) were usually taken twice per month. Gaps occurred in May, July and October 2016, January and February 2017, and January, May, and September 2018 due to technical malfunction or meteorological conditions. Outside Sevastopol Bay, Station 1 (St. 1) samples were collected less frequently: in July, August, October, and December 2018. We collected samples using a Juday plankton net (with a mouth area of 0.1 m² and a mesh size of 150 µm) from the whole water column: 40–0 m from St. 1, 10–0 m from St. 2, and 9–0 m at St. 3. Samples were collected in the morning and fixed with formaldehyde solution (4% final conc). They were processed in the laboratory using the methodology for zooplankton (Postel, Fock & Hagen, 2000; Aleksandrov et al., 2014). The P. marinus nauplii and copepodites were identified according to Uye & Onbe (1975).

We preserved the additional sample taken from St. 3 in December 2018 in 95% ethanol immediately after collection and stored it at 4 °C. Fifteen P. marinus individuals, CIV–CV, were extracted from the sample for genetic analysis.

Morphological studies

We selected specimens for morphological study from zooplankton samples collected at Stations 2 (two females) and 3 (three females and three males) on October 8, 2018. The formalin-fixed specimens were immersed in a 1:1 solution of glycerin and water on glass slides, and were then measured and dissected. All operations were performed using a LOMO MBR-10 stereomicroscope. All line drawings were made using glycerine-mounted specimens and a camera lucida on a Leica DM LS2 compound microscope at magnifications 200 ×, 400 ×, and 1, 000 ×. We stained specimens with a solution of 1% chlorazol black E (SBE) dissolved in 70% ethanol for better visibility of the setae and aesthetascs.

Abbreviations used in the P. marinus drawings and descriptions were: Pr, prosome; Ur, urosome; Ur1–5, urosome somites; A1, first antenna (antennula); P5, fifth leg; Enp, endopod; Exp, exopod; and Exp1–3, exopod segments.

Molecular genetic analysis

We obtained nucleotide sequences of the cytb gene fragment and nuclear ribosomal internal transcribed spacer (ITS1, 5.8 S, ITS2) using PCR-based sequencing procedures for P. marinus specimens.

DNA extraction. DNA was extracted using a DNK-EHKSTRAN Kit (Syntol, Moscow, Russia). We selected single specimens from the ethanol-preserved samples and incubated them overnight in 100 µL of lysis buffer (Syntol) with 5 µL of Syntol Proteinase K and 1 µL of 2-mercaptoethanol at 56 °C. After lysing, the specimens were crushed with a vortex pestle for 20 s and we carried out DNA extraction according to the DNK-EHKSTRAN Kit manufacturer’s protocol. The elution volume was 30 µL, and the DNA was stored at −20 °C.

PCR amplification and sequencing. The total 20 µL reaction mix volume for both DNA fragments was prepared as follows: 5xPCR Mix with MgCl₂ (Evrogen, Moscow, Russia), 0.6 µM of each primer, and 2 µL of template DNA.

The cytb fragment was amplified using the primers UCYTB151F (5′-TGT GGR GCN ACY GTW ATY ACT AA-3′) and UCYTB270R (5′-AAN AGG AAR TAY CAY TCN GGY TG-3′) (Merritt et al., 1998). The reaction conditions included initial denaturation at 94 °C for 2 min, 35 cycles of denaturation at 94 °C for 30 s, annealing at 50 °C for 30 s, extension at 72 °C for 1 min, and a final extension at 72 °C for 4 min (Ohtsuka et al., 2018).

The ITS2 fragment was amplified using the primers ITS4 (5′-TCCTCCGCTTATTGATATGC-3′) and ITS1 (5′-TCCGTAGGTGAACCTGCGG-3′) (White et al., 1990). The reaction conditions included initial denaturation at 95 °C for 3 min, 40 cycles of denaturation at 95 °C for 30 s, annealing at 45 °C for 30 s, an extension at 72 °C for 30 s, and a final extension at 72 °C for 5 min.

We electrophoretically separated amplicons using a 1% agarose/TBE buffer gel with ethidium bromide, and visualised them using an ultraviolet (UV) transilluminator. PCR products were sequenced in both directions using the standard BigDye Terminator Cycle Sequencing Ready Reaction Kit (Applied Biosystems Inc., Foster City, CA, USA) and an ABI PRISM 3500xL analyser (Evrogen, Moscow, Russia).

Sequence analysis. We aligned the rDNA fragments using the BioEdit software program (Hall, 1999) and a P. marinus speciment from GenBank as a reference sequence for ITS2 (GenBank KT808252), and then we manually refined the alignment. The multiple alignment was run using ClustalW (Thompson, Higgins & Gibson, 1994) in MEGA7 (Kumar, Stecher & Tamura, 2016) software. We deposited all nucleotide sequences generated during this study in the GenBank database.

Evolutionary analysis was conducted using MEGA7. As shown earlier, ITS1 is more variable than ITS2 (Sabia et al., 2017) and is recommended for inter-population surveys, while ITS2 is usually used more for species discrimination. We decided to conduct phylogenetic analysis for both rRNA regions of ITS1 and ITS2. We determined the best substitution model using the Bayesian information criterion (BIC) and Akaike information criterion (AIC) model tests for best fit (Posada & Crandall, 2001; Raftery, 1995). The Kimura-2-parameter model and a discrete Gamma distribution (five categories, +G) were used for phylogenetic reconstruction using the maximum likelihood (ML) method (Kimura, 1980). We also determined the evolutionary distances inside the Pseudodiaptomus genus using the Kimura 2-parameter model (Kimura, 1980). The rate variation across sites was modeled with a Gamma distribution (shape parameter = 1). The analysis involved 10 Pseudodiaptomus species nucleotide sequences, and all positions with less than 95% site coverage were eliminated.

We aligned cytb sequences from our sample against sequences found in GenBank using BLAST (Altschul et al., 1990).

Morphological examination

Description of female. Total body length 1.25–1.29 mm (mean 1.27 mm; n = 5).

Prosome (Figs. 2A, 2B) elongated, about 2.5 times longer than wide, with broadly rounded anterior part in lateral view (Fig. 2A) covered with rare short hairs (not shown in Fig. 2B). Pedigerous somites 4 and 5 fused. Posterior corners produced into sharp spines directed outward. Rostrum with paired filaments (Fig. 2C). Eyes clearly visible, both in lateral and dorsal projections.

Urosome four-segmented (Figs. 2A, 2B, 2D). Pr/Ur about 1.8. Genital double-somite (Ur1) similar to that shown in a scanning electron micrograph by Soh et al. (2001: p. 213, Fig. 8B): asymmetrical, strongly swollen ventral side (Figs. 2A, 2D), lateral surfaces with asymmetrical patches of fine spinules, rows of stronger spinules on left lateral surface and proximal part of ventral protrusion, and one pair of thin setae near lower edge of ventral protrusion (Figs. 2E, 2F). Ur2 with rather long hairs perpendicular to right lateral surface (Figs. 2E, 2F). Ur3 longer than Ur2 and Ur4. Postero-dorsal margins of Ur1–Ur3 with rows of flat triangular teeth increasing in size from first to third somites (Figs. 2A, 2B, 2D).

Caudal rami symmetrical and diverge (Fig. 2B). Each ramus about three times longer than wide, with plumose inner side, and six setae (one lateral, one dorsal and four distal). All caudal setae plumose, except for dorsal seta with rare hairs. Longest caudal seta about 16% female total length.

Antennulae symmetrical, 22-segmented (segments 6 and 7 partly fused; counted separately), reaches about mid-length of Ur2 (Fig. 2A). Number and location of setae, spines, and aesthetascs on each segment shown in Fig. 2G and Table 1. First segment has groups of small spines and long hairs on ventral surface (Fig. 2H). Segments 1–4 bear tiny hair-like seta (Fig. 2G; not mentioned in Table 1). Segment 20 has two distal modified setae; shorter of them has thickened and elongated plumose base (Fig. 2G, arrow 1), and second seta with several pairs of recurved spines (Fig. 2G, arrow 2).

Fifth legs symmetrical, uniramous (Figs. 2I, 2J). Coxa with row of spines on anterior mid-surface (Fig. 2I) and tiny spines in posterior view (Fig. 2J). Basis with medial plumose seta on posterior surface and one or two small spinules on distal outer corner. Exp three-segmented. Exp1 elongated, with one distal outer spine. Exp2 with one outer spine and well-developed pointed serrated process on inner side, which about equal in length to Exp2. Exp3 pointed, slightly longer than serrated process of Exp2, with small inner spiniform process near base.

Description of male. Total body length 1.05–1.11 mm (mean 1.08 mm; n = 3).

Prosome (Figs. 3A, 3B) similar to that of female, but slightly more slender, about 2.6 times longer than wide. Posterior corner spines directed backwards. Rostrum with paired filaments (Fig. 3C). Eyes located close to each other (Fig. 3A).

Urosome five-segmented (Figs. 3A, 3B). Pr/Ur about 2.0. Distal margins of Ur1–Ur4 with rows of flat triangular teeth similar to those in female, but around the entire perimeter of somite, excluding Ur1 (teeth only on dorsal margin). Size of teeth increases from Ur1 to Ur4. Ur2 with two paired groups of small spines on ventral surface (Fig. 3D).

Caudal rami (Fig. 3A) similar to those of female, but most caudal setae relatively longer. Longest seta about 25% male total length.

Antennulae asymmetrical. Left A1 similar to that of female (not shown in Fig. 3). Right A1 geniculate (between segments 18 and 19), 21-segmented (segments 6–7 and 10–12 appear to be partly fused; counted separately) (Figs. 3E, 3F; Table 1). Ventral surface of first segment has no small spines or long hairs (Fig. 3E). Segments 1–3 have tiny hair-like seta (not mentioned in Table 1). Segment 10 has a large hooked spine. Segment 20 with two pairs of setae, distal and medial; one seta from each pair strongly shortened, and long medial seta (Fig. 3E, arrow 1) similar to shorter modified seta on segment 20 of female A1 (see Fig. 2G, arrow 1).

Fifth legs (Figs. 3G, 3H) strongly asymmetrical with well-developed endopods on both left and right P5. Coxa similar to that of female, with row of spines on anterior surface (Fig. 3G) and tiny spines in posterior view (Fig. 3H). Right basis larger than left, with row of thin spinules on proximal margin in anterior view. Left basis with fine proximal hairs in anterior view, and with short pointed protrusion proximally on inner lateral surface. Both right and left basal segments with medial plumose setae on posterior sides and with rows of spines on outer lateral surfaces.

Left Exp (Fig. 3G) bi-segmented. Exp1 short, with one strong spine on outer distal corner. Exp2 elongated, narrowed distally, with three thin spinules on distal half of inner margin, and two spines on outer margin (lateral spine much larger than distal one); spinulous margin between these spines.

Left End uni-segmented (Fig. 3G), elongated, lob-shaped, narrowed distally, reaches more than mid-length of Exp2.

Right Exp (Figs. 3G, 3H) three-segmented. Exp1 has stout distolateral spine forked near mid-length with small process between rami; outer ramus longer and thinner than inner. In posterior view (Fig. 3H), Exp1 has short thick spine near base of forked spine, small medial seta, and patch of thin hairs. Exp2 elongated, with long, straight, partly serrated spine on outer distal corner and small medial seta on posterior surface. Exp3 has long sickle-shaped outer process, short inner process, and two small setae near their bases.

Right Enp (Figs. 3G, 3H) has two rami extending beyond distal edge of right basis. One of them bifurcated at apex, slender, and slightly longer than another ramus. Thicker ramus ends with four to five spiniform processes (Figs. 3G, 3I), two of which sharp-pointed and others with slightly rounded tips. Two pointed processes may seem fused into one depending on the angle of view.

Molecular genetic analysis

The rRNA sequences, including the ITS1, 5.8S, and ITS2 regions obtained from most of the individuals collected from the Black Sea, revealed a notable intraindividual heterogeneity observed using the electropherograms obtained by direct and reverse sequencing. We detected the presence of several allelic variants by double peaks. Because divergent rRNA copies affect phylogenetic reconstruction (Karlep, Reintamm & Kelve, 2013), we excluded all heterogenetic sequences from the analysis. Only three sequences had one allelic variant, and these sequences were deposited in GenBank (accession numbers MN555816–MN555818). The length of these sequences varied: PMB7 –  706 bp (ITS1, 5.8S, ITS2), PMB11 –  416 bp (ITS2), and PMB14 –  414 bp (ITS2).

During our calculations, we treated the beginning of the rRNA sequences (which did not have information about ITS1) as missing data. The total length of alignment was 742 bp. We conducted the rDNA-based phylogeny with some Pseudodiaptomidae and Diaptomidae species sequences from GenBank (their accession numbers are presented in Fig. 4), and we used Pseudodiaptomus specimens from the Black Sea and cyclopoid copepod Oithona similis (KF153700) as the outgroups.

All Pseudodiaptomus individuals collected from the Black Sea were grouped with two other P. marinus individuals in a separate clade (ML, 91%) found in Genbank. We identified relatively high genetic heterogeneity within this species. The phylogenetic tree demonstrated 3% divergence between individuals from South Korea and individuals from the Black Sea and Lake Faro (in the Mediterranean basin). The genetic distance within the P. marinus species was the same as the distance between different Pseudodiaptomus species (Table 2, shown in bold).

We obtained fifteen sequenced fragments for the mitochondrial cytb gene with lengths ranging from 344 to 374 bp and deposited them in the GenBank (accession numbers MN561264–MN561278). Based on the blast similarities in GenBank, we identified all haplotypes as P. marinus (Ohtsuka et al., 2018). Three haplotypes were found in a total of 15 cytb sequences from the Black Sea: 10 individuals demonstrated 99.73% identity (100% query) with AB920868, three individuals demonstrated 99.73% identity (100% query) with AB920869, and we found 99.73% identity (100% query) with AB920872 in two individuals. These sequences corresponded to the H1, H2, and H5 haplotypes (Ohtsuka et al., 2018). All haplotypes were shared between Japanese coastal waters and the San Francisco Estuary, and only two (H1 and H2) were detected from Haeui Island, Korea.

Seasonal and interannual variations in abundance

P. marinus was first found in the Black Sea during a routine survey of Sevastopol Bay at the end of September 2016. We recorded 23 ind m⁻³ of this alien species at each developmental stage, including six females at the centre of the bay (St. 3) and only two nauplii at its mouth (St. 2). In November 2016, the abundance of P. marinus had increased significantly, reaching up to 1,373 ind m⁻³ at St. 3 and 103 ind m⁻³ at St. 2 (Fig. 5), which were the highest densities that year. The following month, the population abundance had decreased sharply, and only one female was found at St. 3.

In 2017, P. marinus was not found until September, except for one young copepodite found in June at St. 3 (Fig. 5). The abundance was low (6 ind m⁻³) in September and increased to 73 ind m⁻³ in October, which was significantly lower than the density in October 2016, which saw the highest recorded abundance. Therefore, we collected samples every week in November 2017 to determine the annual population maximum. However, the abundance varied between 12 and 46 ind m⁻³ and did not exceed the levels measured in October. In December, the abundance decreased to 9 ind m⁻³.

In 2018, P. marinus was not found in samples during the first half of the year, similarly to 2017. Single copepodites appeared in June and were found in July and early August at both stations. We observed a small rise in abundance at the end of August (33 ind m⁻³) at St. 3. The peak in abundance was recorded in early October when the density increased to 5,184 ind m⁻³ at St. 3 and 2,157 ind m⁻³ at St. 2, which was the highest P. marinus abundance during the whole research period. Copepodites, nauplii, and adults of both sexes, including egg-carrying females, were found in this sample. Additionally, males were recorded for the first time in the studied period. In November and December, the abundance dropped down to 21 ind m⁻³ and 197 ind m⁻³ at St. 3, and 8 ind m⁻³ and 39 ind m⁻³, at St. 2 , respectively.

Altogether, we found P. marinus in 12 of the 26 samples collected from the centre of Sevastopol Bay (St. 2), and in 19 of the 26 samples collected at its mouth (St. 3). Outside the bay (St. 1), P. marinus was not found during the entire studied period.

Our year-round observations of P. marinus in Sevastopol Bay during 2017–2018 showed that it did not appear during the first half of the year in the pelagic community (Fig. 5). Single young copepodites and nauplii appeared in June when the temperature reached 21−22 °C. We recorded considerable levels from September to the end of the year, with the maximum abundance recorded in October. The water temperature was 17.3 °C in 2017 and 19.6 °C in 2018, which were lower than the optimum temperature (20−25 °C) for this species’ egg production (Uye & Kasahara, 1983). However, favourable temperatures (20−22 °C) were recorded two weeks earlier (from the 20th of September to the 6th of October) that could support intensive population growth. The development of P. marinus from eggs to adults was shown to be 13 days (Huang, Zhu & Liu, 2006) so the species abundance can rapidly increase under the favourable conditions (Huang, Zhu & Liu, 2006; Deschutter et al., 2018).

Taxonomic identification and variability; vector of introduction

P. marinus belongs to the Ramosus species group and hickmani-subgroup, which include nine morphologically close species (Walter, Ohtsuka & Castillo, 2006). P. marinus is distinct from these species because of its combination of the following features (based mostly on the structural details of male P5):

(a) distolateral spine on right Exp1 bifurcated near mid-length, with the outer ramus longer and thinner than the inner ramus (typical only for P. marinus);

(b) left Enp narrowed distally (not rounded);

(c) well-developed right Enp with two rami extending beyond the distal edge of right basis, thicker ramus ends with three to six spiniform processes.

The drawings and descriptions of P. marinus can be found in Sato (1913: pl. 7, Figs. 69–71), Brodskii (1948: p. 63–64, 96, 116, Table 18, Figs. 1–7), Brodskii (1950: p. 322, Fig. 225), Tanaka (1966: p. 44, Fig. 4 (from Razouls et al., 2005–2020)), Kos (1985: p. 236, Figs. 6, 7), Nishida (1985: p. 132, Fig. 4f–j), Walter (1986: p. 146–147, Fig. 7L–O), Fleminger & Kramer (1988: p. 539–540, Fig. 4), Soh et al. (2001: p. 213, Fig. 8B), and Brylinski et al. (2012: p. 580, Fig. 3). In contrast, specimens from the coastal waters of the Indian Ocean (Grindley & Grice, 1969; Pillai, 1976) and eastern Australia (Greenwood, 1977) are not P. marinus, but distinct and closely related species (Walter, 1986; Walter, 1987; Fleminger & Kramer, 1988).

After comparing diagnostic features, we observed that Pseudodiaptomus specimens from the Black Sea resembled P. marinus Sato, 1913. Some differences, namely the perpendicular hairs on the right lateral surface of Ur2, a group of long hairs on the ventral surface of the first A1 segment in females, and a modified seta with a thickened plumose base on segment 20 of the right A1 in males, are probably not specific to Black Sea specimens but have not been previously recorded. However, we did note the variability in the number of spiniform processes at the end of the P5 right Enp in Black Sea males. Of the three examined males, two specimens had four processes (see Figs. 3G, 3H) and one specimen had five processes (see Fig. 3I). Walter (1986) re-examined individuals from various locations across Japan (including Nemuro Bay, Hokkaido, near the type locality of P. marinus) and from Ala Wai Canal, Hawaii, and found connections between the number of these processes and the geographical location of the re-examined individuals. Fleminger & Kramer (1988) found a similar variation to Walter’s (1986), but in individuals from the same bay (Mission Bay, Southern California). They concluded that “these variations appear to vary randomly and may not characterise geographically different stocks within the species” (Fleminger & Kramer, 1988: p 539). Our results were similar to those of Fleminger & Kramer (1988). However, further morphological studies are required for a clearer understanding of this phenomenon and the variability of other distinctive P. marinus features.

Currently, molecular genetic analysis is one of the most efficient ways to identify and confirm new invasive species, their region of origin, and possible dispersal pathways. Using the genetic variability of ITS2 and cytb, we determined that all analysed copepods from Sevastopol Bay belonged to the P. marinus species.

Moreover, our results revealed intraspecific genetic differentiation for P. marinus. Sabia et al. (2017) recorded the first evidence of such divergence and suggested that a ‘new form’ had inhabited Lake Faro, but no related studies followed. Our ribosomal ITS2 sequences showed small differences between Lake Faro and Black Sea specimens (a genetic distance of about 1.1% only), which allowed us to consider that these samples belonged to one group. At the same time, the ITS2 sequence from South Korea differed by 3.2–4.3%. The nucleotide diversity within P. marinus from South Korea is incompletely known due to sampling limitations. In copepods, ITS2 is considered less variable than ITS1 or COI (Makino, Knox & Duggan, 2010; Sabia et al., 2017) and the degree of interspecific differences varies from <1% to 55% for this marker. For the Eucalanidae family, no variation has been observed within species and only 0.2–3.4% divergence has been detected between the closest species (Goetze, 2003). We observed the same situation for the Diaptomidae family, in which phylogenetic reconstruction revealed small genetic interspecies differentiation (Fig. 4). Our interspecific genetic distance analysis for the Pseudodiaptomus genus showed an average divergence levels as high as 55% (Table 2). However, the P. japonicus, P. inopinus, and P. poplesia species group differed by 2.1–3.6%. (Table 2, shown in bold), which was even lower than the difference across P. marinus geographical groups. We confirmed the genetic differences for these three species using ITS1 (12–14%) and COI (17.6–26.7%) (Eyun et al., 2007; Soh et al., 2012). Therefore, closely related or recently diverged species may have quite low genetic differentiation that ranges between 2% and 5% at ITS2, making it likely that we found two cryptic species with 3.2–4.3% divergence (Table 2, shown in bold red).

When observing the ITS2 nucleotide diversity, we expected that individuals from Lake Faro and the Black Sea formed one group (the Mediterranean group) that differed significantly from the South Korean specimens. Low genetic distance values may be associated with the resolution of the molecular marker discussed above. We could not deduce the origin of the Mediterranean group from our genetic data because of insufficient information about ITS2 genetic diversity. Therefore, an increased number of samples from different parts of the distribution range and more variable molecular markers are required.

Previous studies have analysed the genetic and haplotype diversity of P. marinus populations from the coastal regions of Japan, Korea, and the San Francisco Estuary based on mtDNA cytb-sequenced data (Ohtsuka et al., 2018). Of the 39 haplotypes detected, three (H1, H2, and H5) accounted for more than half of the population in six of the seven studied localities. The only exception was Haeui Island, Korea, where 11 haplotypes were apparently endemic and only two (H1 and H2) were shared between the Korean, Japanese, and San Francisco Estuary localities. However, these haplotypes had a low frequency in Haeui Island that was different than other areas where they made up a significant part of the P. marinus population. Thus cytb sequence investigations, like the ITS2 ones sited above, reveal the specificity of the P. marinus population from Korean region. This fact allow to assume that another form (or cryptic species) is located just in South Korean waters not in Lake Faro as Sabia and coauthors have supposed (Sabia et al., 2017) and the Mediterranean P. marinus groups thus may be close to the population from other regions of the Western Pacific.

The same ubiquitous haplotypes H1, H2, and H5 (Ohtsuka et al., 2018) were found in the Black Sea at frequencies of approximately 67, 13, and 20%, respectively. The Black Sea had the lowest number of haplotypes of the studied localities (Ohtsuka et al., 2018). The low haplotype diversity of Sevastopol Bay P. marinus was consistent with the general pattern of genetic depletion in a newly introduced population, which is usually associated with strong genetic drift or bottleneck effects (Geburzi & McCarthy, 2018).

Using mtDNA cytb diversity analysis, Ohtsuka et al. (2018) confirmed that there were multiple P. marinus introductions in the San Francisco Estuary in the 1990s from several Japanese localities, although introductions from Korea could not be excluded. P. marinus was most likely introduced into the Black Sea from the Mediterranean Sea, where it had widely spread during the previous decade, and from where vessels had most often came into Sevastopol Bay. However, mtDNA cytb sequences for P. marinus populations in Mediterranean or South Atlantic waters cannot be found in the NCBI. It is currently impossible to conduct reliable genetic analysis to determine the donor region of the P. marinus population in the Black Sea.

Ballast water from trans-oceanic ships have been suggested as the main vector of P. marinus introduction into European coastal marine environments (Lučić et al., 2015; Sabia et al., 2015; Deschutter et al., 2018; Uttieri et al., 2020). Another suggested distribution route for this copepod is aquaculture (Zagami et al., 2018; Sabia et al., 2014), although there are no aquaculture farms inside Sevastopol Bay. It is most likely that P. marinus was introduced into the Black Sea via ballast water from ships, similarly to other non-indigenous copepods Acartia tonsa and O. davisae (Gubanova, 2000; Altukhov, Gubanova & Mukhanov, 2014), since P. marinus was first found in Sevastopol Bay and did not spread elsewhere until at least the end of 2018. The previous invader, O. davisae, was also initially found within the bay and began to spread to other coastal areas of the Black Sea four years later.

Main causes of successful spread and establishment

P. marinus has been found in more than 10 countries (Uttieri et al., 2020; Suzuki, Nakayama & Tanaka, 2013). Great environmental adaptability and some biological features have allowed this species to survive in new areas. This copepod thrives in habitats with temperatures ranging from 6.3 to 31.5 °C and 2.5 to 38.5 ppt salinity (Sabia et al., 2015; Uttieri et al., 2020). According to an experimental study by Svetlichny et al. (2019), the adult P. marinus has a salinity tolerance ranging between 5.0 and 44.0 ppt, and a temperature tolerance ranging from 8.0 °C and 27.0 °C. P. marinus is epibenthic during the day, but swims up the water column at dusk and acts as a pelagic species during the night (Uye & Kasahara, 1983; Sabia et al., 2015). Therefore, both the pelagic and benthic feeding modes are available to P. marinus, and copepod can act as either herbivores or detritivores, respectively (Uye & Kasahara, 1983; Sabia et al., 2015). The species’ ecology and biological plasticity allow it to adapt to the Black Sea’s low salinity and low winter water temperatures and successfully compete (or avoid competition) with native species. Between 2016 and 2018, the salinity in Sevastopol Bay ranged from 17 to 18 ppt, and the temperature ranged from 6.2 to 26.2 °C. The lowest temperature was reported in December 2016 and the highest in August 2018.

Normally, the invasion of a new species is preceded by a series of changes in an ecosystem (Alimov et al., 2004). The imbalance of the Black Sea ecosystem was initially caused by eutrophication, pollution, and overfishing in the 1970s and 1980s. In turn, this preconditioned the invasion of predatory ctenophores in the 1990 and 2000s, which resulted in profound changes in the basin ecosystem (Shiganova et al., 2001; Gubanova et al., 2001; Kideys, 2002; Finenko & Romanova, 2000; Oguz, Fach & Salihoglu, 2008; Gubanova et al., 2019).

Seasonal dynamics in the new environment

Previous studies on P. marinus’ seasonal and interannual variability in abundance have shown that, since its first appearance in 2016, P. marinus has successfully survived in the new environment of Sevastopol Bay. We collected this species at all stages of its development, including ovigerous females and nauplii within the studied period of 2016 to 2018 (Fig. 5). Therefore, we concluded that the new invasive species reproduced there and established a self-sustaining population. We also observed identical seasonal patterns with one pronounced abundance peak in the autumn during 2017 and 2018. However, the abundance was noticeably higher in 2018 where we observed the species in plankton from July to December, while in 2017, the species was not found until September, with the exception of one individual in June, Our findings most concern St. 3, which was located in the centre of the bay. In the mouth of the bay at St. 2, P. marinus appeared more rarely, particularly during 2017, and with a lower density in comparison to St. 3. However, the frequency of P. marinus occurrence and abundance increased in 2018 at both stations in the bay. At the same time, we did not find specimens outside the bay throughout the studied period, although it is possible that P. marinus may spread outside the bay in the following year, similarly to O. davisae, another non-indigenous copepod (Altukhov, Gubanova & Mukhanov, 2014).

The density of P. marinus was mostly insignificant during the observed years. There were exceptions during the autumn peaks when its abundance reached 1,373 ind m⁻³ (in 2016), 73 ind m⁻³ (in 2017), and 5,183 ind m⁻³ in 2018 (the highest one). O. davisae was a dominating species in the copepod community, representing 70–90% of total copepod abundance along with Paracalanus parvus, Acartia clausi, and A. tonsa (Fig. 6). This has been the typical composition of copepod species for autumn in Sevastopol Bay since 2006 (Gubanova et al., 2019). The new member of the community, P. marinus, did not significantly affect the species ratio and made up less than 1% of the copepod community, except during its maximum abundance in autumn 2016 and 2018 when its percentage amounted to 9 and 11%, respectively, and exceeded that of the Acartia species (Fig. 6).

P. marinus is reported as a rare species at most sites in the Mediterranean Sea. Its abundance is usually lower than 200 ind m⁻³ and can even be lower than 20 ind m⁻³ (Uttieri et al., 2020). Sevastopol Bay’s abundance of more than 5,000 ind m⁻³ is one of the highest in the Mediterranean basin. However, according to our data, P. marinus seasonal dynamics are characterised by one very sharp increase in abundance that lasts for a short period of time. Regular and frequent sampling should be conducted to catch such short-term increase in the species’ number.

O. davisae has occupied a free ecological niche and has replaced O. nana in the Black Sea, while A. tonsa has ousted Acartia latisetosa due to a number of competitive advantages (Gubanova, 2000; Gubanova et al., 2014). Since P. marinus is not as abundant yet, its total effect on the native copepod community has not been revealed. P. marinus’ ability to switch between pelagic (night) and demersal (day) life modes may provide it with competitive advantages in the Black Sea community. However, this hypothesis has to be thoroughly tested in future investigations.