Assessing Environmental Control on Temporal and Spatial Patterns of Larval Fish Assemblages in a Marine Protected Area

: The northern coast of the Iberian Peninsula is an important spawning and nursery area for several marine ﬁsh species, some of which are economically exploited by ﬁsheries and under management plans. Larval stages of ﬁsh are highly sensitive to environmental change and anthropogenic pressures, and Marine Protected Areas (MPA) can help mitigate the impacts on ﬁsh populations. This study investigated the environmental drivers of the temporal and spatial patterns of the larval ﬁsh assemblages inhabiting a small coastal MPA along the NW Portuguese Iberian Coast. Seasonal surveys were conducted over two years at nine sampling stations distributed throughout the MPA to collect larval ﬁsh samples and water parameters. Results showed that a total of 39 different ﬁsh taxa were identiﬁed. In terms of abundance, reef-associated species, such as Parablennius gattorugine (54.6%), and marine species that use estuaries as nursery areas, such as Ammodytes tobianus (15.7%) and Clupeidae n.i. (8.8%) dominated the larval ﬁsh assemblages. The larval ﬁsh assemblages were characterized by a strong temporal pattern that, according to CCA analyses, was related to the temporal variability of water temperature, pH, chlorophyll α , TPM, and also the river ﬂow of an adjacent river. This study showed that 47% of the ﬁsh larvae belonged to commercially exploited species, highlighting the importance of this MPA as a nursery area for the early life stages of the ﬁsh population. Overall, these new ﬁndings emphasize the role of MPAs in ensuring the connectivity of ﬁsh species between marine-estuarine habitats and enhancing the conservation of ﬁsh populations.


