Coexistence of Two Copepods, Recorded for the First Time, in NW Iberian Shelf: The Case of Oithona atlantica and the Allochthonous O. davisae

: This paper compiles the data regarding the first occurrence of Oithona davisae and O. at-lantica in NW Spain, which is supported by morphological and molecular analysis. Additionally, we investigated the seasonal dynamics of the invasive O. davisae , revealing that its abundance is conditioned by upwelling-downwelling patterns in the R í as Baixas of Galicia. Temperature was the most correlated factor, with higher abundances in upwelling relaxation-downwelling events. More studies in long-term zooplankton dynamics and molecular analysis are needed to determine if O. davisae is displacing other native species of the same genus, such as O. atlantica , in Galician waters.


Introduction
Copepods are crucial components in marine planktonic food webs, functioning as the main transfer node between primary producers and fish, and represent a crucial role in nutrient recycling and export [1,2].Due to their short generation times, copepods can respond rapidly to environmental changes and have been identified as indicators of climatic variability [3][4][5].
The Galician coastal and shelf waters are highly productive ecosystems that support important commercial activities related to marine resources exploitation and aquaculture [6].The particular hydrodynamic conditions in the estuary of the Ría de Arosa favor successful shellfish production, mainly represented by filtering mollusks that feed on the planktonic community [7,8].Seasonal zooplankton dynamics in Galician waters have been described on a monthly basis by a zooplankton monitoring program (RADIALES), carried out by the Spanish Institute of Oceanography (IEO) in the nearshore waters of A Coruña and the proximity of the Ria of Vigo [9,10].However, field works in the Ria de Arosa on the planktonic assemblages are scarce and mainly focused on the phytoplankton community [8,[11][12][13].In this area, the zooplankton taxonomic composition has been described only in 1978 [14], and more recent studies have been directed to Bivalvia veliger [15] and Cirripedia larvae [16][17][18].
Three species of the genus Oithona are considered native to the NW Iberian Peninsula: O. plumifera, O. nana, and O. similis [19].All were recorded in the monitoring program RA-DIALES.Previously, Corral and Álvarez-Ossorio [14] described O. nana and O. helgolandica (unaccepted status, now accepted name O. similis) in the Ria of Arosa.
Oithona atlantica (Farran, 1908) (Copepoda, Cyclopoida) is considered an epimesopelagic and cosmopolitan species, present in the entire Atlantic Ocean from the Arctic to the Antarctic, including the Mediterranean Sea [20].Based on its morphological characteristics, it is closely related to O. longispina and often confused with O. plumifera.
Oithona davisae (Ferrari F.D. and Orsi, 1984) is known to be indigenous to the coastal waters of Japan and the China Sea [21].This species, sometimes initially misidentified as O. brevicornis, has colonized several American and European estuaries and bays in the last decades.O. davisae was originally described by Ferrari and Orsi [22] in Californian estuarine waters and also reported in the Bay of Varna [23], the Wadden Sea [24], Puget Sound [25], the Marmara Sea and the Golden Horn Estuary [26], and Venice Lagoon [27,28].This species is known as an invasive species to the Black [29] and Mediterranean Seas [30].
In Iberian Peninsula waters, Oithona davisae has been recorded in Barcelona Port, in the north-western Mediterranean [31], and in the Estuary of Bilbao in the Cantabric Sea [32].Currently, it is considered a non-indigenous species (NIS) and has been included in the European Alien Species Information Network (EASIN).Alien species can cause negative impacts on the ecosystems they colonize, displacing native species by increasing predation, competing with them for food and habitat resources [33] or introducing allochthonous pathogens in the colonized area.Increasing the knowledge of the zooplankton communities' taxonomic composition and their seasonal variance could be useful to address the early detection of NIS and evaluate their influence on the ecosystems.
To resolve the status of this copepod in Galician waters, we aimed to identify the Oithona species collected in the Ria of Arosa, applying morphological and molecular techniques, and to describe the seasonal variability of the allochthonous species O. davisae.

Sampling
Oblique tows were performed to study the presence and variability of zooplankton species close to shellfish beds in three areas of the Ría de Arosa (NW Spain): Lombos do Ulla, Corón, and Sarrido (Figure 1), covering depths from 4 to 12 m.A total of 13 surveys were conducted from November 2015 to March 2017.Due to some logistic problems, there were no samplings in December 2015, January 2016, March 2016, or February 2017.Fifteenminute oblique tows were carried out during the day using a 72 cm diameter double bongo plankton net of 200 µm mesh size and equipped with a flowmeter.The mesh size was chosen to capture mesozooplankton species that could potentially act as intermediate hosts for Bivalvia parasites.The mesh size used only collected adults of the genus Oithona.

Morphological Analyses
For microscopic examination and body measurements, Oithona spp.individuals preserved in formaldehyde were selected and put in a glycerin drop in depression slides.Individuals of O. davisae were identified following the descriptions of Ferrari and Orsi [22] and Temnykh and Nishida [21].O. atlantica (Farran, 1908) was identified following the description of Nishida [34].Females were measured under a stereomicroscope (Nikon eclipse 80i, Tokyo, Japan) attached to NIS Elements D softwareVersion 4.13 Nikon, (Tokyo, Japan), using a magnification of 20×/0.4objective lens for O. davisae and 10×/0.25 for O. atlantica.Total length of adult females was measured from the head to the posterior end to the furcal ramus, excluding setae.Confocal images were acquired using an upright Samples of each bongo were preserved differently: one was preserved in ethanol (96%) for molecular analyses, and the twin was preserved in 10% formaldehyde-seawater solution for 24 h and then preserved in ethanol (70%) for morphological studies.
Folsom splitter was used to obtain representative subsamples.Counts in the subsample were transformed into individuals per cubic meter for both total and relative abundance within all individual taxa.At least 500 individuals per sample were identified under a stereomicroscope.Identification was made to the lowest possible taxonomic level.Daily temperature, salinity, and fluorescence data were obtained from a network of oceanographic stations in the same study area and at a depth up to 5 m, provided by the Technological Institute for the Marine Environment Monitoring of Galicia.Upwelling index time series was provided by IEO (www.indicedeafloramiento.ieo.es)accessed on 18 May 2021.River discharge data were available at the River Gauging Network (www.meteogalicia.es)accessed on 19 May 2021.Monthly and daily average and Pearson correlations were calculated for each variable.

Morphological Analyses
For microscopic examination and body measurements, Oithona spp.individuals preserved in formaldehyde were selected and put in a glycerin drop in depression slides.Individuals of O. davisae were identified following the descriptions of Ferrari and Orsi [22] and Temnykh and Nishida [21].O. atlantica (Farran, 1908) was identified following the description of Nishida [34].Females were measured under a stereomicroscope (Nikon eclipse 80i, Tokyo, Japan) attached to NIS Elements D softwareVersion 4.13 Nikon, (Tokyo, Japan), using a magnification of 20×/0.4objective lens for O. davisae and 10×/0.25 for O. atlantica.Total length of adult females was measured from the head to the posterior end to the furcal ramus, excluding setae.Confocal images were acquired using an upright Leica Microsystems (Wetzlar, Germany) SPE laser scanning microscope.Endogenous autofluorescence was detected by exciting the samples with a 488 nm laser, and the emission signal was between 500 and 633 nm.Plan-apochromatic 10× and 20× objectives were used.All images shown correspond to a maximal intensity projection of a confocal z-series.Images were processed with Adobe Photoshop CS (California, EEUU) for Figure 2. Adult females of both species were differentiated using various diagnostic characters: (1) dorsal head profile (the rostrum is dorsally not visible in O. davisae and pointed in O. atlantica) (Figure 2A); (2) lateral rostrum appearance (Figure 2B); (3) number of outer marginal spines on the exopod of swimming legs; and (4) the size of the adult (Figure 2C).

Molecular Analyses
Genomic DNA purification was performed on individual copepods, employing QI-Aamp DNA Micro Kit (QIAGEN, Hilden, Germany), according to the manufacturers' protocol for Isolation of Genomic DNA from Tissues.DNA quality and quantity were checked in a spectrophotometer Nanodrop ® ND-2000 (Thermo Fisher Scientific, Waltham, MA, USA).To identify Oithona species, a fragment of large subunit ribosomal RNA (28S rRNA) gene was amplified using the primers 28SF1: 5 ′ -GCGGAGGAAAAGAAACTAAC-3 ′ and 28SR1: 5 ′ -GCATAGTTTCACCATCTTTCGGG-3 ′ , as described by Ortman [35].PCRs were performed in a total volume of 25 mL containing 1 µL of genomic DNA, PCR buffer at 1x concentration, 0.3 µM primers, 0.2 mM nucleotides, and 0.025 U.µL −1 DreamTaq DNA Polymerase (Thermo Scientific, Waltham, MA, USA).A negative control (no DNA) was also used.The PCR assays were carried out in a Tgradient thermocycler (Biometra, Jena, Germany), under the following reaction conditions [36]: 4 min at 95 • C, 35 cycles of 40 s at 95 • C, 40 s at 50 • C, and 1.5 min at 72 • C, followed by 7 min at 72 • C. PCR products were separated on a 2% agarose gel in Tris acetate EDTA buffer, stained with Red Safe, and scanned in a GelDoc XR documentation system (Bio-Rad Laboratories, Hercules, CA, USA).Subsequently, they were cleaned for sequencing using ExoSAP-IT™ PCR Product Cleanup (Applied Biosystems, Waltham, MA, USA) for 15 min at 37 • C, followed by inactivation for 15 min at 80 • C. Sequencing was performed with a specialized service (StabVida, Caparica, Portugal), and the chromatograms were analyzed using ChromasPro 2.1.8.(Technelysium Pty Ltd., South Brisbane, Australia).All generated sequences were searched for identity using BLAST (Basic Local Alignment Search Tool) through web servers of the National Center for Biotechnology Information (NCBI, Bethesda, MD, USA).All the Oithona 28 rDNA sequences reported in this work have been deposited in the GenBank under the accession numbers: ON114103, ON114104, and ON114106-ON114118.

Molecular Analyses
Genomic DNA purification was performed on individual copepods, employing QIAamp DNA Micro Kit (QIAGEN, Hilden, Germany), according to the manufacturers' protocol for Isolation of Genomic DNA from Tissues.DNA quality and quantity were checked in a spectrophotometer Nanodrop ® ND-2000 (Thermo Fisher Scientific, Massachusetts, USA).To identify Oithona species, a fragment of large subunit ribosomal RNA (28S rRNA) gene was amplified using the primers 28SF1: 5′-GCGGAGGAAAAGAAACTAAC-3' and 28SR1: 5′-GCATAGTTTCACCATCTTTCGGG-3′, as described by Ortman [35].PCRs were performed in a total volume of 25 mL containing 1 µL of genomic DNA, PCR buffer at 1x concentration, 0.3 µM primers, 0.2 mM nucleotides, and 0.025 U.µL −1 DreamTaq DNA Polymerase (Thermo Scientific, Waltham, MA, USA).A negative control (no DNA) was also used.The PCR assays were carried out in a Tgradient thermocycler (Biometra, Jena, Germany), under the following reaction conditions [36]: 4 min at 95 °C, 35 cycles of 40 s at 95 °C, 40 s at 50 °C, and 1.5 min at 72 °C, followed by 7 min at 72 °C.PCR products were separated on a 2% agarose gel in Tris acetate EDTA buffer, stained with Red Safe, and scanned in a GelDoc XR documentation system (Bio-Rad Laboratories, Hercules, CA, USA).Subsequently, they were cleaned for sequencing using ExoSAP-IT™ PCR Product Cleanup (Applied Biosystems, Waltham, MA, USA) for 15 min at 37 °C, followed by in- A sequence of the Cyclopoida species Paracyclopina nana (KR048803) was also included as outgroup.Alignments were performed using Clustal W [37], included in MEGA 7 [38].jModelTest Ver.2.1.10software under the Akaike information criterion (A.I.C.) allowed us to infer the best DNA substitution model for the maximum likelihood (ML) tree [39].The tree was computed with 1000 bootstrap replications.Evolutionary divergence between the 15 28S rDNA sequences obtained from Oithona species from Galicia and the other 96 sequences of Oithona species deposited in GenBank (Table 1) was also calculated by pairwise distance estimation using the maximum composite likelihood method in MEGA 7 software.

Morphological Analyses
Herein we report, for the first time, the coexistence of two species of Copepoda Cyclopoida of the family Oithonidae, O. davisae and O. atlantica, within the study area.The average total length was 0.57 ± 0.03 mm and 1.14 ± 0.09 mm for O. davisae females (N = 30) and O. atlantica females (N = 10), respectively.Male specimens could not be measured due to low abundance; they were almost absent in the study period.
A maximum likelihood (GTR + G model) phylogenetic tree was constructed with 412 sites and showed that the sequences of O. davisae from Ría of Arousa were included in the same clade with O. davisae sequences from the North Sea, with bootstrap values of 100% (Figure 3).However, the species O. atlantica constituted a polyphyletic group, including in the clade sequences of O. plumifera, O. frigida, and O. similis.Our sequences of O. atlantica were placed together, forming a subclade with other O. atlantica sequences from the Atlantic Ocean and Norwegian Sea with low bootstrap values (65%), but phylogenetically distant from two sequences from Azores (KP033174-5) and one sequence from the Atlantic Ocean (JF419543).Sequences of O. atlantica from Azores were grouped in a clade with a 93% bootstrap, being close to the clade of O. similis.The sequence JF419543 of O. atlantica did not cluster in the same clade as other sequences for this species and for the same geographic area, such as JF419541 and JF419542; rather, it was closer phylogenetically to the clade that included the species O. frigida, O similis, and O. atlantica from the Azores, but with a low bootstrap value.Of the eight species of Oithona included in the phylogenetic analysis, seven appeared to be monophyletic supported with high bootstrap values, and only O. atlantica appeared to be polyphyletic.
the Atlantic Ocean (JF419543).Sequences of O. atlantica from Azores were group clade with a 93% bootstrap, being close to the clade of O. similis.The sequence JF4 O. atlantica did not cluster in the same clade as other sequences for this species and same geographic area, such as JF419541 and JF419542; rather, it was closer phy ically to the clade that included the species O. frigida, O similis, and O. atlantica f Azores, but with a low bootstrap value.Of the eight species of Oithona include phylogenetic analysis, seven appeared to be monophyletic supported with hig strap values, and only O. atlantica appeared to be polyphyletic.The pairwise distances calculated between Oithona species from Galicia and other Oithona species deposited in Genbank are shown in Table 2.The analysis involved 424 positions.Intraspecific distances between the five sequences obtained of O. atlantica from Galicia ranged from 0 to 0.012, which is in accordance with those for the O. atlantica sequences deposited (0-0.019).In the study area, the intraspecific distance of O. davisae were lower than those of O. atlantica (0-0.005),being in the same range as those deposited in GenBank (0-0.002).
Thus, the species with the highest intraspecific genetic distances were O. atlantica, O. dissimilis, and O. nana.Interspecific distance confirmed that O. atlantica from Galicia is more similar to O. atlantica, O. plumifera, O. frigida, and also to O. similis, with values lower than 0.043, while the genetic distances with other species were above 0.103.Likewise, O. davisae from Galicia showed interspecific distances as low as 0-0.005 with O. davisae deposited sequences, whereas the values obtained when other Oithona species are compared were higher than 0.117.

Seasonal Abundance Variation
The analysis of seasonal variation was focused on the inner part of the Ría de Arosa (Lombos do Ulla shellfish bed), due to the presence of oceanographic stations in the same area.We examined the relationship between abiotic factors and the abundance of copepods, trying to untangle the limiting factors.
O. atlantica was detected only in October in this area (Figure 4) at very low abundances (10.39 ind.m −3 ).In this area, this species did not present a clear seasonal variation.Abundance of O. davisae was mainly concentrated from May to October 2016, coinciding with higher temperatures in the study area.Its abundance peaked in June (120 ind.m −3 ) and September (108 ind.m −3 ), with an average temperature over 17 • C. On the other hand, density of O. davisae was halved from July, and its individuals disappeared in August.The abundance reached another maximum in September, and then decreased slightly in October.From that month on, the population decreased sharply, reaching minimum abundances, coinciding with low temperatures in the area (below 14 • C).In august, when O. davisae disappeared from the area, the fluorescence was minimal (0.15 mg.m −3 ) (Figure 5b).Salinity fluctuated greatly due to rainfall and river discharge in this inner zone of the estuary (Figure 5c-e).
The greatest presence of O. davisae occurred during the upwelling favorable season in this region (March-April to September-October) (Figure 6a).The upwelling index showed that the periods of maximum abundance of O. davisae were characterized by periods of upwelling relaxation (June, July and October) or downwelling (May and September) in days prior to sampling (Figure 6b).O. davisae did not appear in August, when the upwelling index peaked (5.30 × 10 3 m 3 s −1 km −1 ).Regarding rainfall regime (Figure 6c), when the presence of this species began to be relevant in May, there were several events of intense rains, reaching maximums of 25.2 L.m −2 a few days before sampling.In September, with previous downwelling and high abundance of O. davisae, the highest rainfall reached 27.2 L.m −2 .The greatest presence of O. davisae occurred during the upwelling favor in this region (March-April to September-October) (Figure 6a).The upwe showed that the periods of maximum abundance of O. davisae were character riods of upwelling relaxation (June, July and October) or downwelling (Ma tember) in days prior to sampling (Figure 6b).O. davisae did not appear in Au the upwelling index peaked (5.30 × 10 3 m 3 s −1 km −1 ).Regarding rainfall regime when the presence of this species began to be relevant in May, there were sev of intense rains, reaching maximums of 25.2 L.m −2 a few days before sampli tember, with previous downwelling and high abundance of O. davisae, the h fall reached 27.2 L.m −2 .The Pearson correlation of abundance to each parameter showed that the presence of O. davisae is significantly positively correlated with temperature (r = 0.749).Our data revealed positive correlations with the oceanographic parameters analyzed, but there was an almost null correlation with salinity (r = 0.076).In relation to fluorescence, the correlation was at the limits of significance (r = 0.442).

Expansion of the Two Copepod Species in the Ría de Arosa
Surveys in the Ría de Arosa showed the presence of O. davisae in shellfish beds of Corón (depth 13.8 ± 0.43 m), an area located in the center of the estuary, and in Sarrido (depth 2.97 ± 0.40 m), an outer area of that Ría.Abundance in these shellfish beds was lower than in the inner zone of the estuary (Lombos do Ulla) at depths of 7.15 ± 1.03 m.
Both species (O.davisae and O. atlantica) presented different spatial distributions.While O. davisae was the dominant species of Oithonidae in the innermost zone of the estuary, O. atlantica was more prevalent in the outermost zones.O. atlantica was detected in the inner part of estuary (Lombos do Ulla) only in October (10.39 ind.m −3 ) and in other shellfish beds beginning in May.Its greatest abundance was recorded in Corón in May (41.04 ind.m −3 ) with another peak of abundance in October (29.30ind.m −3 ).In the outermost area (Sarrido), both species presented low abundances (<6 ind.m −3 ).

Discussion
Herein, we confirm the coexistence of invasive copepod Oithona davisae and O. atlantica in the Ría of Arosa (NW Spain), based on morphological and molecular features.Previous studies of the zooplankton community in Galician waters did not record these species [19].The absence of studies of the mesozooplankton community in the Ria of Arosa does not allow us to define when they first appeared and whether O. atlantica is a native species.O. plumifera, morphologically a very similar species to O. atlantica, is included in the Project Checklist of marine species from Galician waters [19].Thus, Oithona plumifera atlantica (Farran, 1913) was even considered a synonymized name [34] of O. atlantica.Additionally, considering that the Galician waters are located within the distribution range of O. atlantica, we could infer a potential confusion between O. plumifera and O. atlantica, since most of the studies in these waters were based on morphological characters, with an absence of molecular studies.
The species O. nana, O. similis, and O. plumifera, described so far in the Ría of Vigo and the waters of A Coruña (by the monitoring program RADIALES), were not detected in our area and time of study.Due to the sampling design in shallow shellfish beds, the mesh size used heavily underestimates small Oithona, and the short study period prevents us from concluding that they have been displaced by the presence of the invasive species O. davisae.
DNA-based analyses have become consolidated as part of biodiversity studies and have proven to be powerful tools for assessing the diversity of marine zooplankton species.For molecular analysis, we decide to use the large-subunit (28S) rRNA gene, because it has proven to be a useful marker to identify and discriminate among Oithona species [36,40,41] and to identify invasive species [24].DNA extraction from Oithona species was difficult due to the small size of the individuals.Maximum likelihood analysis of the 28S rRNA gene points out that O. davisae from Galicia formed a well-supported monophyletic clade together with sequences of O. davisae from the North Sea, as described in previous studies [24,40].Likewise, the phylogenetic tree showed that the sequences of O. atlantica obtained in this study were grouped with those of O. atlantica from the Norwegian Sea and Atlantic Ocean in a polyphyletic group, which also included O. plumifera, O. similis, and O. frigida species.Our tree revealed several well-supported clades for Oithona species but clearly did not solve the phylogeny of O. atlantica.Cornils et al. [40] described a similar phylogenetic topology, and it was supported by the interspecific genetic distances, which were significantly lower between the species of the O. atlantica group than between O. atlantica and the other Oithona species.Intraspecific genetic distances obtained for both O. davisae and O. atlantica sequences also agree with those calculated in other studies [36,40].
The existence of cosmopolitan species, living In at least two oceans, is widely accepted in marine plankton ecosystems due to the lack of apparent geographic barriers and potentially broad oceanic dispersal.In recent years, it has been shown that some of the Oithona species are circumglobally distributed and constitute a complex of cryptic species [40], demonstrating that O. similis is a complex of at least seven distinct mitochondrial lineages.Ueda et al. [41] also suggested that O. dissimilis is a species complex containing more than two cryptic species.For O. atlantica, the absence of a monophyletic clade and its mixture with other species seem to suggest that a species complex could be housed under its denomination, so it would be interesting to carry out population genetic studies for this species worldwide.
The main objective of the project that involved this zooplankton sampling was to detect which zooplankton organisms acted as intermediate hosts for the Bivalvia parasite.For this reason, the sampling design involved the use of 200 µm mesh size, necessary for the collection of the mesozooplankton targeted (0.2 to 20 mm in size).The optimal mesh size for the estimation of O. davisae would be 100 µm [42], and thus the mesh size used in this work was not completely appropriate for a correct and complete study of the life cycle of O. davisae.Nauplii, juveniles, and adult males are not retained in a 200 µm mesh size, which would explain why these smaller individuals were absent in the samples.However, we consider that its presence throughout the study period and its temporal variability are of interest since it is an invasive species that has not been described in this study area.
Regarding the temporal distribution, O. davisae was present throughout the study period (November 2015-March 2017).The highest abundances were observed in the warm season (May-October).The maximum abundance in June coincides with previous studies by O. davisae in other geographic areas such Venice Lagoon [28].Other studies, in Sevastopol Bay [29] and Estuary of Bilbao [32], showed maximum values in autumn, and in the Aegean Sea [43], higher abundances were observed during summer months (July, September) and in winter (February).A second peak of abundance was recorded in late summer (September, 108.9 ind.m −3 ).In our study, the coldest season showed a quasi-absence of this species.We have to consider that we cannot compare the abundances obtained in our study with other works ( [44] in southern Black Sea with maximum of 49,761 ind.m −3 at the coastal station; Pansera et al. [28] in Venice Lagoon with records of 6490 ind.m −3 at the inner station), because the mesh size used is different and our abundances are underestimated.Seawater temperature seems to be the main oceanographic parameter that influences the abundance of O. davisae (r = 0.749), with peaks of abundance in periods of higher temperatures and minimums when the temperature was less than 14 degrees.This is consistent with the thermophilic character of O. davisae.An increase in seawater temperature favors the appearance and settlement of thermophilic invasive species.
Fluorescence also appears to be related to increasing population density.During the present study, fluorescence reached two maximums, in June (9.47mg.m −3 ) and in October (6.36 mg.m −3 ), coinciding with its peaks of abundance.These results suggest that there is a correlation between fluorescence and maximum abundance of O. davisae (r = 0.442).This coincidence seems to be related to the omnivorous feeding of this species.Previous studies [45][46][47][48] showed that Oithona species feed mainly on ciliates and dinoflagellates and that O. davisae avoids diatoms [49].
Monitored oceanographic conditions allow us to shed light on the seasonal dynamics of O. davisae.Although the greatest presence of this invasive species is concentrated in the warm season, it is necessary to reduce the time scale (monthly to daily) to untangle its dynamics in the study area, especially in species with such short life cycles and rapid growth.It was observed that higher densities during the warm period correspond with upwelling relaxation-downwelling in the Ría of Arosa coupling temperature and upwelling index data in the days prior to sampling.A large number of studies have analyzed the hydrodynamic characterization of the Galician Rías and this particular coastal upwelling system [50,51].Additionally, previous studies have investigated the influence of upwelling-downwelling events in seasonal plankton dynamics [12,[52][53][54][55].The presence of O. davisae presented herein seems to be directly related to upwelling relaxation or downwelling events, and consequently with seawater temperature.Considering various factors (upwelling index, temperature, river discharge, and rainfall), the data suggest that retention and accumulation of seawater in the innermost area is favored during rainfall events (May, June, and September; Figure 6c).This event slows down the seawater renewal, and therefore favors warming up the surface water.These conditions of high seawater temperatures are optimal for the growth and reproduction of O. davisae [56].Furthermore, upwelling relaxation and downwelling events favor the dominance of dinoflagellates in phytoplankton communities in Rías of Galicia [57,58].Conversely, upwelling events favor the prevalence of diatoms.As we have already mentioned, dinoflagellates constitute the main food of O. davisae.Other species of the same genus, such as O. plumifera, are considered an upwelling-indicative species [52].Our results suggest that O. davisae, despite being a species adapted to a wide range of temperatures [29,59,60], prefers upwelling relaxation and downwelling events.Furthermore, in recent years the average intensity of upwelling and the length of the upwelling season have been reduced [61], which may favor the settlement of these invasive species from warmer areas.
Our results showed that the presence of O. davisae in terms of abundance is limited by physical (upwelling-downwelling, seawater temperature) and biological variables such as food availability (dinoflagellates) and predation by gelatinous plankton [56] and by fish larvae [62].According to the available literature, O. davisae is more abundant in inner areas [28,42].In our series, we observed the expansion of this invasive species from the inner zone of the estuary (Lombos do Ulla) to other more external areas of the Ría of Arosa (Corón and Sarrido), but to a lesser extent.Consequently, the detection of O. davisae in these outermost areas seems to indicate that this species is adapting to different hydrographic conditions, according to the high tolerance of this species to variability of environmental conditions [29].This high tolerance can become a danger to the ecosystem because this species seems to be able to establish itself successfully in a new environment [24].Besides this, the short sampling period and the mesh size used do not make it possible to determine whether this invasive species has displaced native species of the same genus Oithona.
In this study, O. davisae and O. atlantica presented an opposite pattern.O. atlantica was abundant in more external areas with higher densities in the middle area (41.04 ind.m −3 in May and 29.30 ind.m −3 in October).The reason why the abundances are higher in the middle zone (Corón) than in the outermost zone (Sarrido) of Ría de Arosa, may be that the outermost zone where the samplings were carried out was shallow (2.97 ± 0.40 m).If we observe its geographical distribution, O. atlantica could be considered a typical species of these waters.Furthermore, other studies [63] have shown that the copepod size composition can decrease with eutrophication and in confined coastal waters [64].This fact could also explain the greater abundance of O. davisae in the inner part of the Ría de Arosa.
In the sampling area, ballast water and bivalve transplants are probably the main routes of entry of invasive species.Organisms included in this ballast water and originating from a distant port of origin are represented by bacteria, protozoa, and invertebrates, mainly as resting stages able to survive travel conditions [65] and establish a new reproductive population in the new ecosystem.Additionally, bivalve transplants activities are also frequent, since the three sampling areas are important cockle and clam shellfish beds in Galicia.In fact, bivalve seed is imported from countries such as France or Italy, Mediterranean countries where O. davisae is settled.Moreover, the introduction of a new species entails additional implications if we assimilate the holobiont concept, which includes the host and all its symbionts.It is well known that copepods act as intermediate hosts in the biological cycle of numerous parasitic species [66][67][68][69][70][71][72][73].The co-introduction of parasitic pathogens in wild shellfish production systems can cause unexplained mortality phenomena in populations of bivalve mollusks, with its consequent impact on native species, ecosystems, and ecosystem services [74][75][76].Undertaking interdisciplinary studies that include the assessment of the impacts of the invasion by any particular given holobiont is a very challenging topic in socio-ecological system management.

Figure 3 .
Figure 3. Maximum likelihood (GTR + G model) analysis of 28S rRNA gene showing the ta position of Oithona species from Ría de Arousa.Numbers at branch nodes indicate boots

Figure 3 .
Figure 3. Maximum likelihood (GTR + G model) analysis of 28S rRNA gene showing the taxonomic position of Oithona species from Ría de Arousa.Numbers at branch nodes indicate bootstrap confidence values in percent.Sequences obtained in this study: ON114103-ON114104, ON114106-ON114118.

Figure 4 .
Figure 4. Seasonal variations of O. davisae abundance and water temperature.Figure 4. Seasonal variations of O. davisae abundance and water temperature.

Figure 4 .
Figure 4. Seasonal variations of O. davisae abundance and water temperature.Figure 4. Seasonal variations of O. davisae abundance and water temperature.

Figure 4 .
Figure 4. Seasonal variations of O. davisae abundance and water temperature.

Figure 6 .
Figure 6.Upwelling index: (a) interannual variability, (b) UI monthly variation (upwelling index) during the period of maximum abundance and average daily temperature (temperature records), (c) rainfall.

Figure 6 .
Figure 6.Upwelling index: (a) interannual variability, (b) UI monthly variation (upwelling index) during the period of maximum abundance and average daily temperature (temperature records), (c) rainfall.

Table 1 .
28S rDNA sequences of Oithona species from Galicia and others deposited in GenBank used to estimate the genetic divergence.N = number of sequences.

Table 2 .
Genetic divergence between 28S rDNA sequences of Oithona atlantica/O.davisae from Galicia (Gal) and sequences of other Oithona species deposited in GenBank.Intraspecific distances are also shown.