Introduction
Marine protected areas (MPAs) are created with the purpose of protecting natural habitats and conserving biodiversity and abundance levels within the protected areas [1]. They are considered important management tools to help conserve and protect marine ecosystems from negative impacts of human activities, such as over-exploitation of fish stocks, underwater mining, and tourism-related pressures, or as a tool for species conservation [1][2][3][4]. MPAs may have different levels of protection, with different limitations of human activities. For example, in some MPAs, all fishing activities are prohibited in the so-called "no-take zones", while in other MPAs, only some restrictions on fishing activities may be applied [5]. MPAs are used as an important tool for protecting commercially important fish species by improving abundance and mean fish size inside the MPA [3,6,7]. In addition, MPAs positively contribute to fish larvae abundance inside protected areas, which will also improve abundances in adjacent unprotected areas due to fish larvae and egg dispersal [1,[8][9][10]. However, these effects are dependent on the MPA size and may take some years to be noticeable [2,11].
Globally, MPAs represent 6.35% of the ocean's total area; however, only 1.89% are defined as no-take MPAs, where no fishing, mining, or other extractive activities are allowed.

Study Area
The present study was undertaken in the Parque Natural Litoral Norte (PNLN), an MPA located on the NW Portuguese coast. Created in 2005 to protect the coastal region against anthropogenic impacts such as illegal construction on the dunes, the park extends 16 km from the Neiva River to the south of Apúlia. The PNLN extends 2.5 nautical miles from the intertidal areas until 40 m depth, at the most, offshore limits. Along the marine domain of the PNLN, the seabed is mostly composed of rocky reefs, with small soft-substrate areas between the reefs. It is located at the northernmost limit of the Canary Upwelling System, where seasonal upwelling, important to the pelagic ecosystems and local productivity, can be observed during spring and summer [39,40]. Another relevant feature is the presence of the Western Iberia Buoyant Plume, which originates from the discharge of freshwater by rivers along the NW Iberian coast [41,42], which also contributes to coastal productivity. At the center of the PNLN is the Cávado River estuary, and further up north is the Lima estuary ( Figure 1). Both estuaries play an important role in this coastal region, providing productivity to this region and functioning as important nursery areas for many fish species, some with high economic value [43,44]. Apart from its conservation status, there is still little information about the marine area of the PNLN, namely its biodiversity and major ecological functions, making this study the first approach to understanding the dynamics of PNLN planktonic assemblages.
However, its role for fish species is still not well understood, and the benefits of this MPA for enhancing fish stocks are still not fully known [38]. This study aims to fill this knowledge gap by investigating, for the first time, the larval fish assemblages of the PNLN, namely: (i) describe seasonal and spatial patterns of several abiotic parameters; (ii) characterize the PNLN larval fish assemblages; and (iii) investigate the influence of environmental parameters on the PNLN larval fish assemblages.

Study Area
The present study was undertaken in the Parque Natural Litoral Norte (PNLN), an MPA located on the NW Portuguese coast. Created in 2005 to protect the coastal region against anthropogenic impacts such as illegal construction on the dunes, the park extends 16 km from the Neiva River to the south of Apúlia. The PNLN extends 2.5 nautical miles from the intertidal areas until 40 m depth, at the most, offshore limits. Along the marine domain of the PNLN, the seabed is mostly composed of rocky reefs, with small softsubstrate areas between the reefs. It is located at the northernmost limit of the Canary Upwelling System, where seasonal upwelling, important to the pelagic ecosystems and local productivity, can be observed during spring and summer [39,40]. Another relevant feature is the presence of the Western Iberia Buoyant Plume, which originates from the discharge of freshwater by rivers along the NW Iberian coast [41,42], which also contributes to coastal productivity. At the center of the PNLN is the Cávado River estuary, and further up north is the Lima estuary ( Figure 1). Both estuaries play an important role in this coastal region, providing productivity to this region and functioning as important nursery areas for many fish species, some with high economic value [43,44]. Apart from its conservation status, there is still little information about the marine area of the PNLN, namely its biodiversity and major ecological functions, making this study the first approach to understanding the dynamics of PNLN planktonic assemblages.

Sampling and Data Collection
A dedicated sampling strategy was designed to characterize for the first time the larval fish assemblages of PNLN and investigate the role of environmental parameters on the temporal and spatial dynamics of fish larvae in both commercially and non-commercially exploited fish species. Environmental parameters studied included abiotic parameters, such as temperature, salinity, pH, dissolved oxygen concentration, oxygen saturation, turbidity, nitrate (NO 3 ), nitrite (NO 2 ), phosphate (PO 4 ), ammonium (Nh 4 ), silicate (Si), total particulate matter (TPM), particulate organic matter (POM), river flow of the nearby Cávado estuary, upwelling index, and biotic parameters such as chlorophyll a and zooplankton abundance, both used as proxies of larval fish prey.
Seasonal sampling was performed (winter, spring, summer, autumn) from Autumn 2017 to Summer 2019 at nine sampling stations homogenously distributed along three transects perpendicular to the coast (North, Center, and South) of the PNLN (Figure 1). Sampling campaigns were performed once for each season, as follows: Autumn 2017 (October), Winter 2018 (February), Spring 2018 (May), Summer 2018 (July), Winter 2019 (January), Spring 2019 (March), and Summer 2019 (June). All nine sampling stations were sampled during daylight on the same day. Ichthyoplankton samples were collected at each sampling station, and zooplankton samples were taken at the station located in the middle of each transect. For ichthyoplankton sampling, a 1 m diameter, 4 m long, and 500 µm mesh size net was used. Subsurface (1-2 m depth) tows were performed at a constant velocity of ca. 1 ms −1 for 10 min; this sampling methodology was chosen so it could be compared with other studies conducted in nearby areas (e.g., [45,46]). Zooplankton samples were collected using a 150 µm mesh size net for 1 min, also at a depth of 1 to 2 m. After collection, all samples were immediately fixed in 4% buffered formalin and stored. The volume of filtered water was measured by a flowmeter (Hydro-Bios) attached to each plankton net. Due to climactic and navigability constraints, it was not possible to sample in Autumn 2018.
To get a better understanding of the environmental conditions along the water column at each sampling station, vertical profiles of several abiotic parameters were also measured, along with the collection of surface and bottom water samples, with a Van Dorn bottle. Water samples were stored in coolers, transported to the laboratory, and processed for further analytical quantification of nutrients (NO 3 , NO 2 , PO 4 , Nh 4 , Si), chlorophyll a, TPM, and POM. In situ measurements of physical-chemical parameters of the water column (temperature, salinity, pH, dissolved oxygen concentration, oxygen saturation, and turbidity) were performed at each sampling station using a multiparameter probe (YSI EXO1 Sonde) along the water column producing a vertical profile for each measured parameter. Due to a multiparameter probe misfunction, no environmental data were available for the summer of 2019.
Upwelling index values were provided by the Spanish Institute of Oceanography (http://www.indicedeafloramiento.ieo.es accessed on 17 February 2023) and were generated using the sea level pressure of the FNMOC (https://www.metoc.navy.mil/fnmoc/ fnmoc.html accessed on 17 February 2023) to determine the geostrophic wind, at a position of 41 • N latitude by 10 • W longitude. Positive values correspond to upwelling events, while negative ones indicate downwelling. A variable designated as "Upwelling2m" was created, corresponding to the mean value of the upwelling value during the sampling time and the month previous to the sampling time to better explain the upwelling conditions affecting fish larval assemblages. River flow data for the Cávado River was obtained from Caniçada dam, using the Portuguese Sistema Nacional de Informação de Recursos Hídricos (https://snirh.apambiente.pt/ accessed on 17 February 2023).

Laboratory Processing
Ichthyoplankton samples were sorted and identified to the highest taxonomic classification possible using ichthyoplankton identification keys [47][48][49][50] under a stereo zoom microscope (Nikon SMZ 800). The number of individuals was standardized to the number of fish larvae per 100 m 3 of filtered water. Zooplankton samples were counted and identified to the main taxonomic groups using a Bogorov counting chamber and sub-sampled in 2 mL samples, and then standardized to the number of individuals per m 3 of filtered water. Chlorophyll a concentration was determined spectrophotometrically by filtering water samples using 45 µm cellulose acetate membrane filters and 90% acetone for posterior extraction [51] with cell homogenization using the SCOR-UNESCO [52] trichromatic equation. Quantification of silicate, nitrite, phosphate, and ammonium concentrations were determined by spectrophotometry using the methodology described by Grasshoff et al. [53]. Nitrate values were measured using an adaptation of the spongy cadmium technique by Jones [54]. For TPM and POM analysis, water samples were filtered through pre-combusted GF/F glass-fiber filters, dried overnight at 100 • C for TPM quantification, and later combusted at 500 • C for POM quantification [55].

Data Analyses
Larval fish diversity was measured by the Shannon-Wiener index (H') [56] and was calculated using the diversityfunction in R from the vegan package, and the equitability was measured by Pielou's evenness index (J ) [57] and was calculated using the following equation where s represents the total number of species: These indices were calculated using all identified taxa, namely species, genus, and family (e.g., Labridae n.i).
To better understand the connectivity between PLNL and nearby estuaries, species were assigned into different ecological guilds based on their estuarine use [58], namely: Estuarine Species (ES), Marine Migrants (MM), Marine Stragglers (MS), and OTHER when no information was found. Fish species were also classified as exploited (hereafter designated as commercial species) or not exploited (hereafter designated as non-commercial species) by commercial fisheries, according to official Portuguese fisheries data [59].
Larval fish assemblages descriptors (abundance, diversity, and ecological guilds) and environmental variables along the PNLN during the different sampling times were mapped using QGIS 3.10.8. For the temporal and spatial characterization of the different environmental variables, continuous layer maps were created using a deterministic method, the inverse distance weighting interpolation.
Taking into consideration that at temperature latitudes, larval fish assemblages are characterized by strong seasonal patterns (e.g., [24]), temporal and spatial patterns of the larval fish assemblages, and water characteristics were assessed by two-way ANOVAs. To assess differences in the abundance and diversity of larval fish assemblages and also on the environmental variables, two-way ANOVAs were performed, studied, and investigated: one with sampling time and transect as fixed factors; and another with sampling time and distance from the coast as fixed factors. Although there were no replicates of planktonic tows, in order to have at least three replicates for the statistical analysis ANOVA, the three tows of each horizontal transect (north, center, south) were considered as replicates for factor Distance and tows from each vertical transect (d1, d2, d3) were considered as replicates for factor Transect. Sampling time had seven groups (Autumn 2017, Winter 2018, Spring 2018, Summer 2018, Winter 2019, Spring 2019, and Summer 2019), transect had three groups (North, Center, and South), and distance had three groups (1, 2, 3) from closest to furthest from the coast. To study the temporal and spatial effects of the environmental variables, two two-way ANOVAs were performed using the same design described previously, with the exception of the environmental parameters where no data could be obtained for Summer 2019 (Oxygen concentration, Oxygen saturation, and Turbidity), in this case, the fixed factor sampling time was composed by only 6 groups. Two-way ANOVAs were also used to investigate temporal and spatial patterns between commercial and non-commercial species, considering sampling time, transect, and distance as fixed factors. Since none of the ANOVA assumptions were met, the temporal differences of the Cávado River flow were investigated using the non-parametric Kruskal-Wallis test. Upwelling2m differences were tested using a one-way ANOVA with sampling time as a fixed factor. Prior to statistical analysis, variables were tested for homogeneity (Levene's test) and checked for normality. Chlorophyll a and NH 4 data were transformed [Ln(x)] in order to stabilize the variance and to fit data to a normal distribution, fulfilling at least one of the ANOVA assumptions. For those cases, ANOVA results were only accepted where significance levels were <0.01; otherwise, ANOVA results were accepted at significant levels of <0.05. Furthermore, in the event of significance, a posthoc Tukey HSD for unequal sample sizes was used to determine which means were significantly different (p < 0.05) [60]. Kruskal-Wallis test, ANOVAs, and posthoc analysis were performed using the kruskal.test, aov and TukeyHSD functions in R from the native stats library.
To investigate spatial and temporal variations in the structure of the larval fish assemblages, a 2-way analysis of similarity (ANOSIM) [61], based on a Bray-Curtis similarity matrix, was used. In one ANOSIM analysis, sampling time and transect were considered fixed factors, and in another ANOSIM analysis, sampling time and distance from the coast were considered fixed factors. To assess the contribution of each species to the differences in fish larvae assemblages along the different sampling times and transects studied, similarity percentages (SIMPER) [62] were calculated. Fish larvae data were square-root transformed, and only species with more than 0.1% of abundance were considered to avoid rare species interference. Non-metric multidimensional scaling (MDS) based on the Bray-Curtis similarity matrix [63] was carried out using non-transformed data. ANOSIM, SIMPER, and MDS were performed using PRIMER v.6 (Plymouth Routines Multivariate Ecological Research).
The influence of environmental variables on the fish larvae assemblages between the different sampling times was investigated using canonical correspondence analysis (CCA) [64] with the software CANOCO v4.5 (Microcomputer Power, Ithaca, NY, USA). Larval abundances were transformed (square-root), and only species with an abundance higher than 0.05% and a frequency of occurrence higher than 2% were included in the analysis to avoid the undue effect of very rare species. Biplot scaling with a focus on the interspecies distances option was used, and the significance of the canonical model was given by a Monte Carlo test [65]. Inter-set correlation coefficients were used to assess the importance of the environmental variables; when inter-set ≥ 0.4, variables were considered biologically important [23]. The environmental variables used in CCA were: temperature, salinity, oxygen saturation, dissolved oxygen, pH and turbidity, chlorophyll a, nitrate, nitrite, ammonium, phosphate, silicate, TPM and POM, upwelling, and river flow. For temperature, salinity, oxygen saturation, dissolved oxygen, pH, and turbidity, the mean of the vertical profile of the water column was considered, while for chlorophyll a, nitrate, nitrite, ammonium, phosphate, silicate, TPM, and POM, it was considered the mean of surface and bottom values.

Environmental Variables
Overall, water column parameters were relatively similar throughout the study area and had temporal variations ( Figure 2, Table 1, Supplementary Material Figures S1 and S2, Table S1). The water temperature had a significant temporal variation (ANOVA F = 128. 93 Table 1). Temperature also varied significantly between transects (ANOVA F = 3.39; p < 0.05), with the Center transect exhibiting significantly higher temperatures than the North (Tukey HSD p < 0.05). Regarding distance to the coast, no significant differences were observed (ANOVA F = 1.11; p = 0.34). Along the water column, Temperature was stratified, with a stronger temperature gradient in summer and the inverse in winter (Supplementary Material Figure S1). Salinity varied significantly between sampling times (ANOVA F = 20.64; p<0.001) and transects (ANOVA F = 4.00; p < 0.05), with higher levels in the South transect and during Autumn 2017 and Spring 2019 while, lower levels were registered during Spring 2018 and Winter 2019 (Tukey HSD p < 0.01) ( Figure 2b, Table 1). Salinity had almost no vertical stratification in summer; inversely, in the winter, there was a brackish water layer at the surface, and salinity increased with depth (Supplementary Material Figure S1). Overall, water oxygenation was good (Table 1). Oxygen saturation and dissolved oxygen concentration showed temporal variations (ANOVA F = 24.62; F = 21.04; p < 0.001), with higher values in summer 2018 (Tukey HSD p < 0.05) (Figure 2c,d, Table 1). No significant spatial differences were observed. Along the water column, both oxygen saturation and dissolved oxygen showed higher values near the surface, with values decreasing with depth (Supplementary Material Figure S1). pH varied significantly between sampling times (ANOVA F = 65.42; p < 0.001), with higher values during Spring and Summer (Tukey HSD p < 0.01). pH also showed significant differences between the South and Center transects (ANOVA F = 3.89; p < 0.05), while no differences were observed with distance to the coast (ANOVA F = 0.27; p = 0.76) (Figure 2e, Table 1). Overall, turbidity of the water column was low (Table 1, Figure 2f); however, it varied significantly between sampling times (ANOVA F = 11.03; p < 0.001) (Figure 2f) since significantly higher turbidity values were registered in Autumn 2017, Winter, and Spring 2018 (Tukey HSD p < 0.05). Turbidity also showed significant variation related to distance to the coast (ANOVA F = 12.37; p < 0.001), with the near-shore sampling stations showing significantly higher turbidity (Tukey HSD p < 0.01). No differences were observed related to the different studied transects (ANOVA F = 011; p = 0.87). Turbidity values tended to be higher at the surface and lower in the middle/bottom of the water column (Supplementary Material Figure S1). Nutrient concentrations were generally low (Table 1). Nitrates (NO3) had a significant temporal variation, with higher values during Autumn 2017, Spring, and Winter 2019 (Tukey HSD p < 0.05), with no significant differences in nitrate values between transects or distance from the coast. Ammonium (Nh4) values did not vary temporally or spatially; however, a peak was observed near shore during the Summer of 2018, reaching 13.20 ± 11.38 (Supplementary Material Figure S2). Silicate (Si)  Nutrient concentrations were generally low (Table 1). Nitrates (NO 3 ) had a significant temporal variation, with higher values during Autumn 2017, Spring, and Winter 2019 (Tukey HSD p < 0.05), with no significant differences in nitrate values between transects or distance from the coast. Ammonium (Nh 4 ) values did not vary temporally or spatially; however, a peak was observed near shore during the Summer of 2018, reaching 13.20 ± 11.38 (Supplementary Material Figure S2). Silicate (Si) concentrations were significantly higher during Winter 2019, and the Center transect had significantly higher silicate values throughout the sampling times. There were no significant differences in silicate values along the distance from the coast (Tukey HSD p < 0.05) (Supplementary Material Figure S2). Nitrite (NO 2 ) levels showed significant temporal differences with lower values during the Summer of 2019 (Tukey HSD p < 0.05) (Supplementary Material Figure S2). Between transects, values of NO 2 were significantly higher at the south transect compared to the north transect (Tukey HSD p < 0.05). Phosphate (PO 4 ) had significant differences between sampling times; levels were non-significant (Tukey HSD p > 0.05) when comparing the same seasons from different years with the exception of winter samples, where Winter 2018 had significantly higher values than Winter 2019 (Tukey HSD p < 0.001) (Supplementary Material Figure S2); no spatial variation was observed. Both TPM and POM had significant temporal variation, with higher concentrations in summer and spring (Tukey HSD p < 0.05); no significant spatial variation was observed.
Cávado River flow had a significant temporal variation (p < 0.001), with higher values of river flow occurring during Winter samplings. Upwelling2m had no significant variation between sampling times (F = 0.41; p = 0.85); also, all sampling times occurred during upwelling events with the exception of Summer 2019 (−31.11 m 3 s −1 km −1 ) (Figure 3b).

Larval Fish Assemblages
During the study, a total of 4170 fish larvae were collected, 39 different taxa from 16 families were identified, and 47% of the identified taxa were commercial species ( Table 2). The most abundant families were Blenniidae (59.7%), Ammodytidae (15.7%), and Clupeidae (8.8%), totaling 84.2% of the fish larvae collected. The three most abundant taxa were Parablennius gattorugine (54.6%), followed by Ammodytes tobianus (15.7%), and Clupeidae n.i. (8.8%), responsible for 79.1% of all fish larvae collected. Regarding the frequency of occurrence (FO), P. gattorugine was the most frequently observed taxa, occurring in 66.7% of the total 63 samples collected. Clupeidae n.i. was the second most frequent taxa with a FO of 46.0%, followed by A. tobianus and Parablennius ruber, which appeared in 30.2% of the total samples collected. Of the 39 identified taxa, 12 taxa were rare (i.e., only observed one time), corresponding to a FO of 1.6% (Table 2). Unidentified fish larvae (n.i.) accounted for 7.9% of the total larvae collected, mostly due to a bad conservation state and/or yolk-sac larvae phase. variation was observed. Both TPM and POM had significant temporal variation, with higher concentrations in summer and spring (Tukey HSD p < 0.05); no significant spatial variation was observed.
Cávado River flow had a significant temporal variation (p < 0.001), with higher values of river flow occurring during Winter samplings. Upwelling2m had no significant variation between sampling times (F = 0.41; p = 0.85); also, all sampling times occurred during upwelling events with the exception of Summer 2019 (−31.11 m 3 s −1 km −1 ) ( Figure  3b).

Larval Fish Assemblages
During the study, a total of 4170 fish larvae were collected, 39 different taxa from 16 families were identified, and 47% of the identified taxa were commercial species ( Table 2). The most abundant families were Blenniidae (59.7%), Ammodytidae (15.7%), and Clupeidae (8.8%), totaling 84.2% of the fish larvae collected. The three most abundant taxa were Parablennius gattorugine (54.6%), followed by Ammodytes tobianus (15.7%), and Clupeidae n.i. (8.8%), responsible for 79.1% of all fish larvae collected. Regarding the frequency of occurrence (FO), P. gattorugine was the most frequently observed taxa, occurring in 66.7% of the total 63 samples collected. Clupeidae n.i. was the second most frequent taxa with a FO of 46.0%, followed by A. tobianus and Parablennius ruber, which appeared in 30.2% of the total samples collected. Of the 39 identified taxa, 12 taxa were rare (i.e., only observed one time), corresponding to a FO of 1.6% (Table 2). Unidentified fish larvae (n.i.) accounted for 7.9% of the total larvae collected, mostly due to a bad conservation state and/or yolk-sac larvae phase.    The larval fish assemblages of the PNLN included the three different ecological functionals with twenty-three taxon identified as Marine Stragglers (MS), five as Marine Migrants (MM), six as Estuarine Species (ES), and four had no group assigned (Table 2). MS was mainly composed of representatives of the Blennidae family and A. tobianus, accounting for 13 species classified as non-commercial and 10 as commercial species. MM included only commercial taxa such as Clupeidae n.i., Diplodus sargus, and Atherina presbyter ( Table 2).
Pomatoschistus microps and P. minutus were the main representatives of ES, and none of these species were considered commercial (Table 2)   ANOSIM results revealed significant differences in the structure of the PNLN larval fish assemblages between sampling times (Table 3), with groups clearly separated (R = 0.836 p < 0.001). Pairwise results showed that all sampling times had a significantly different assemblage composition (p < 0.05), with groups clearly separated (R > 0.5), with the exception of Autumn 2017 and Summer 2018, when the structure of the larval fish assemblages was not considered significantly different (R = 0.370). In contrast, the structure of the larval fish assemblages did not vary spatially (R < 0.05) ( Table 3). SIMPER results showed that the three most abundant species were the main species responsible for the temporal dissimilarities. Higher abundances of P. gattorugine occurred during Spring in 2018 and 2019 and were the main species responsible for dissimilarities related to these sampling times, while A. tobianus and Clupeidae n.i. were responsible for dissimilarities related to Winter 2018 and Summer 2019, respectively.    MDS plots revealed clear temporal variation with samples clustering according to their sampling time (Figure 6a). Spring samples, characterized by high abundances of P. gattorugine, were clearly separated from other sampling times (Figure 6b), while Winter 2018 where A. tobianus abundance was at its highest abundance was isolated from all the other sampling times (Figure 6c

Environmental Control of Fish Larval Assemblages
The influence of environmental variables on the larval fish assemblages was investigated using canonical correspondence analysis (CCA), which showed that species were distributed along the first two axes of the CCA. The first axis had an eigenvalue of 0.6, and the second axis also had an eigenvalue of 0.6, with both axes exhibiting high species-environment correlations of 0.9. A Monte-Carlos permutation test showed that the effect of the environmental variables on the explained distribution of the CCA axes was significant (F = 2.73; p < 0.01). According to the inter-set correlation coefficients, temperature, chlorophyll a, and TPM were positively related to the first CCA axis, while river flow was negatively correlated. PO4 and river flow were positively correlated with the second CCA axis, while pH and TPM were negatively correlated (Table 4).

Environmental Control of Fish Larval Assemblages
The influence of environmental variables on the larval fish assemblages was investigated using canonical correspondence analysis (CCA), which showed that species were distributed along the first two axes of the CCA. The first axis had an eigenvalue of 0.6, and the second axis also had an eigenvalue of 0.6, with both axes exhibiting high speciesenvironment correlations of 0.9. A Monte-Carlos permutation test showed that the effect of the environmental variables on the explained distribution of the CCA axes was significant (F = 2.73; p < 0.01). According to the inter-set correlation coefficients, temperature, chlorophyll a, and TPM were positively related to the first CCA axis, while river flow was negatively correlated. PO 4 and river flow were positively correlated with the second CCA axis, while pH and TPM were negatively correlated ( Table 4).
The ordination of the samples showed that they tended to cluster in accordance with sampling times, revealing a strong temporal gradient defined by the CCA axes (Figure 7a). Winter 2018 and 2019 samples with higher levels of river flow and PO 4 clustered on the negative side of the first CCA axis and on the positive side of the second CCA axis. Spring 2018, Spring 2019, and Summer 2019 were characterized by higher temperature, chlorophyll a, TPM, and pH levels clustered on the positive side of the first CCA axis and the negative side of the second CCA axis (Figure 7a). Species ordination showed that ES clustered in the positive part of the first CCA axis since they were positively correlated with temperature and TPM. However, they were separated along the second CCA axis, with Pomatoschistus minutus correlated with Up2m and Pomatoschistus microps correlated with pH. MM clustered in the center of the plot indicating a high/absent correlation with all variables. MS were separated along the two CCA axes (Figure 7b), with the majority clustering in the positive part of the first CCA axis and the negative part of the second CCA axis, correlated with temperature, chlorophyll a, TPM, and pH. A group of species composed of A. tobianus, Callionymus lyra, and Taurulus bubalis clustered in the negative part of the first CCA axis and were positively correlated with river flow. Commercial species tended to cluster on the positive part of the first CCA axis and the negative part of the second CCA axis, positively correlated with temperature, chlorophyll a, TPM, and pH. Labridae n.i. was the only commercial taxon clustering on the positive part of the second CCA axis and correlated with Up2m (Figure 7c). Table 4. Inter-set correlations of environmental variables with the first two CCA axes. Based on (square-root) transformed abundance of fish larval assemblages of the PNLN along the different sampling times. bold * Inter-set ≥ 0.4 corresponds to biologically important variables.

Discussion
This study aimed to investigate the larval fish assemblages of the PNLN to help understand the role of this coastal MPA for fish, with a special focus on commercially important species. The two most abundant species, P. gattorugine and A. tobianus, were non-commercial species, but the third most abundant taxa, Clupeidae n.i., was a commercial taxon. Although Clupeidae n.i. larvae were only identified at a family level, considering other studies performed along the NW Portuguese coast [31,45,66,67], we can infer that most of these larvae are Sardina pilchardus. S. pilchardus spawns and recruits along the Iberian coast, especially in the NW area [68], with spawning usually occurring from late autumn to late winter, although it may occur throughout the year [68][69][70][71]. In our study, the peak in Clupeidae abundance occurred during summer, indicating a spring-

Discussion
This study aimed to investigate the larval fish assemblages of the PNLN to help understand the role of this coastal MPA for fish, with a special focus on commercially important species. The two most abundant species, P. gattorugine and A. tobianus, were noncommercial species, but the third most abundant taxa, Clupeidae n.i., was a commercial taxon. Although Clupeidae n.i. larvae were only identified at a family level, considering other studies performed along the NW Portuguese coast [31,45,66,67], we can infer that most of these larvae are Sardina pilchardus. S. pilchardus spawns and recruits along the Iberian coast, especially in the NW area [68], with spawning usually occurring from late autumn to late winter, although it may occur throughout the year [68][69][70][71]. In our study, the peak in Clupeidae abundance occurred during summer, indicating a spring-summer spawning period since most Clupeidae larvae were recently hatched larvae aged 1-2 weeks old [71]. To support the hypothesis of PNLN being used by S. pilchardus as a spawning area, we compared the peak abundance observed in PNLN with other areas, and values were generally lower than those reported for other coastal protected areas [14,72], and for non-protected coastal areas both oceanic [73][74][75] and estuarine areas [46,71]. In some studies, peak values of S. pilchardus abundance were similar to those observed during this study [67,76,77]. However, these abundance comparisons must be tempered by the different sampling methodologies. In fact, the depth of ichthyoplankton collection is an example of a methodological approach with relevant influence on S. pilchardus larvae abundance since this species performs diel vertical migrations occurring at the surface during the night and sinking during the day [78]. Considering the current overexploited status of the S. pilchardus Iberian stock [79], the present results highlight the importance that coastal MPAs can play in the early life stages of economically important species.
Another interesting result was the fact that all marine migrant species observed in the PNLN are commercially exploited species, highlighting the important role of PNLN for marine species that use estuaries as nursery habitats. Indeed, PNLN may function as a connectivity pathway [29,30] between the oceanic areas and the adjacent Cávado and Lima estuaries. Marine migrant species have complex life cycles, with early life stages as larval and/or juveniles migrating from marine environments to estuarine zones, where they stay until joining adult populations in the ocean [31] and references therein. Fisheries pose additional stress to these vulnerable populations, and connectivity between different habitats is crucial to support the reestablishment of adult fish stocks [29][30][31]. Overall results indicate that the PNLN is an essential fish habitat for several commercially exploited species. Hence, this first study of PNLN larval fish assemblages indicates that the PNLN is an essential fish habitat for several commercially exploited species.
The species composition of the PNLN larval fish assemblages was similar to those reported for other coastal habitats of the Atlantic Iberian coast, which are typically composed of representatives of Blennidae, Gobbiidae, Clupeidae, Sparidae, Ammodytidae, and Scombridae families [31,45,72,73,75,80]. In addition, in other MPA areas with a shallow rocky substrate similar to those present in PNLN, the ichthyoplankton abundance and diversity tends to increase during the spring months and are typically dominated by one or two species of fish larvae from bottom dweller species (e.g., Parablennius pilicornis, Pomatoschistus pictus) [15,72,80]. This is similar to what was observed in PNLN since both peaks in abundance occurred during spring and were of P. gattorugine, a reef-associated bottom-dwelling species. In a non-protected coastal area close to the PNLN, fish larvae composition and temporal dynamic were also similar to what was observed in this study's results, with higher abundance and diversity during the spring months and fish larval assemblages being dominated by P. gattorugine and A. tobianus [45].
Larval stages of fishes are highly vulnerable to environmental parameters, both biotic and abiotic [24,81,82]. CCA analysis showed that temperature was an important driver for most of the species. In fact, temperature is usually regarded as one of the most important abiotic parameters for successful larval fish survival since it affects prey availability and also increases the levels of activity and metabolism [24,35,75,83]. During the study period, temperature exhibited the typical temporal patterns of temperate regions, with higher values during spring and summer and lower values during winter. Water column temperature also presented a typical temporal pattern with clear stratification during spring and summer months, with a high temperature near the surface decreasing with depth. During winter, the inverse was observed as temperature increased with depth and was slightly vertically stratified. Total particulate matter and pH were also identified by CCA as important environmental controls of fish larvae. These variables also revealed a clear temporal pattern, increasing during spring-summer when fish larvae were more abundant and diverse. TPM concentration may be related to estuary discharges located above the PNLN; these discharges have a southbound direction derived from the currents present in the region [41,42]. TPMs positive correlation with some fish larval species may also be related to chemical sensory cues used by fish larvae to find suitable nursery habitats and by adult fish species to find optimal spawning sites [84][85][86]. Although pH is a conservative variable, results showed that pH increased in the spring-summer period. This might have been associated with the occurrence of phytoplankton blooms, which can lead to changes in the pH of the water due to the release of oxygen or the uptake of carbon dioxide, and also with the decrease of river inputs from adjacent Lima and Cávado Rivers, whose water have lower pH [87].
River flow and PO 4 were also correlated with a group of species (A. tobianus, C. lyra, and T. bubalis). River flow values were higher during winter months when these species were more abundant. PO 4 values were inside the ranges that are usually observed in this coastal region [88], with the exception of Autumn 2017 and Winter 2018, which showed higher PO 4 values (max value = 1.71 µM L −1 ). Such PO 4 high values may be related to agricultural runoff since, despite its status as an MPA, there are extensive agricultural areas along the PNLN terrestrial area [87]. Nonetheless, the correlations associated with river flow and PO 4 may be related to species-specific adult reproductive strategies, such as the case of A. tobianus, which usually spawns during winter months in this region [45].
Biotic factors such as zooplankton abundance and chlorophyll a may be related to peaks in larval abundance [42,89,90]. CCA results did not identify zooplankton abundance as an important environmental control of the PNLN larval fish assemblages, which could be explained by the sampling methodology of these organisms that were only sampled at the middle sampling points of each transect, leaving several sampling points without accurate zooplanktonic information. However, several species were correlated with chlorophyll a, most of them commercial species (e.g., Clupeidae n.i.) [83]. Apart from these factors, the survival of fish larvae can also be associated with other factors such as larval behavior (i.e., passive or active swimming along coastal tides), spawning location, and predation that play an important role in survival and composition of larval fish assemblages [44,77,[91][92][93][94].
Given the proximity to the Cávado River estuary, it was expected to observe the influence of this ecosystem not only in the abundance and diversity of the larval fish assemblages but also in the water characteristics, mainly at the Center transect that is close to the Cávado River mouth. However, there was no spatial pattern for almost all environmental parameters studied, not even during winter months when river flow reached the highest values, and consequently, the river plume of the Cávado estuary could be more evident. No spatial pattern was observed in the abundance or diversity of the larval fish assemblages that were overall homogenous throughout the study area, except in Spring 2018, when abundance was higher at the central transect. These results indicate that connectivity between the coastal area and the Cávado estuary can be compromised, with important negative impacts for marine species that use estuaries as feeding or nursery areas. In fact, currently, the river mouth of the Cávado estuary is almost closed due to sediment deposition and restricted to a tiny channel [87], which may compromise oceanriver connectivity.
The PNLN showed clear temporal variation in both larval fish assemblages and environmental conditions, with clear temporal changes in species abundance; one of the peaks in abundance was related to a commercially exploited taxa (i.e., Clupeidae n.i.); however, the other peaks observed were for non-commercial species. Almost all the environmental parameters had a temporal variation; temperature, as expected, had a significant correlation with most commercially exploited species as well as chlorophyll a. When comparing the larval fish assemblage results with other similar rocky substrate areas, with the same level of protection or with no protection, the PNLN presented similar fish larvae composition and abundance. These results appear to indicate that the PNLN fish larval assemblages are not influenced by the conservation status as MPA, although results also highlighted its importance for some vulnerable commercially exploited species such as S. pilchardus. Given the apparent importance of the PNLN as a spawning area and/or nursery habitat for commercially important species, it would be useful to enforce conservation measures of PNLN and extend the PNLN limits, especially northwards, in a way that would include the Lima estuary; this would ensure the protection of this important fish nursery while also increasing this MPA size which has been highlighted as an effective measure to increase the protective role of MPAs [2,11]. Finally, the increased sand sedimentation of the Cávado River mouth should also receive special attention from PNLN managers to ensure that connectivity of early life stages of fish, especially those that use estuaries as nursery areas, are not compromised.

Conclusions
The present study showed that PNLN larval fish assemblages were mostly composed of reef-associated species (e.g., Blenniidae) and marine commercially exploited species that use estuaries as nurseries areas, indicating that PNLN is an essential fish habitat for marine fish species. Larval fish assemblages were characterized by a strong temporal pattern associated with seasonality of environmental parameters such as temperature, pH, PO 4 , chlorophyll a, TPM, and Cávado River flow. In contrast, no clear spatial pattern was evident in both environmental parameters and fish larvae indicators, revealing a spatial homogeneity throughout the PNLN area. The present study contributed with relevant baseline data knowledge of PNLNs role for fish, namely for the early life stages, that can help future management plans to improve the MPA efficiency to buffer negative impacts of fisheries and improve the connectivity to early life stages of commercially important species.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.