Table 1 displays the estimation outputs, where the significance level is *α* = 5% (as default). Almost all predictors are significant (except 'visibility'), displaying a *p*-value close to zero and lower than *α*, and the residuals are null, highlighting the robustness and validity of the results. Looking closely at the estimates, two findings should be dealt with carefully: (i) according to the predictor 'visibility', it would not be significant across units over time even if its related sum of squared dissimilarities and Pseudo-F statistic were null and highly large, respectively; (ii) the predictor 'oktas' displays, on average, larger *p*-values and lower Pseudo-F than the other significant ones. These findings highlight that multicollinearity problems would matter. For instance, the factor 'visibility' could be strongly correlated with one or more predictors within the model, rejecting the assumption in (ii) (exchangeable variables, such as sites, observations, and factors, are independent). To deal with this, we ran a new PERMANOVA test dropping 'oktas', the factor displaying lower Pseudo F statistic, and obtained the expected result (Table 2). The variable 'visibility' becomes significant with lower SS and higher Pseudo-F, compared to before. According to the other predictors, the estimated results improve as well. Thus, the estimates are now valid and unbiased.

**Table 2.** PERMANOVA analysis between the elements 'years' and 'effects' across units, dealing with multicollinearity problems. Here, the labels stand for 'degrees of freedom' (df), 'sum of squared dissimilarities' (SS), 'Pseudo F statistic' (Pseudo-F), and 'associated *p*-value' (Pr (>F)). The significance level is: (\*\*\*) significance at 1%.

The five main considerations are, in order: (i) there is a relevant heterogeneity among environmental effects. Thus, they are time-varying and change significantly year by year; (ii) heterogeneity also matters across units, and thus these effects heterogeneously affect sharks' prey preferences over time; (iii) non-linear relationships between predictors and outcomes need to be quantified through appropriate econometric algorithms; (iv) the environmental factors all have a strong effect on WSs' prey preferences; (v) potential environmental problems can be addressed in a particular time. For instance, we are interested to investigate why WSs would prefer the tuna bait rather than the seal-shaped decoy in 2008 and 2011. The last two points are addressed in the second and third step, respectively.

#### *3.3. LOGIT Regression Model*

The second step, the LOGIT regression model (Table 3), highlights the main environmental factors causing the WSs to prefer the seal-shaped decoy instead of the tuna bait.

Table 3 displays the estimation outputs, where the significance level is *α* = 5% (as a default). The five main findings are, in order: (i) all the results displayed in the first step are significant; (ii) all predictors, except 'temperature', positively affect the variable of interest 'prey'. For instance, better sea conditions increase the seal-shaped decoy preferences across sightings. The same accounts for high tide and better underwater visibility. On the contrary, higher water temperature increases the tuna bait preferences across sightings; (iii) according to the magnitude of the model, the factor that most affects the outcome is 'tide' (2.10), followed by 'visibility' (0.83), 'temperature' (−0.66), and 'sea' (0.65). The constant has not been considered, displaying an effect on the outcome, with all predictors being fixed (*Xi* = 0); (iv) 'visibility' and 'oktas' predictors are strongly negatively correlated (−43%) and show a significant linear relationship (Figure 2). Thus, by dropping 'oktas' within the system, any (potential) collinearity problems are dealt with, and we obtain better predicted outcomes; (v) the variable 'year' is also significant and positive, highlighting

the presence of relevant endogeneity issues affecting sharks' predatory behavior and heterogeneous effects given the unexpected environmental effects.

**Table 3.** Logit regression model accounting for all the five environmental factors across units over time. Here, 'Coefficients' refers to the factors within the model; 'Estimate' denotes *bk* (the estimated *βk*); 'SE' stands for standard error; 'z-value' denotes the test statistic obtained for each predictor; and 'Pr (>|z|)' refers to the associated *p*-value in a two-sided hypothesis test (where the null accounts for non-significance). The significance levels are: (\*\*) significance at 5%; and (\*\*\*) significance at 1%.

**Figure 2.** Scatter plot by regressing 'okta' (independent variable) on 'visibility' (dependent variable). The estimation output is negative and significant, displaying a consistent relationship between the variables. The scatter plot shows a large component of outliers affecting the output due to omitted variables and unobserved heterogeneity (endogeneity issues) dealt with the logistic function. The significance level used in the analysis is *α* = 5% (as default).

#### *3.4. Sample Marginal Effects*

The third and final step, the sample marginal effects (Table 4), evaluates major environmental factors affecting feeding behavior.

**Table 4.** Sample marginal effects for each observation unit, given *n* observations, are accounted for. Here, 'Coefficients' refers to the factors within the model; 'dF/dx' denotes the partial derivatives displaying the marginal effects of the predictors (*x ik*) on *yi* ('prey'); 'SE' stands for standard error; 'z-value' denotes the test statistic obtained for each predictor; and 'Pr (>|z|)' refers to the associated *p*-value in a two-sided hypothesis test (where the null accounts for non-significance). The significance level is: (\*\*\*) significance at 1%.

Table 4 displays the estimation outputs, where the significance level is *α* = 5% (as default). Three relevant findings are addressed: (i) the main factors affecting the outcome correspond to the ones found in the logistic function (second step). More precisely, in order of importance, we have: (1) 'tide', (2) 'visibility', (3) 'temperature', and (4) 'sea'; (ii) the marginal effect values appear sensible. For instance, a one-unit change associated with an observation increases the probability of preferring the seal-shaped decoy by 61% with high 'tide', 55% with better 'visibility', 51% with less water temperature ('temp.'), and 49% with better 'sea' conditions; (iii) the factor 'year', though it affects the outcome with lower associated probability, is significant and must be included in the analysis when studying environmental problems.

Finally, we look more deeply into these sample marginal effects, referring to the average individual effects to investigate the same effects in every period (Table 5). We focus on 2008 and 2011, representing the time periods in which white sharks preferred the tuna bait rather than the seal-shaped decoy [10]. In 2011, according to the total sightings, tide events were exclusively characterized by low tide and the highest water temperature was registered. The same results were found in 2013, with two key differences: (i) in 2011, the total sightings were more than double (63 vs. 29), highlighting the highest sample marginal effects according to the predictor 'tide'; (ii) the water temperature (the third most significant outcome-affecting factor) was, on average, balanced in 2013 compared to the other time periods, whereas in 2011 the temperature held at 18 ◦C. In 2008, the most sightings were characterized by low tide (25 on 40). Similar effects were observed in 2009, showing the same visibility and water temperature. However, in 2009, tide ranges were similar (12 observations occurred during low tide and 10 during high tide out of 22 total sightings), showing that, once again, 'tide' denotes the most significant/important factor in terms of affecting sharks' predatory behavior.

**Table 5.** Summary and descriptive table for the average individual (total sightings) over time. Here, 'Tot. Sigh.' stands for total sightings, and 'W. Temp.', 'U. Visibility', and 'Sea Cond.' refer, in average ± SD, to water temperature (◦C), underwater visibility (m), and sea conditions, respectively. The ranges would correspond to the minimum and maximum values displayed in the table.

#### **4. Discussion**

This study demonstrates that prey discrimination of bait-attracted WSs in Gansbaai is affected by environmental factors that drive the sharks to detect and choose prey depending on the conditions. Passive preys (tuna bait and seal-shaped decoy) do not show spatial autocorrelation, as demonstrated by the Moran's Test, confirming that the odor from the tuna bait did not serve as chum within the 10 m range of the experiment. In Gansbaai, as stated by Micarelli et al. [10] from 2008–2013, the majority of WSs transient population, immatures included, were attracted by the seal-shaped decoy rather than by the tuna bait, except in the two years of 2008 and 2011. As observed by the first PERMANOVA test (Table 1), 'visibility' was not significant and 'oktas' displayed larger *p*-values and lower Pseudo F than other variables due to a multicollinearity problem, since underwater visibility is strongly correlated (−43%) with the light levels, as stated by the linear regression model in Figure 2. For that reason, light levels were omitted in the second PERMANOVA test (Table 2). As suggested by Strong [12], one of these aspects is a consequence of the other: low light levels (represented by high oktas values) reduce the surface underwater visibility that disadvantages the seal's ability to detect the WS coming from below.

We observed a relevant heterogeneity among environmental effects that undergoes large changes year by year and affect sharks' prey choice over the time. According to the sample marginal effects, the factors influencing the prey discrimination and the choice of the seal-shaped decoy for a one-unit change associated with the observation (Table 4) are, in decreasing order, tide ranges (61%), underwater visibility (55%), water temperature (51%), and sea conditions (49%).

Tide range is the most important environmental factor that influences prey choice. At Seal Island in South Africa, Hammerschlag et al. [2] stated that attack frequency to Cape fur seals increased significantly during high tides within 400 m of the island and over a depth range of 5–31 m where seals swim over shallow reefs. During high tide, the minimum depth from which sharks could approach seals was closer to shore, increasing the availability of prey and the possibility of predation. Additionally, at the Farallon Islands in California, WSs' attacks increased with high tide that reduced the haul-out area for elephant seals, forcing them into the water and increasing predation events [11,16,27]. In fact, during high tide events, the number of Cape fur seals around our boat anchored at 100–150 m from Dyer Island increased, inducing sharks to use their visual acuity both to attack them and our seal-shaped decoy. For that reason, we propose that high tide could induce WSs from Gansbaai to turn their attention towards seal-shaped decoys rather than towards the tuna bait, consequently obtaining more caloric sources over low-energy ones, according to the theory of optimal foraging [28,29].

It is true that most of the years (2008, 2009, 2011, and 2013) were characterized by low tide, but it is also necessary to observe other environmental parameters that influence prey choice, such as underwater visibility and water temperature.

Underwater visibility is the second environmental factor that influences the sharks' discrimination of prey and, on average, was good during all the years (at a visibility range registered from 1 to 5 m during all the years, minimum average visibility was 2.55 m in 2008 and 2009; maximum average visibility was 3.55 m in 2013). Higher average visibility and, consequently, low oktas values seemed to facilitate the attack on the seal-shaped decoy in Gansbaai. Our findings do not support the Seal Island results of Hammerschlag et al. [2] that observed a higher predatory success of WSs on Cape fur seals during low light levels and low underwater visibility. This inconsistency could be linked to two reasons: (i) stalking and ambush with low visibility are important for successful prey capture in many fishes, as well as sharks, including *C. carcharias* [12,30–32]. However, a depth range of 26–30 m around Seal Island is optimal for sharks to remain undetected and stalk seals from below with low visibility and enough vertical distance to build the perfect moment for attack and surface strike, as observed by Hammerschlag et al. [2]; on the other hand, around the Dyer Island Nature Reserve and in the Shark Alley channel that divides Dyer Island from Geyser Rock, the maximum depth is ca. 8–10 m and the vertical distance is not enough to surprise seals by exploiting the low visibility. For that reason, most of the attacks to passive preys at Gansbaai were horizontal [15] and in the presence of a good average underwater visibility. Sometimes, as observed by Huveeners et al. [17] in Neptune Island (South Australia), sharks hunting near the surface orient the sun directly behind them during approach, directly illuminating prey and increasing their detectability. (ii) As demonstrated by Micarelli et al. [10], vision plays an important role in mature and immature sharks with a dietary shift in their feeding pattern, and a good average visibility could help them to more consistently secure higher-energy prey sources, such as our shaped decoy.

Water temperature is the third environmental factor that influences prey selection and plays an important role in the physiology of ectotherms, regulating their internal functions, distribution, foraging, reproduction, prey distribution, and abundance [14]. Unlike most sharks and rays, the white shark is an endothermic animal within the family Lamnidae, and water temperature has an influence on its prey distribution rather than on the shark's physiology [19]. Warmer water resulting in blooms of diatoms is associated with an increase in the number and diversity of teleosts, which are prey of white sharks [14]. Different tuna species, used as passive prey to attract white sharks, are present in South African waters

and around Gansbaai, such as Albacore or Longfin tuna *Thunnus alalunga* (Bonnaterre, 1788), Yellowfin tuna *Thunnus albacares* (Bonnaterre, 1788), Big-Eye tuna *Thunnus obesus* (Lowe, 1839)*,* and Skipjack tuna *Katswonus pelamis* (Linnaeus, 1758) [33]. All of these species, except the Albacore, are mostly found off of the Southern Cape Coast of South Africa from spring to summer (www.deep-sea-fishing.co.za, accessed on 15 April 2022) when waters are warmer; their higher availability and abundance during these seasons and during our observation period in which the water was warmer likely influenced the tuna bait choice. Although WSs generally prefer to feed on juvenile Cape fur seals that leave Geyser Rock to forage offshore for the first time during the autumn and winter season, warmer water induces phytoplankton blooms that reduce the underwater visibility, hindering the shark's ability to detect the odorless seal-shaped decoy and encouraging the tuna bait choice.

According to what has been said, warmer water temperatures registered in 2011 (18 ◦C on average) were associated with only low tide events (63), one of the lowest average underwater visibilities (2.79 m), and the highest number of interactions (63) that could have influenced the choice of the tuna bait. Additionally, in 2008 sharks preferred the tuna bait, even if the average water temperatures were the lowest registered (13.25 ◦C), but most of the interactions were characterized by low tide events (25 out of 40 total interactions) and by the lowest average underwater visibility (2.55 m), forcing the sharks to choose the tuna bait. During the other years (2009, 2010, 2012, and 2013) sharks preferred the seal-shaped decoy. In 2009, tide ranges were similar between high and low (10 vs. 12, respectively). The average water temperature in 2009 was also slightly higher than 2008 (14 ◦C vs. 13.25 ◦C, respectively) despite having the lowest average underwater visibility (2.55 m), justifying the choice, since both total sightings were the lowest registered during the entire study (22) and tide range remained the main influencing factor. In 2010, mostly high tide events were recorded (35 out of 54 total interactions) with a cool average water temperature (15 ◦C) and good average underwater visibility (2.83 m), justifying the seal-shaped decoy choice. In 2012, the situation was like 2009, since the number of low and high tide events were identical (16); the average water temperature was cool (15.2 ◦C), the underwater average visibility was good (3 m), and total sightings were few (32). In 2013, high tide events were not observed, the average water temperature was cool (15 ◦C) and total sightings were few (29), but average underwater visibility was the best out of all the years (3.55 m), probably justifying the seal-shaped decoy choice.

Ultimately, sea condition is the variable that affects prey selection the least. In any case, sea condition was, on average, slightly rough during all the years, except for 2013 and 2011 where it was calm. Wind intensity has a direct effect on the creation of rough conditions, and, therefore, could propel pinnipeds' chemical stimuli (e.g., excreta, blood, organic compounds) from the Dyer Island Nature Reserve to the open sea, helping the sharks to locate the area where the seals are concentrated, as observed in the study of Strong [12] and Hammerschlag et al. [2] for Seal Island. At the same time, pinnipeds are forced to swim against the current with rough sea, producing sounds and losing their ability to maintain subsurface vigilance, which gives the sharks greater possibilities to strike [12]. However, it is also true that we used an odorless seal-shaped decoy to attract the sharks. A calm or slightly rough sea also has a direct effect on underwater visibility that, in our study, is the second most important environmental factor to influence prey discrimination, helping the sharks in Gansbaai to be able to visually detect the seal decoy more successfully. Notably, 2011 saw an average calm sea condition (1.94), but the other environmental factors, such as tide ranges events, underwater visibility, and water temperature, had a greater influence on the choice of the prey, which, in this case, was tuna. On the other hand, we registered the roughest average sea condition in 2008 (2.45 in terms of mean), alongside the lowest underwater visibility, the lowest water temperature, and a higher number of low tide events, justifying the tuna bait selection.

#### **5. Conclusions**

Environmental factors influence the WSs' sensory ecology to choose between passive prey types during interactions. The choice of the seal-shaped decoy is strictly related to high tide events, good underwater visibility, cooler water temperature, and calm or slightly rough sea. The last three conditions highlight the importance of the mature and immature WSs' visual ability in detecting the prey. WSs are initially attracted by the olfactory trace emitted by the tuna bait, but they subsequently shift their focus to the more calorically rich seal prey, helped by good environmental conditions. High tide, on the other hand, gives the sharks more chances to get closer to the island and reduces the prey haul-out area, increasing the possibility of predation on seals and the seal decoys alike. Some environmental factors were not taken into consideration during the study, such as wind direction, currents, and distance from the island, and could, therefore, be considered as topics of future inquiry at this site.

**Author Contributions:** Conceptualization, F.R.R., P.M. and E.S.; methodology, F.R.R., P.M. and A.P.; software, G.G. and A.P.; validation, P.M., F.R.R., E.S., A.P. and M.M.; formal analysis, A.P., P.M. and F.R.R.; investigation, P.M., F.R.R. and E.S.; data curation, F.R.R., P.M. and A.P.; writing—original draft preparation, F.R.R., M.M. and P.M.; writing—review and editing, F.R.R., P.M., E.S., A.P., G.G. and M.M.; funding acquisition, P.M. and E.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** The study and experimentation protocols were reviewed and approved in accordance 166 with the Directive 2010/63/EU.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** https://www.researchgate.net/project/Great-White-Shark-Carcharodoncarcharias-Behaviour-Ecology-and-cotoxicology/update/61312b952897145fbd6df39c accessed on 10 April 2022.

**Acknowledgments:** We are grateful to the CSS team members that carried out 11 expeditions for their indirect financial support of this research, and thanks are also due to Shark Diving Unlimited for the logistical assistance and all other field assistance with data collection. All research protocols were approved by the University of Calabria and Siena.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

### **Courtship and Reproduction of the Whitetip Reef Shark** *Triaenodon obesus* **(Carcharhiniformes: Carcharhinidae) in an Ex Situ Environment, with a Description of the Late Embryonic Developmental Stage**

**Sérgio Ricardo Santos 1,2,\*, Veronica Takatsuka 1, Shayra P. Bonatelli 3, Nicole L. L. Amaral 3, Matheus F. Goés <sup>1</sup> and Rafael F. Valle <sup>1</sup>**

**Simple Summary:** The reproduction of key reef species is still largely unknown due to difficulties in documenting all elements and steps involved. Sharks are particularly affected by this scarcity of information due to being long-lived species, and witnessing courtship, gestation, and birth is still mostly limited to fortuitous encounters by divers or specimens captured by fishers. Still scarcely described in the literature, our study reports the successful reproduction of *Triaenodon obesus* in an ex situ environment, which offers an opportunity to observe all steps of the reproduction in detail. Furthermore, we offer the first description of the late embryonic developmental stage based on ultrasound imagery.

**Abstract:** Elasmobranchs represent a group of species under considerable anthropic pressure because of the scale of industrial and artisanal fisheries and the loss of essential areas for nursery and feeding, which are causing substantial population losses around the world. Reproduction in an ex situ environment enables a healthy population to be built and maintained in networks of public aquariums, increasing our knowledge of elasmobranch reproductive biology and offering the opportunity for reintroductions in areas where native populations have been removed. The study reports two successful pregnancies of the whitetip reef shark *Triaenodon obesus*, considered a vulnerable species by the International Union for the Conservation of Nature. Copulation and gestation data are provided, including ultrasound recordings of the late stage of embryo development. Ultrasonography was performed with the GE Logiq and convex transducer and revealed a fetus with defined fins and organogenesis, with definition of eyes, gills, liver, a heart with individualized chambers, partially defined kidneys, and a well-defined spiral intestine. A cartilaginous skeleton forming a posterior acoustic shadow was detailed, as well as a moving fetus with a biparietal diameter of 6.47 cm and a heart rate of 62 Beats Per Minute on spectral Doppler. This is the first successful reproduction of *T. obesus* in an aquarium in Brazil.

**Keywords:** ultrasound; elasmobranch; gestational; aquarium; conservation

#### **1. Introduction**

Declines in elasmobranch populations related to diverse and continuous anthropic impacts have been widely reported in the ecosystems of several parts of the world [1–4]. These losses represent only a fraction of the actual damage suffered by populations in these ecosystems, where elasmobranchs represent important components of the food web [5–8]. Fishery

**Citation:** Santos, S.R.; Takatsuka, V.; Bonatelli, S.P.; Amaral, N.L.L.; Goés, M.F.; Valle, R.F. Courtship and Reproduction of the Whitetip Reef Shark *Triaenodon obesus* (Carcharhiniformes: Carcharhinidae) in an Ex Situ Environment, with a Description of the Late Embryonic Developmental Stage. *Animals* **2022**, *12*, 3291. https://doi.org/10.3390/ ani12233291

Academic Editor: Martina Francesca Marongiu

Received: 30 October 2022 Accepted: 22 November 2022 Published: 25 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

impacts are also underestimated, given the scarcity of data for most marine ecosystems, particularly in tropical waters and off the coasts of developing countries, with numerous unmonitored fishing communities associated with small-scale fisheries [9]. The growing popular appeal of the conservation of marine environments, the public demand to establish the extent of environmental changes, and the worsening of the conservation status of several flagship species of elasmobranchs have promoted the growth of sustainable management plans in the last decades. A compromise between governments and nongovernmental organizations (NGOs) was established based on commitments assumed since the Rio 1992 convention [10,11], yet the environmental policies and the species protection efforts implemented in the following decades have still been unable to revert or mitigate the population declines of several species. Special interest has been devoted to the conservation of coral reefs and their inhabitants, considering the fragility of these ecosystems in the face of climate change and a wide variety of human impacts, with the whitetip reef shark *Triaenodon obesus* (Rüppell, 1837) being a known representative.

Despite being a coastal demersal species with a limited dispersal capacity and home range, *T. obesus* is the most widely distributed reef shark species in the world, found from the Red Sea to Micronesia, as well as from the coast of Mexico to Panama in the Eastern Pacific, and its presence is suspected even in the southwestern Atlantic [12–14]. This wide distribution, however, does not mean that the species is not threatened. The continuous pressure exerted by human impacts on the ecosystems they inhabit, whether through direct activities such as fishing [15] or through environmental degradation and loss of habitat [16,17], exposes populations of this species to conditions hostile to their survival. The situation is worsening in the context of climate change and its effects on marine ecosystems, within which coral reefs are especially impacted areas [16,18]. As a direct consequence of this situation, the decline seen in the last three generations of whitetip reef shark populations led to the classification of their conservation status as 'vulnerable' by the Red List of the International Union for Conservation of Nature (IUCN) [19].

*Triaenodon obesus*, despite its wide distribution, is assessed as vulnerable due to population decline and habitat loss in many areas of its range [20,21], but gaps in knowledge persist. The whitetip reef shark is a placentotrophic viviparous species, with a matrotrophic nutrition of the fetus [22,23]. That is, *T. obesus* is a live-bearing shark with the development of the embryo sustained by nutrients produced by the mother and delivered by a placental connection, a feature only identified within five families of Carcharhiniformes (Carcharhinidae, Sphyrnidae, Hemigaleidae, Leptochariidae, and Triakidae) [22]. Recently, ultrasound studies on elasmobranchs in captivity were published, focusing mainly on determining gestation and general anatomy [24–27]. However, there is a scarcity of studies focused on the description of the ultrasonographic aspects of fetal organs and their organogenesis in sharks. Since *T. obesus* is considered vulnerable, the detailed monitoring of pregnancy is of paramount importance for the preservation of the species.

While the high cost of field studies hinders the in situ observation of the species, alternative sources are available. The limited home range, average size of adult specimens, and ease of access and exchange among public aquariums have made the species a common sight in public exhibitions around the world [28]. The need for data that support the handling of *T. obesus* under human care and the knowledge about the biology of the species that can be acquired in public aquariums and research centers, especially by monitoring aspects that are difficult to track in its natural environment, such as reproduction, make ex situ studies fundamental in reversing its population decline. The present study describes the reproductive success at the Rio de Janeiro Marine Aquarium (AquaRio) in 2022 with the birth of a male and a female specimen of whitetip reef shark, as well as presenting continuous monitoring of the courtship, gestation, parturition, and postpartum processes. Due to the scarcity of ultrasound data of pregnancy for the species, the objective of the work also included a description of the fetal structures in the final stage of gestation.

#### **2. Materials and Methods**

The Rio de Janeiro Marine Aquarium's breeding stock consists of more than 2600 specimens of 306 species of marine and freshwater organisms, with an emphasis on 67 individuals of 11 species of Elasmobranchii. The whitetip reef shark *T. obesus* was originally represented by five specimens (three females, two males) from the coast of Indonesia that were transported to Brazil in August 2017. At the time, all specimens were smaller than 90 cm TL (total length), which identified all individuals as juveniles, considering the size at first maturation reported in the literature (♂= 104 cm TL, ♀= 105 cm TL) [29]. After a period of quarantine, the sharks were introduced into the 3.5 million-liter marine exhibition tank, with a maximum depth of 7 m, where they have been kept for the past four years. The oceanic tank is equipped with a mechanical filtration system comprising 20 sand filters, 9 skimmers (Altamar 25c) with ozone injection (O3R and Altamar, 40 g/h), and a degasser (1200 gal/h). Reverse bottom biological filtration is complemented by weekly siphoning. Temperature control is achieved via a heat exchanger. The diet of whitetip reef shark adults is offered daily and includes *Euthynnus alletteratus* (Rafinesque, 1810), *Coryphaena hippurus* Linnaeus, 1758, *Sardinella brasiliensis* (Steindachner, 1879), *Cynoscion* spp., *Anchoa* spp., *Pseudupeneus maculatus* (Bloch, 1793), and *Dorytheutis* spp. Feeding is unrestricted three times a week at four points in the oceanic tank where specimens of *T. obesus* have been conditioned to go.

Two observation strategies were used on our protocol: The first one comprises two programmed daily visits on the entire perimeter of the oceanic tank performed by a veterinarian or an aquarist. The ethogram records any deviation in behavior from the known baseline for each elasmobranch. The second source of observations derives from opportunistic records (i.e., video, photographs, or reports) from environmental educators, a group comprised of a dozen professionals that are continuously present around the oceanic tank for ten hours a day. When isolated, the pregnant female is closely monitored throughout the day by a team of aquarists and veterinarians.

#### *Ultrasound Examination*

The ultrasound examination of the first parturient whitetip reef shark female was performed using a Logiq e model (General Electric®, Boston, MA, USA) equipped with a convex transducer ranging in frequency from 3 to 5 MHz. Latex gloves coated with a generous amount of gel were used to protect the transducer from possible damage caused by contact with salt water and shark skin and to promote adequate transmission of the ultrasound wave. The animal was physically restrained using the tonic immobility maneuver and positioned in dorsal decubitus throughout the examination, according to the structure being analyzed and with the patient's cooperation. A state also known as animal hypnosis, it induces loss of muscle tone and equilibrium while the shark remains unresponsive to major stimulation. It allows safer and less stressful procedures and, unlike chemical anesthetics, allows immediate recovery and minimal disruption to respiration [30,31]. An ultrasound scan of the entire coelomic cavity, the diagnosis of pregnancy, and the detailing of the fetus were performed. Anatomical descriptions follow Crow and Brock [32] and Mylniczenko [23].

#### **3. Results**

#### *3.1. Courtship, Copulation, and Pregnancy*

The presence of bite scars on the female's fins was noticed in the first days of June 2021, and the first courtship with copulation was recorded on 2 June 2021. The event lasted for 3 min 55 s from the moment of the bite until the removal of the clasper and separation of participants (Video S1). The courtship began with the male's pursuit of the female until the male was anchored by biting the female's left pectoral fin, starting a spiral descent to the bottom of the tank with slow or no swimming. Still in the water column, the male inserted its right clasper into the female's cloaca, remaining attached until the end of copulation. The position of the couple varied throughout the process, with the snout keeping in contact

with the bottom and the body angle varying between 90◦ and 45◦ with the substratum. The swimming of both stopped, and the male moved its pelvic region toward the female. Copulation was closely observed by the second male of the squad and by another female, and physical contact with the mating couple was recorded, removing them from their initial position in the water column and both coming to lay their bodies directly over the substrate. Although they suffered this interference, the second male did not bite the female, nor was there any attempt to insert its clasper. Upon completion of copulation, the male and female did not immediately separate. The clasper remained inserted in the cloaca while both began to swim in opposite directions, forcing the withdrawal of the clasper and the couple's separation.

The monitoring of the behavior of the specimens of *T. obesus* in the oceanic tank also allowed observation of the female's escape strategies from unaccepted copulations. Several events were observed in the ocean tank, two of which were recorded on video and are described here (Video S2). The beginning of courtship followed the expected course and started with the chase by the male. The first strategy was to escape by fast swimming, avoiding the direct approach to the exposed pectoral fin and the necessary bite to anchor. If anchoring was successful, the female moved erratically and swam at different speeds while the male stopped swimming. Finally, the female made use of the substrate's complexity to release the male's bite, initiating a quick escape soon afterwards before the male could restart its approach.

#### *3.2. Internal Morphology*

Upon ultrasound examination, it was possible to detail a fetus with a well-defined brain, tail, and pectoral and dorsal fins, with well-defined gills and coelomic cavity (Figures 1, 2a–d, and 3a–d). The fetal liver was characterized by a hypoechoic, homogeneous parenchyma with regular contours, also showing the gallbladder with anechogenic content (Figure 2b). The kidneys were partially defined and presented a slightly heterogeneous hypoechogenic parenchyma. The intestinal spiral with the characteristic spiral pattern was evident, as was a defined parietal stratification (Figure 3a). The cartilaginous skeleton was hyperechoic and formed a moderate posterior acoustic shadow artifact (Figure 3b). It was possible to detail the cranial conformation, define the orbital fossa housing the anechoic eyeball with a hyperechoic capsule (Figure 2c), and evaluate the biparietal diameter, which measured approximately 6.47 cm (Figure 3c).

The vertebrae of the spine presented themselves as hyperechoic, individualized structures, homogeneously spaced, with an anechoic central area, possibly the spinal canal (Figure 3b). In the cross-section, the spinous processes of the vertebrae were also detailed. The gills were defined as parallel slits of alternating echogenicity (hyperechogenic and anechogenic) (Figure 2c), which, on color Doppler mapping, showed an evident signal, indicating the presence of fluid passage. The heart presented individualized chambers and an evident blood flow signal when we used the color Doppler tool (Figure 2d). With spectral Doppler, it was possible to measure the heart rate, which presented a value of 62 beats per minute at the time of the examination (Figure 3d). Considering that the average gestation period observed by Schaller [24] is 387 days and knowing that the neonate described above was born 12 days after the ultrasound examination, we calculated that it was approximately 375 days old. Therefore, the description of organogenesis by this method characterized the final third of the pregnancy.

**Figure 1.** Basic internal anatomy of a carcharhinid shark (*Carcharhinus melanopterus*), showing the location of principal organs. Adapted from Crow and Brock [32].

**Figure 2.** Ultrasonography of a female *Triaenodon obesus* (Rüppell, 1837) in late pregnancy: fetal orientation (brain/tail) and pectoral fins (**a**), hypoechogenic liver (**b**), eyes and gills (**c**), and heart with individualized chambers (**d**).

**Figure 3.** Ultrasonography of a female *Triaenodon obesus* (Rüppell, 1837) in late pregnancy: intestinal spiral and parietal stratification (**a**), cartilaginous skeleton, forming a posterior acoustic shadow artifact (**b**), biparietal diameter of 6.47 cm (**c**), and cardiac blood flow of 62 beats per minute (**d**).

#### *3.3. Birth and Feeding*

A *T. obesus* female of 145.0 cm TL and 20.1 kg TW (total weight) successfully completed her reproductive process with the birth on 7 March 2022 of a male of 58.0 cm TL and 1.15 kg TW. A second female of 136.4 cm TL and 19.0 kg TW gave birth to a female neonate of 56.4 cm TL and 1.0 kg TW on 7 May 2022. This was the first occurrence of the reproduction of a whitetip reef shark in an aquarium in Brazil. Both offspring were active, displaying constant swimming even though they rejected food for the first three and seven days, respectively. The first neonate fed for the first time on the third day, consuming 15 g of tilapia (*Oreochromis* sp.). The diet in the first trimester consisted of Brazilian sardine *Sardinella brasiliensis*, tilapia, manjubas *Anchoa* spp., shrimp *Penaeus* spp., and squid *Doryteuthis* spp. Over the initial three months, teleosts composed 76% of the diet, followed by shrimp (21%) and squid (3%). As for the second neonate born in 2022, feeding only started after seven days, and any food offered before then was rejected. The items offered in the diet were the same as those of the male, with tilapia accepted first, while in the following feedings, the male consumed mostly manjubas and sardines. The pregnant female isolated during the quarantine showed a marked reduction in her feeding during gestation, fasting for up to ten days, with erratic intervals in the following days until giving birth. Following veterinary recommendation and internal postpartum protocol, the recovering female was medicated via intramuscular injection with Ceftazidime 30 mg/kg and sodium methylprednisolone succinate 2 mg/kg, every 72 h, totaling five applications, due to a slight increase of leukocytes. Afterwards, she returned to her normal appetite and was transferred to the ocean tank.

#### **4. Discussion**

#### *4.1. Courtship, Copulation, and Pregnancy*

These records, made between 2021 and 2022, allowed us to evaluate the process of courtship with copulation and the gestation of the whitetip reef shark in an ex situ environment. The observations were compared with events already described in the literature and coincided with the copulations recorded previously. Whitney et al. [20] described the estimated period in the natural environment for the mating season and birth, based on photo identification made off the coast of Kona, Hawaii (USA). The authors pointed out a possible overlap of the birth period with the mating period, estimating that gestation spans one year. In the case reported here, the records of mating scars were made in June 2021, while births occurred between March and May, close to but shorter than the period reported for the species by Whitney et al. [20]. Comparing the mating record obtained in 2021 and the date of the first birth in 2022, the gestation period for the female born ex situ would have been 278 days, or just over 9 months. However, if the filmed copulation and the birth of the second neonate are considered, the gestation period would have been 339 days, which differs little from that observed for the species. Schaller [24] reported 5 births of *T. obesus* at the Steinhart Aquarium (San Francisco, CA, USA) between 2001 and 2004 and observed a longer gestation period, between 355 and 422 days, based on the first mating scar and the date of birth. Based on specimens captured off Indonesia by a commercial fishing fleet, White [33] observed two pregnant females of sizes close to our females (158.1 and 140.6 cm TL), with four formed embryos gestated by the first and three embryos found in the second, and estimated the total gestation time to be greater than six months.

Whitney et al. [34] observed three courtship events in the Coco Islands (Costa Rica), two of which involved only one male and one female. The beginning was recorded at mid-water, with the descent to the substrate provoked after the bite on the female's pectoral fin and, consequently, the inhibition of swimming. The courtship concluded with copulation in only one of the events. The courtship attempts we observed in our oceanic tank were all initiated in mid-water. Similar to the observations made by Whitney et al. [34], most attempts did not result in copulation. Incomplete copulation occurred owing to the female's detachment from the pectoral fin bite and the male and female's inability to keep swimming in parallel and in close enough contact to allow a second bite to restart the descent to the bottom. Despite being in an ex situ environment, the female's nonacceptance and escape strategies observed at AquaRio confirmed previous observations made in a natural environment [34] that indicated the success of copulation being largely dependent on acceptance from the female. Pratt Jr. and Carrier [35], in a review of reproductive behavior in elasmobranchs, noted that biting is apparently a universal characteristic of the group, although for larger species some degree of cooperation or acceptance by the female is necessary for actual copulation to occur. Copulation events observed or recorded within other elasmobranch species from our aquarium (*Ginglymostoma cirratum* and *Stegostoma tigrinum*) support this finding, especially for an event involving a female of *S. tigrinum*. In this observed event, after a tail bite, the female moved to the substrate and exposed its belly, even though the male did not conclude the copulation. In the case of *T. obesus*, the nonacceptance of courtship was associated with accelerated swimming, even after the pectoral bite, and erratic movement, which made it difficult for the male to chase. The cloaca was protected by the female's body movement and the structural complexity of the habitat was used to enact the escape or disarm the bite.

Tricas and Le Feuvre [36] described two mating events for the species from observations on Molokini Island, Hawaii (USA), in a shallow (7 m) protected reef area. The description began with the male anchored to the female by biting the pectoral fin, the

male's right clasper inserted into the cloaca, and the bodies positioned with the head downward and the rest of the body suspended in the water column at an angle of 45◦. The two copulation events observed in our oceanic tank coincided with the events described by Tricas and Le Feuvre [36], occurring shortly after the spiral descent. Whitney et al. [34] expanded on the findings of Tricas and Le Feuvre [36], having observed a larger set of courtships, copulations, and the functioning of the siphon bags. Whitney et al. [34] reported the presence and behavior of other males that did not participate directly in copulation but remained present and sometimes participated in the initial moment of copulation, biting the female's pectoral fin to provoke the descent to the bottom, but separating soon after without copulating. In the two courtships with copulation recorded by us, the second male was observed circling the main couple, not having participated with any bite on the pectoral fin. In the first event, in 2021, the male's physical contact caused the couple to leave their suspended and inclined position in the water column, but there was no insertion of the clasper by the second male, who moved away from the couple shortly afterwards. The most recent event did not present any interference from the second male other than constantly circling the couple during mating.

Studies on changes in behavior and feeding during pregnancy and in the period before birth are scarce in the literature for most species [37,38]. The whitetip reef shark is characterized by greater activity at night, remaining at rest for most of the day [39,40]. After monitoring the activity of adults by recording swimming acceleration and the intensity of breathing in neonates, Whitney et al. [39] pointed out that the latter assumed the pattern observed for adults soon after birth. The circumstances of the two neonates born in quarantine and in the ocean tank agreed with the observation made by Whitney et al. [39], with the same circadian pattern as adults being verified in neonates. Even so, the female born in quarantine and observed in more detail demonstrated more active diurnal swimming in the first day after birth before assuming the nocturnal pattern similar to that described for the species in the following days.

The feeding during the pregnancy of the first female of *T. obesus* was intermittent, being generally characterized by the acceptance of food every two days. It is important to note that the transfer process from the oceanic tank to the isolation in the quarantine tank was followed by a period of ten days in which the pregnant female remained within the expected activity pattern for an adult, but without accepting any food. Although we cannot assume this to be a feeding pattern of a pregnant whitetip reef shark, the disturbance of the transfer may have provoked enough stress to cause an interruption in the feeding of a pregnant female, despite it being a routine activity that does not cause the same reaction in nonpregnant females.

Given that it was the first pregnancy for the two females in question, in addition to the known relative sizes of the pregnant female and the embryo in *T. obesus*, the occurrence of only one offspring each could have been caused by the restricted space for their development. However, based on a visual census conducted in 21 reefs of the Great Barrier Reef (Australia) and records for 76 females, Robbins et al. [41] and Robbins [1] identified pregnancies of 1 to 4 offspring with an average of 2.07 neonates and showed no significant correlation with the size of the pregnant female or her age, with this being considered a particularly low level of fecundity compared to other species of Carcharhinidae.

#### *4.2. Internal Morphology*

Many of the studies devoted to the embryonic development of sharks have been carried out using oviparous species, given the ease of obtaining fertilized eggs to monitor the different stages of development [42]. In viviparous or ovoviviparous species, monitoring is more complex owing to restrictions on obtaining pregnant females found in the natural environment [43–46]. Research heavily relies on accidental captures by fishing of pregnant females at different stages of embryonic development to compose a sequence of embryonic stages [38,39]. Pregnancy monitoring in an ex situ environment offers the opportunity to monitor the same females from mating to birth, recording behavioral and

physiological changes throughout the reproductive cycle [24]. The monitoring, however, depends on routine exams that allow the identification of pregnant females at the beginning of embryonic development. Maintenance of the usual behavioral pattern and the minor morphological alteration of the female during the initial and intermediate stages of gestation did not call attention to our whitetip reef shark until the last months before birth. While ultrasound exams are dedicated to the diagnosis and treatment of individuals suspected to be ill, monitoring by routine scheduled exams is necessary for viviparous species of elasmobranchs to enable the identification of pregnant females and recording of the early stages of embryo development.

Sharks of the Triakidae (*Mustelus* spp.) and Carcharhinidae (*Rhizoprionodon* spp.), placentotrophic species with annual reproduction, have concomitant pre-ovulatory vitellogenic follicles and well-developed embryos [45,46]. Although the ultrasound examination did not allow the visualization of follicles for the female at the end of gestation, a new courtship and attempt at copulation by the males started a few days after birth and the female's return to the ocean tank. Considering that the average gestation period observed by Schaller [24] is 387 days and that the first neonate described here was born 12 days after the ultrasound examination, we assumed that the female was pregnant for approximately 1 year. Therefore, the description of organogenesis via this method characterized the final third of pregnancy. This study offers the first images of the development of a healthy embryo in the late gestation period of *T. obesus*.

Viviparous placental development has been recorded only in Selachii of the order Carcharhiniformes, in the families Leptochariidae, Triakidae, Hemigaleidae, and Carcharhinidae [22,47]. The literature describing the developmental stages for placental species is scarcer than that available for oviparous lecitotrophic species [42] or viviparous aplacental species with oophagia, histotrophy, and adelphophagy [48,49]. Although the stages for *T. obesus* have not yet been described, the whitetip reef shark embryo described herein, based on an ultrasound taken 12 days before birth, showed full structural development, with developed fins and skeleton, as well as internal organs without apparent anomalies. Behavior, feeding patterns, and routine veterinary examinations in the months following birth confirmed that both newborns were fully developed, healthy individuals.

#### *4.3. Birth and Feeding*

The ex situ reproduction of a shark species classified as vulnerable, such as *T. obesus* [19], highlights the need to integrate public aquariums in the consolidation of a viable ex situ population, according to criteria already outlined by Buckley et al. [50], and an essential step for this is the reproductive success and genetic variability of this population. Schaller [24] described four litters of *T. obesus* obtained between 2001 and 2003 in work conducted at the Steinhart Aquarium (USA), providing precise data on courtship, pregnancy, and parturition for the species. Portnoy et al. [51] described the birth of a stillborn whitetip reef shark at the BioPark Aquarium (USA), with no apparent reason for death or signs of trauma or malformation. It is important to note that Portnoy et al. [51] also pointed out that all parthenogenesis cases reported for sharks generated females, which initially suggested that the female pup born at AquaRio could have been a parthenogenetic event. Still, photographic records and recordings obtained in 2021 made it possible to follow courtship and copulation events among the *T. obesus* group, supporting the suspicion that the births were the product of sexual reproduction. The appearance of bite scars in June 2021 and the presence of a male pup also suggested that parthenogenesis did not take place. Furthermore, a DNA analysis is going to be conducted to establish if the same male was responsible for both shark pups. The diet of adult whitetip reef sharks consists mainly of teleosts (Holocentridae, Lutjanidae, Pomacentridae, Scaridae, Acanthuridae, Balistidae, Mullidae), which can comprise 90% of their food, although mollusks and crustaceans have also been reported [29,52]. The food preference displayed by the neonates reflected the described diet for adults, consisting mainly of teleosts, but with a relevant inclusion of penaeid shrimps.

#### **5. Conclusions**

The study compared the data available for the reproduction of *T. obesus* in the literature, based mostly on in situ observations, and the case reported by Schaller [24] for births registered at the Steinhart Aquarium (USA), making it possible to corroborate the courtship stages already reported, as well as to identify female escape strategies and contribute data on gestation duration and the initial development of neonates and their food preference. Furthermore, ultrasonography is a noninvasive examination technique and it allowed us to detail fetal organogenesis and study the viability of the fetus in whitetip reef sharks, with this being the first ultrasonographic description of fetal organogenesis in the species and the first report of reproduction of the species in a Brazilian aquarium.

Completing the reproductive cycle of endangered elasmobranchs in an ex situ environment is another relevant step that may lead to an effective plan for the preservation of these species. While natural populations are under pressure from habitat loss and overfishing, a healthy population kept safe from these impacts will enable a species reintroduction plan to be implemented after minimum living conditions are restored for the species in their original distribution. This strategy is even more vital when it comes to coral reefs, highly productive ecosystems with rich fauna and complex food webs. Existing literature reports several threats already affecting these communities, within which *T. obesus* is an essential predator. Reproductive success will offer hope for Brazilian species of Carcharhinidae, as many share characteristics with the whitetip reef shark, such as being viviparous placentotrophs. The lessons learned from whitetip reef sharks under human care may assist in the development of protocols that help with the proper handling of other carcharinid sharks in ex situ conditions, with the possibility of increasing the chances of successful reproduction of endangered species, such as the highly threatened sandbar shark *Carcharhinus plumbeus* (Nardo, 1827) and the small-tail shark *Carcharhinus porosus* (Ranzani, 1839).

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/ani12233291/s1, Video S1: Courtship and copulation of *Triaenodon obesus* (Rüppell, 1837), recorded at the Rio de Janeiro Marine Aquarium. Video S2: Failed attempt at courtship of *Triaenodon obesus* (Rüppell, 1837), recorded at the Rio de Janeiro Marine Aquarium.

**Author Contributions:** Conceptualization, S.R.S., V.T. and M.F.G.; methodology, V.T., S.P.B. and M.F.G.; validation, S.R.S., V.T., S.P.B. and M.F.G.; investigation, S.R.S., V.T., S.P.B., N.L.L.A. and M.F.G.; resources, R.F.V.; data curation, S.R.S. and V.T.; writing—original draft preparation, S.R.S., V.T., S.P.B., N.L.L.A., M.F.G. and R.F.V.; writing—review and editing, S.R.S.; visualization, S.R.S.; supervision, M.F.G. and R.F.V.; project administration, R.F.V.; funding acquisition, R.F.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Seasub 7 S.A. through the partnership with IMAM-AquaRio.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable, as this research did not involve humans.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors would like to thank the combined efforts of aquarists and veterinarians that help preserve the lives of thousands of fish, crustaceans, and mollusks of different species, among other groups, under a strict daily routine. Their ethological observations and constant surveillance helped to identify the first signs of courtship and pregnancy and build a detailed profile of the behavioral evolution of both pregnant females and neonates. We also thank both anonymous reviewers for their careful reading of our manuscript and their insightful remarks and suggestions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

### *Article* **Growth Characteristics of Long-Nosed Skate** *Dipturus oxyrinchus* **(Linnaeus, 1758) Inhabiting the Northeastern Mediterranean Sea**

**Nuri Ba¸susta 1,\* and Fatih Volkan Ozel <sup>2</sup>**

**Simple Summary:** Skates and rays generally have low fecundity, delayed maturation age and slow growth rates. These life history traits make them very vulnerable to commercial fisheries' activities, although they are not usually target species. Understanding age and growth is important for stock assessment of the long-nose skate. We believe that determining these can be the first step toward correct management actions. Little is known about the growth characteristics of the long-nose skate in the Northeastern Mediterranean Sea. Our study aimed at comparing three different growth models (the von Bertalanffy, the Robertson (logistic), and Gompertz) and also the absolute and relative growth characteristics of this species. Our results show that the relationship between the age of the long-nose skate and its total length is adequately explained by the Robertson (logistic) growth model, with the Gompertz growth model being the second best.

**Abstract:** This study aims to determine the age and growth characteristics of *Dipturus oxyrinchus* living in the Northeastern Mediterranean Sea and to present data that can provide a comparison with previous studies on the same subject. A total of 255 long-nose skates at a total length of 12.2–93.5 cm and weight of 8.34–3828 g were collected as non-target species from a commercial fishing boat. The male−female ratio was determined as 1:1.27. Using the von Bertalanffy equation and the Gompertz or logistic growth models, the growth parameters of *Dipturus oxyrinchus* were estimated as L∞ = 154.0, K = 0.064, t0 = −1.622; L∞ = 104.0, K = 0.35, *I* = 4.99; L∞ = 128.40, K = 0.19, *I* = 4.39 for all individuals, respectively. Maximum absolute growth was calculated as 9.33 cm at 5–6 years of age. Maximum relative growth at 1–2 years of age was estimated as 36.39%. Both absolute and relative growth were minimal in the 11–12 age group. The highest condition factor value was estimated as 0.416 in the 8-year-old group. As a result, the growth data of long-nose skates were obtained for the first time in the Northeastern Mediterranean Sea.

**Keywords:** age; growth; *Dipturus oxyrinchus*; long-nosed skate; Northeastern Mediterranean Sea

### **1. Introduction**

The long-nosed skate, *Dipturus oxyrinchus* (L. 1758), is a demersal species that lives at depths of 70−1230 m on sandy or muddy substrates and, rarely, on rocky and pebbly grounds; it is usually found at depths of 200−500 m [1–6]. The long-nosed skate ranges from Norway to Senegal, from the Northeast Atlantic to the Faroe Islands, Skagerrak (the strait connecting the Baltic Sea to the North Sea), the Canary Islands, the Madeira Islands and the Mediterranean Sea [3] and has little commercial value [4]. Sizes between 60 and 100 cm are common, but the largest recorded individual was 150 cm [1]. *Dipturus oxyrinchus* is globally recognized as a Near Threatened (NT) species by the International Union for Conservation of Nature (IUCN) [7]. Long-nosed skates have been studied satisfactorily by researchers during recent years in other parts of the Mediterranean Sea in terms of age, growth, distribution, systematics, length–weight relationships (LWR) and feeding

**Citation:** Ba¸susta, N.; Ozel, F.V. Growth Characteristics of Long-Nosed Skate *Dipturus oxyrinchus* (Linnaeus, 1758) Inhabiting the Northeastern Mediterranean Sea. *Animals* **2022**, *12*, 3443. https://doi.org/10.3390/ ani12233443

Academic Editor: Martina Francesca Marongiu

Received: 29 October 2022 Accepted: 2 December 2022 Published: 6 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

habits. However, no data on growth parameters are available for this species in the Northeastern Mediterranean Sea [8–12]. The presence of juveniles and mature males of the long-nosed skate was previously reported by Ba¸susta and Ba¸susta [13] from the same region. Griffiths et al. [14] concluded that Mediterranean long-nosed skates may have been genetically isolated from other stocks (e.g., Atlantic). This study aimed to determine the growth characteristics of long-nosed skates living in the Northeastern Mediterranean Sea and compare them with the data reported in other studies conducted in the same region.

#### **2. Materials and Methods**

#### *2.1. Collecting of Samples*

The long-nosed skate individuals were caught by a commercial trawler (F/V NIHAT BABA/31-A-1463) in Iskenderun Bay (Figure 1; 36◦29 200 N; 35◦05 973 E–36◦07 052 N; 35◦17 936 E–36◦07 148 N; 35◦17 978 E–36◦13 720 N; 35◦22 998 E–36◦13 650 N; 35◦23 032 N–36◦16 622 E; 35◦18 509 N) between May 2015 and June 2016. The samples were collected monthly. The bottom trawl gear used was equipped with a 42 mm stretched-mesh size net at the cod-end. Each hauling lasted three hours, with a trawling speed of 2.2–2.9 knots. Approximate sampling depths ranged between 100–150 m, 150–200 m and 200–400 m. After fishing, all samples were transported with ice to the laboratory of the Fisheries Faculty, Firat University. Total length (*L*) was measured as a straight line from the tip of the rostrum to the end of the tail to the nearest mm and body mass (*W*) was weighed witha1g accuracy for each individual.

**Figure 1.** Red dots indicate out of Iskenderun Bay, Turkey, where *Dipturus oxyrinchus* specimens were collected.

#### *2.2. Processing of Vertebral Centra*

A section of 10–12 vertebral centra was removed from the widest portion of the body of 255 long-nosed skate specimens (143 females and 112 males) and subsequently labeled, frozen and stored until further processing. Vertebrae were later thawed and cleaned of excess tissue, rinsed in tap water and then stored in a 70% ethanol solution. Three random vertebrae from each sample were removed from the ethanol and air-dried [15,16]. Smaller centra of less than 5 mm were fixed to a clear glass slide using resin (Crystol bond 509™) and were sanded with a Dremel™ tool to replicate a sagittal cut [17–19]. Vertebral sections (0.6 mm thick) were taken using a Ray Tech™ (Littleton, CO, USA) gam saw for large centra >5 mm in diameter [20]. Vertebral cross-sections were mounted on microscope slides using clear resin (Cytoseal 60; Fisher Scientific, Pittsburgh, PA, USA) [21].

#### *2.3. Age Assessment and Verification*

Vertebral cross-sections were examined under a Leica S8 APO™ (Singapore) microscope using LAS software (Version 4.8.0, Leica Microsystems Limited, Heerbrugg, Switzerland). One growth band was defined as an opaque and translucent band pair that traversed the intermedialia and clearly extended into the corpus calcareum (Figure 2) [22–24].

**Figure 2.** A vertebral cross-section of an estimated 12-year-old *Dipturus oxyrinchus* (total length = 93.5 cm, female) (BM, Birth Mark). White dots indicate opaque bands.

Each vertebral cross section was examined by two readers (reader 1 = NB and reader 2 = FVO). Reader 1 made two nonconsecutive band-counts of sampled vertebral cross-sections without prior information of the long-nosed skate's length or former counts. Reader 2 made two consecutive counts from 50 randomly selected vertebrae sections. Vertebral cross-sections that had an instability of more than two years between each reading were eliminated from further analyses. Count reproducibility was compared by the percent agreement (%PA) and coefficient of variation (%CV) [25], as well as the index of average percent error between the readers (%IAPE) [26]. All were determined using the following mathematical equations:

$$\text{PA} = \frac{\text{No.} \,\text{ag} \,\text{red}}{\text{No.} \,\text{read}} \times 100 \tag{1}$$

$$\text{CV}\_{\text{j}} = 100 \times \frac{\sqrt{\sum\_{i=1}^{R} \frac{(X\_{i\bar{j}} - X\_{\bar{j}})^2}{R - 1}}}{x\_{\bar{j}}} \tag{2}$$

$$\text{IAPE} = \frac{1}{R} \sum\_{i=1}^{R} \frac{|\mathbf{x}\_{ij} - \mathbf{x}\_{j}|}{\mathbf{x}\_{j}} \times 100 \tag{3}$$

where *R* is the number of readings; X*ij* is the count from the *j*th fish at the *i*th reading and *Xj* is the mean age calculated for the *j*th fish from *i* readings.

Pair-wise age-reader comparisons were independently generated by the two readers by making nonconsecutive band counts from a random sample of 50 vertebral sections [27]. All statistical tests were performed with R software version 4.1 [28], with a significance level set at 5%.

#### *2.4. Total Length–Weight Relationship and Growth Modeling*

The total length–weight relationship parameters of long-nosed skate were estimated according to the equation given below [29]:

$$\mathcal{W} = aL^b$$

where *L* = total length (cm); *W* = body mass (g); *a* is a constant of proportionality and *b* is the allometric factor. The deviance of the estimated *b* values for long-nosed skate from the hypothetical value of 3 (i.e., isometric growth) was tested by a *t*-test at the 0.01 significance level [29,30].

The absolute length growth (ALG) and the relative length growth (RLG) of long-nosed skate were calculated as:

$$\text{ALG} = (Lt + \Delta t - Lt) / \Delta t \tag{5}$$

$$\text{RLG } = (Lt + \Delta t - Lt/Lt) \times 100\tag{6}$$

where *Lt* is total length at the start of the time interval and *Lt*+Δ*<sup>t</sup>* is total length at the end of the time interval (Δ*t*) [30,31].

The values of condition factor were obtained with the formula:

$$Kn = \left(\frac{W}{Lb}\right) \times 100\tag{7}$$

where *W*, *L* and *b* are as defined above [32].

The observed length-at-age data of the long-nosed skate were used as the dependent variable and the age as the independent variable, with the three most-used models to describe the growth of fish. The first model employed was the von Bertalanffy [33] growth equation (VBGM). The VBGM was formulated by Beverton and Holt [34] as:

$$L\_t = L\_\infty \left(1 - e^{-K(t - t\_0)}\right),\tag{8}$$

where *Lt* is the expected total length at age *t* years; *L*<sup>∞</sup> is the asymptotic total length; *K* is the growth coefficient or curvature parameter indicating the rate at which long-nosed skates grow toward their *L*<sup>∞</sup> and *t*<sup>0</sup> is the theoretical age at zero total length.

The second growth model used was that of Gompertz [35], an S-shaped growth model (GGM) [36–38]:

$$L\_t = L\_\infty e^{-e^{-G(t-t\_i)}},\tag{9}$$

where *t*, *Lt* and *L*<sup>∞</sup> are the same as in the VBGM; *ti* is the age at the inflection point of the growth curve, i.e., the age at which the absolute growth rate starts to decrease and *G* is the instantaneous growth rate coefficient at age *t*i, where growth becomes asymmetrical.

The Robertson (logistic) model was also used as the last growth model. While *t*, *t*i*, Lt* and *L*<sup>∞</sup> are the same as in the previous models, *K* is a parameter that affects the rate of exponential growth:

$$Lt = \frac{L\_{\infty}}{(1 + e(-K(t - t\_i)))} \tag{10}$$

Akaike's information criterion (AIC) was used to compare the different growth models. A smaller value of the AIC indicates that the observed data are closer to the fitted model [29,39]. The AIC is defined as <sup>−</sup>2 times the maximum value of the log likelihood (*L*ˆ) plus 2 times the number of parameters (*p*) in the model including the estimated variance [40]:

$$\text{AIC} = -2\ln(\hat{L}) + 2p\tag{11}$$

#### **3. Results**

#### *3.1. Sample Composition, Sex and Vertebral Analysis*

A total of 255 long-nosed skates ranging from 12.2 to 93.5 cm in total length and 8.34–3828 g in weight were collected as bycatch or discard species from a commercial fishing vessel between May 2015 and June 2016 (Table 1). The sex ratio (M/F) was determined as 1:1.27. The ratio of females to males was not statistically different from the expected 1:1 ratio between sexes (*p* > 0.05). The size frequency based on total length and age group is presented in Figure 3.

**Table 1.** Descriptive statistics for *Dipturus oxyrinchus* inhabiting the Northeastern Mediterranean Sea.

SE: Standard Error; Min: Minimum; Max: Maximum.

**Figure 3.** Frequency distribution of *Dipturus oxyrinchus* inhabiting the Northeastern Mediterranean Sea, according to their total length and age.

#### *3.2. Age Estimation, Reading Precision and Age Bias*

In this study, estimated ages ranged from 0 to 12 for females and from 0 to 9 for males. Age estimation for males and females by two independent readers did not show a considerable variation (Figure 4a,b). The highest PA with lowest IAPE and CV were found for males than female (Table 2).

**Figure 4.** (**a**) Age bias graphs for *Dipturus oxyrinchus* inhabiting the Northeastern Mediterranean Sea. (**b**) Age bias plots for two readers for aging *Dipturus oxyrinchus* inhabiting the Northeastern Mediterranean Sea. Plots illustrate reference reader's age estimates on the *x*-axis; the mean difference (circles) and distribution of the differences between corresponding ages (vertical lines) are represented on the *y*-axis.

**Table 2.** Summary statistics for the coefficient of variation (CV), percentage of agreement (%PA) and index of average percentage error (IAPE) to determine the precision of age readings of *Dipturus oxyrinchus* inhabiting the Northeastern Mediterranean Sea.

*n* = sample size and R = number of readings.

#### *3.3. Length–Weight Relationships, Growth Patterns and Condition Factor*

#### 3.3.1. Length–Weight Relationships

The total length–weight relationships of *Dipturus oxyrinchus* for female, male and overall are presented in Figure 5. According to these results, positive allometric growth (*b* > 3) was demonstrated for the overall category. Regression results showed that the length and weight of the long-nosed skate was predicted significantly (r = 0.987, r2 = 0.974, F1,242 = 9092.525 *p* < 0.001) in all sexes. It is possible to state that 97% of the increase is due to that of the size of *Dipturus oxyrinchus* in all individuals for the present study. In addition, when *t*-test results related to the significance of the regression coefficients were analyzed (*t*-test = 95.355 *p* < 0.01), the proportion of the individual's length to weight was found to be important.

**Figure 5.** Total length−weight relationships of *Dipturus oxyrinchus* inhabiting the Northeastern Mediterranean Sea.

#### 3.3.2. Growth Characteristics

The growth parameters were estimated for overall and each sex separately from the VBGM, the Robertson (logistic) and Gompertz growth model using nonlinear regression analysis, as presented in Table 3 and Figure 6. The relationship between total length and age was adequately described by the Robertson (logistic) growth model followed by the Gompertz growth model. The von Bertalanffy model performed relatively weakly compared to other growth models based on AIC values.

**Table 3.** The growth parameters for *Dipturus oxyrinchus* inhabiting the Northeastern Mediterranean Sea derived from different growth models.

*L*<sup>∞</sup> = the asymptotic length (cm); *K* = the growth rate (year−1); *t*<sup>0</sup> = the time when *L* = 0 (year); *I* = the age at the inflection point; and AIC = Akaike Information Criterion.

#### 3.3.3. The Relative and Absolute Growth Rates and Condition Factor

The maximum absolute growth was estimated as 9.33 cm with 5–6 years of age. The maximum relative growth was calculated as 36.39% with 1–2 years of age. Both in absolute and relative growth were observed as minimum in the 11–12 age group (Table 4). Average condition factor value of population was estimated as 0.363. The highest and lowest condition factor values were estimated as 0.416 in age group 8 and 0.308 in age group 2, respectively (Table 5).

**Table 4.** The relative and absolute growth rates of for *Dipturus oxyrinchus* inhabiting the Northeastern Mediterranean Sea.

**Figure 6.** The different growth models fitted to the data of *Dipturus oxyrinchus* inhabiting the Northeastern Mediterranean Sea.

**Table 5.** The condition factor for *Dipturus oxyrinchus* inhabiting the Northeastern Mediterranean Sea.

#### **4. Discussion**

The overall length ranges recorded in our study were found to be smaller than those reported by Alkusairy and Saad [11] for the same species in Syrian waters (34.1–100 cm for females and 34.5–81.6 cm for males). Yigin and I¸smen [8] reported a total length of 14.9–100 cm in females in Saros Bay (the North Aegean Sea), while Bellodi et al. [12] reported these values as 15.2–86.5 cm in males and 10.4–117.5 cm in females in Sardinian waters. Finally, Kadri et al. [9] reported 16.5–105 cm in females and 15.5–95 cm in males in the Gulf of Gabès (Southern Tunisia, Central Mediterranean). These values are very close to the values we found in our study. The Robertson (logistic) growth model estimates indicated that the long-nosed skate showed sexual dimorphism with females larger than males. These observations are consistent with other studies for *Dipturus oxyrinchus* in the Mediterranean Sea (e.g., Yigin and I¸smen [8] in the North Aegean Sea, Kadri et al. [9] in Southern Tunisia and Bellodi et al. in Sardinian waters [12]). Sexual dimorphism appears to be a common feature for the Rajidae (e.g., blonde ray *Raja brachyura* [41] and thornback ray *Raja clavata* [42]). The percentage of females and males for all samples was 56.07% and 43.93%, respectively. This was not statistically different from the expected 1:1 ratio between the genders. All genders were equally distributed confirming the pattern proposed by Yigin and I¸smen [8], Kadri et al. [9] and Bellodi et al. [12]. Growth bands were highly legible and visible in cross-sections, with an easily recognizable birthmark. No staining technique was used to determine the age of *Dipturus oxyrinchus*, with a maximum age of 12 noted in females and 9 in males. Differences in age determination after age 9 for both sexes are related to the small number of individuals being sampled. The age estimation process ensured a high level of count repeatability among readers (IAPE = 0.53%; %CV = 0.75; %PA = 90.98) and no signs sensitive to bias were detected among readers. These precision values are acceptable [43]. The b parameter of the length–weight relationship of *Dipturus oxyrinchus* showed positive allometric growth for both genders in the current study. The estimated *b* values for *Dipturus oxyrinchus* by region are shown in Table 6, which are very close to our study's findings. Other b values were reported as 3.539 for the south coasts of Portugal by Borges et al. [44] and 3.40 for North Aegean Sea by Filiz and Bilge [45]. These values differ from those of our study, which may be due to the small sample size of the fish or the fact that samples were made in different seasons. Relative and absolute growth rates and condition factor for *Dipturus oxyrinchus* were calculated for the first time in our study and therefore no comparison with other studies could be made. Absolute growth rates indicate actual growth between two years (ages) in terms of weight or length. Absolute growth rate decreases with age (t) and provides information about the years (ages) when growth is highest. The way the absolute growth rate is calculated depends strongly

on the size the fish has reached. For comparison purposes, the relative growth rate may be more useful. This is used to determine age-related growth rate in natural populations [30].

**Table 6.** Length–weight relationship values for *Dipturus oxyrinchus* from different regions.

*n*: sample size, Lmin-max: minimum-maximum length, Wmin-max: minimum-maximum weight, a: intercept, b: slope of the equation, r2: coefficient of determination.

Bellodi et al. [12] emphasized that the Gompertz function provided the best fit among the four growth models examined. Additionally, Liu et al. [49] stated that multiple model applications should be tested in elasmobranch age and growth studies. They also indicated that the Robertson (logistic) and Gompertz models provide the best fit for small-sized demersal skates/rays living in deep water. In our study, according to the AIC values, the logistic and Gompertz models were found to be more appropriate in describing the growth parameters of *Dipturus oxyrinchus*. These results agree with Liu et al. [49] and Bellodi et al. [12]. The logistic parameters determined in our study show that females attain a slightly larger asymptotic TL∞ (103.54 cm) than males (103.23 cm). In addition, the K values of the long-nosed skate were found to be similar for both genders. These growth rates appear to be similar with other skate species of similar size in the Mediterranean Sea. Bellodi et al. [12] suggested that the best growth model was the Gompertz model and estimated L∞ as 127.55 cm for all genders. This result is very close to the asymptotic value (L∞ = 128.40 cm) calculated with the Gompertz model for all individuals in our study. Yigin and I¸smen [8] reported that, for both genders, the L∞ was 256.46 cm and the K was 0.04, the t0 value was −1.17 year and the maximum age was 9 years.  The largest individual captured in our study (93.5 cm) was smaller than that reported by the above researchers, despite being older. This shows that there may be a mistake in reading the older age rings. This leads to an overcalculation of the asymptotic length. For females, Kadri et al. [9] reported the L∞ as 123.9 cm, the K as 0.08, the t0 as −1.26 and the maximum age as 25; for males, the L∞ as 102.1 cm, the K as 0.12, the t0 as −1.18 and the maximum age as 22. These researchers found that the oldest individual was 105 cm in total length and 25 years old. Again, these values were very high compared to those in our study, which therefore indicated to us that mistakes might have been made in the reading of the age rings of the older individuals. These errors may have caused an underestimation of the L∞ value. Considering that the largest *Dipturus oxyrinchus* caught in nature is 150 cm in TL, the L∞ value should not be lower than this value.

#### **5. Conclusions**

Our work ensured basic growth parameters and the best fit among the three growth models for the long-nosed skate and that the Robertson (logistic) growth model was the best model to describe the species growth. Results of the research showed that the long-nosed skate has life history features similar to other Rajidae species in the Mediterranean Sea and are long-lived and slow growing. The present study has provided the first analysis of growth characteristics and data for *Dipturus oxyrinchus* for conservative management plans in the Northeastern Mediterranean Sea.

**Author Contributions:** Conceptualization, N.B.; methodology, N.B.; software, F.V.O.; validation, F.V.O.; formal analysis, N.B.; investigation, F.V.O.; resources, F.V.O.; data curation, F.V.O.; writing original draft preparation, N.B.; writing—review and editing, N.B.; visualization, N.B.; supervision, N.B.; project administration, N.B.; funding acquisition, N.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Scientific Research Projects Coordination Unit of Firat University, grant number: SUF.15.04.

**Institutional Review Board Statement:** The animal study protocol was approved by the Animal Experimentation Ethics Committee of Firat University, Turkey (protocol code 54/109 and 27.05.2015).

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data associated with this research are available and can be obtained by contacting the corresponding author.

**Acknowledgments:** Specimen collection was permitted through the General Directorate of Fisheries and Aquaculture of the Turkish Ministry of Agriculture and Forestry (permit no: 26.05.2015–2065- 42337). A part of this work was presented as Abstract at 11–13 May 2017, at the Ecology International Symposium in Kayseri, Turkey.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

### *Article* **Long-Term Changes in the Distribution and Abundance of Nine Deep-Water Skates (Arhynchobatidae: Rajiformes: Chondrichthyes) in the Northwestern Pacific**

**Alexei M. Orlov 1,2,3,4,5,\* and Igor V. Volvenko <sup>6</sup>**

**Simple Summary:** Changes in the spatial distribution and abundance of common skate species of the genus *Bathyraja* in the Russian waters of the Northwestern Pacific, where their fishery is considered promising, but still not developed, are traced. Over the past six decades, the boundaries of their ranges have slowly fluctuated near some average annual position and abundance trends have coincided with previously identified ecosystem rearrangements under the influence of climatic and oceanological changes. It will be useful to compare the results obtained with abundance trends in the near future, when fishing for these species is expected to intensify here, as well as with other areas where fishing pressure for these species already exists.

**Abstract:** Based on the analysis of long-term data from bottom trawl surveys (1977–2021), changes in the spatial distribution, position of the boundaries of the ranges and the catch rates of the nine most common deep-sea skates of the genus *Bathyraja* in the Russian waters of the Northwestern Pacific (*B. violacea*, *B. aleutica*, *B. matsubarai*, *B. maculata*, *B. bergi*, *B. taranetzi*, *B. minispinosa*, *B. interrupta,* and *B. isotrachys*) are considered. During the surveyed period, significant changes in the spatial distribution were observed, which are probably due to both subjective reasons (changes in the number of trawling stations, surveyed depths, etc.) and climatic changes. No monotonous displacement of the northern or southern boundaries of the range or its center in one direction was observed in any area of any species during the entire observation period. At the same time, shifts in the boundaries of the ranges of different species in different areas occurred for different decades, i.e., the boundaries of the ranges slowly fluctuated or "pulsed" near some average annual position. In general, from the 1970s to the 1980s, the number of skates grew; from the 1980s to the 1990s, it decreased; from the 1990s to the 2000s, it fluctuated at the achieved level; from the 2000s to the 2010s, it grew again; and from the 2010s to the 2020s, it decreased again. These trends coincide with previously identified ecosystem rearrangements under the influence of climatic and oceanological changes. The identification of links between changes in spatial distribution, range boundaries and catch rates with climatic and oceanological factors require separate additional studies.

**Keywords:** spatial distribution; catch dynamics; latitudinal range; western Bering Sea; Sea of Okhotsk; northwestern Sea of Japan; Sakhalin Island; Pacific waters; Kuril Islands; Kamchatka

**Citation:** Orlov, A.M.; Volvenko, I.V. Long-Term Changes in the Distribution and Abundance of Nine Deep-Water Skates (Arhynchobatidae: Rajiformes: Chondrichthyes) in the Northwestern Pacific. *Animals* **2022**, *12*, 3485. https://doi.org/10.3390/ani12243485

Academic Editor: Martina Francesca Marongiu

Received: 23 October 2022 Accepted: 6 December 2022 Published: 9 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Representatives of the order Rajiformes are widely distributed in the North Pacific Ocean. They play an important role in marine ecosystems [1], being mainly top predators consuming many valuable commercial invertebrates and fish [2–9].

Skates have biological features common to most cartilaginous fish (slow growth rates, late sexual maturation and low reproductive rates), which makes their stocks extremely vulnerable to fishing [10]. However, they are commercially important fishing targets in many countries [11], particularly in Southeast Asia, where skate wings are used for consumption [12]. The flesh of skates is rich in nutrients and contains an almost complete set of essential amino acids. Their livers, being rich in vitamin A, though less so than those of sharks, can be used for the extraction of oil for veterinary and medical purposes. Skates have strong and thick skin that may be used for the manufacture of many leather articles. The flesh of skates is suitable for processing to fish meat jelly, which is used in the preparation of various Japanese national dishes [12–15].

The biomass of skates within Russian waters of the northwestern Pacific amounted to 677 thousand tons [16], and in some areas they account for up to 10% of the total biomass of groundfish [17]. Their fishing in Russian waters has long been considered promising [18,19], but has not yet been developed. In the Far Eastern Fisheries basin of Russia in the 1990s, only 11.4–11.7 thousand tons were recommended for harvesting; and in the 2000s, 7.2–11.9; in the 2010s, 11.2–14.0; and in the 2020s, 11.2–11.3 thousand tons (first within estimated total allowable catch—TAC, then within possible catch—PC or recommended catch—RC). Skates were not separated by species, but their actual total catch never reached the recommended values, since these fish got on board Russian fishing vessels as bycatch only, were not retained and utilized and were subsequently discarded [20–22]. The situation has changed in recent years after some Russian companies started to export skate wings to China and the British Virgin Islands at prices between 0.9 and 1.3 USD per kg, while they also may be profitably marketed in Japan and Korea [23,24].

According to long-term monitoring data [25,26], at least 28 species of skates are regularly found within the Russian waters and adjacent areas of the North Pacific, and all of them have commercial value [27,28]. To date, there is published information about the general features of the spatial distribution and dynamics of catches of the most common and some rare species in the Pacific waters off the North Kuril Islands and southeastern Kamchatka [17,29–31]. Recently, data on the spatial distribution and dynamics of catches of Alaska skate *Bathyraja parmifera*, Okhotsk skate *Bathyraja violacea* and Aleutian skate *Bathyraja aleutica* throughout the North Pacific have also been described [32–34]. Additionally, there is information on the distribution and biology of the bottom skate *Bathyraja bergi* in the Russian waters of the Sea of Japan [35]. Nevertheless, the features of spatial distribution, dynamics of skate catches and their long-term changes have not been considered in detail for Russian waters of the northwestern Pacific. Meanwhile, insufficient information prevents the evaluation of correct values of the TAC, PC or RC differentiated by species and area, and the subsequent rational exploitation of skates stocks in the Russian waters of the northwestern Pacific, which during the period of increasing intensity of their fishing may negatively affect the condition of populations of these fish [36].

The goal of this paper is to identify changes in the spatial distribution and catch rates of the nine most common skates within the Russian Exclusive Economic Zone (EEZ) of the Northwestern Pacific (the northwestern Sea of Japan and the Sea of Okhotsk, western Bering Sea and Pacific waters off the Kuril Islands and Eastern Kamchatka) deep-water skates of the genus *Batyraja* based on the analysis of long-term data of bottom trawl surveys.

#### **2. Material and Methods**

Material was obtained from database [37] supplemented with data from the recent bottom trawl surveys conducted by the Pacific Branch of the Russian Federal Research Institute of Fisheries and Oceanography through 2021. We selected from it trawl stations located within the Russian EEZ (there were 39,420 of them, see Figure 1) and calculated the occurrence of all skate species over the research period (1977–2021). Species encountered less than 100 times are not included in the list. *B. parmifera* and *Bathyraja smirnovi* were also excluded from it, since until recently, information about the distribution of these two morphologically similar species was very contradictory [27,32,38–44]; as a result, research teams confused them with each other. Thus, nine of the most common species were identified, which became the objects of this study (Table 1).

**Figure 1.** Location of trawl stations in the surveyed area. The numbers indicate the geographical names mentioned in the text: 1 — Bering Sea, 2 — Sea of Okhotsk, 3 — Sea of Japan, 4 — Pacific Ocean, 5 — Navarin Cape, 6 — Koryak coast, 7 — Goven Cape, 8 — Kamchatka, 9 — Sakhalin, 10 — Kuril Islands, 11 — Paramushir Isl., 12 — Onekotan Isl., 13 — Simushir Isl., 14 — Urup Isl., 15 — Iturup Isl., 16 — Gulf of Anadyr, 17 — Olyutor Bay, 18 — Karagin Bay, 19 — Kronotsk Bay, 20 — Tatar Strait, 21 — Terpeniya Bay, 22 — Aniva Bay, 23 — Peter the Great Bay, 24 — Shirshov Underwater Ridge.

**Table 1.** Number of records of nine skate species within the Russian EEZ, 1977–2021.

Note: Species are sorted in descending order of total occurrence. Areas with records of particular species <15 (marked with \*) excluded from the analysis.

Identification of skates to species level during past surveys by research teams was conducted using available guides [45–47]. During recent surveys, scientists for this purpose additionally used publications [48,49]. Scientific (Latin) and common names of skates used in this paper are provided according to the most authoritative sources [50,51].

To analyze long-term changes in their distribution and abundance, we divided the time scale into six decades: (1) 1970s, (2) 1980s, (3) 1990s, (4) 2000s, (5) 2010s and (6) 2020s. The sample sizes thus obtained are shown in Table 2.

**Table 2.** Description of the data used for the analysis of spatial distribution and relative abundance of nine skate species in the Russian waters of the northwestern Pacific, 1977–2021.

Note: Surveyed area was calculated as the sum of areas covered by trawl hauls. Area of a trawl haul was calculated by multiplying trawl horizontal opening by hauling distance.

For each species in each decade and area, we found the maximum (max), minimum (min) and average (avg) latitudes of all catches to determine the northern and southern boundaries of the range, and its center in decimal degrees. The abundance of the species at each capture point, as well as the average for decades and areas, was estimated by catch per unit of effort (CPUE). For unit of effort, we took the area of a trawl haul (in square kilometers), which was calculated by multiplying the trawl horizontal opening by hauling distance, i.e., the unit of measure of CPUE is ind./km2.

In addition to the decadal changes in the four variables, which are described in detail in this paper, we estimated their total trends for all time using four simple regression models: (1) Linear: *y=a+b* · *x*, (2) Logarithmic: *y=a+b* · *ln(x)*, (3) Exponential: *y = exp(a + b* · *x)* and (4) Multiplicative: *y=axb*, where *<sup>x</sup>* is the decade number and *<sup>y</sup>* is the variable being approximated. In each particular case the best model was chosen by the minimum of the residual and *p*-value, and the maximum of the correlation coefficient (*r*) modulo. We hoped that the *r* sign would help classify trends as positive or negative. However, the variations in the variables turned out to be not monotonous, poorly approximated by smooth lines. Accordingly, the *r* values in most cases turned out to be non-significant at the 95% confidence level. All graphs, their equations, and *r* and *p* values are given in the Supplement, and the results are shown in the last columns of Tables 3 and 4.

**Table 3.** Trends in shifts of the northern (max) and southern (min) boundaries of a range and its center (avg) of nine skate species in the western Bering Sea (B), Sea of Okhotsk (O), northwestern Sea of Japan (J), and Pacific Ocean (P): "+"—northward shift, positive correlation (green); "–"— southward shift, negative correlation (brown); "0"— no shift (yellow); "?"— unknown, no records of species; "ns"— non-significant correlation (yellow).

#### **Table 3.** *Cont.*

Note: Correlation coefficients (*r*), *p*-values, regression equations and plots are presented in the Supplementary Figures S1–S27.

**Table 4.** Trends in CPUE changes of nine skate species in the western Bering Sea (B), Sea of Okhotsk (O), northwestern Sea of Japan (J), Pacific Ocean (P), and all areas combined (T): "+"—increasing, positive correlation (green); "–"—decreasing, negative correlation (brown); "0"—no changes (yellow); "ns"—non-significant correlation (yellow).

Note: Correlation coefficients (*r*), *p*-values, regression equations and plots are presented in the Supplementary Figures S28–S36.

#### **3. Results**

#### *3.1. Okhotsk Skate*

*Spatial distribution* (Figure 2). During the research period, the spatial distribution of the Okhotsk skate has changed significantly. In the first decade, it was caught only in the waters off the central and northern Kuril Islands and southern Kamchatka. The maximum catches were taken off the southwestern Kamchatka in the Sea of Okhotsk and off the central part of the East Kamchatka coast in the Pacific Ocean. In the 1980s, this species was found most widely within the surveyed area, i.e., from the Southern Kuril Islands in the south to Navarin Cape in the western Bering Sea in the north. In the Sea of Okhotsk, it was observed almost everywhere except in the southern part of the sea, and the maximum catches were recorded off both coasts of Kamchatka and the northern Kuril Islands. In the next decade, the distribution area of the Okhotsk skate was reduced. It was recorded only in the western Bering Sea, off the northern Kuril Islands, western and southeastern Kamchatka. Maximum catches were recorded in the Sea of Okhotsk off the central part of the West Kamchatka coast. In the 2000s, the distribution area of the Okhotsk skate expanded again. It began to occur off the eastern Sakhalin, in the northern Sea of Okhotsk, and off the southern and central Kuril Islands, but ceased to occur near the East Kamchatka coast. The maximum catches during this period were recorded off the southern and northern Kuril Islands and in the western Bering Sea. The spatial distribution of the Okhotsk skate in the 2010s did not fundamentally differ from that in the 2000s. However, the maximum catches shifted to the waters of the western and southeastern coasts of Kamchatka, and also remained off the northern Kuril Islands. At the beginning of the current decade, the Okhotsk skate disappeared from catches in the northern Sea of Okhotsk and off the eastern Sakhalin, but it was still observed in the western Bering Sea, off the northern Kuril Islands and western Kamchatka with maximum catches in the waters off the northern Kuril Island.

*Boundaries of the ranges* (Figure 3, Supplementary Figures S1–S3). In general, the distribution center of the Okhotsk skate shifted significantly to the north during the research period. At the same time, the trends in different areas differed significantly. Whereas the northern distribution boundaries in the western Bering Sea and Sea of Okhotsk have shifted northward, in the Pacific waters of the Kuril Islands and eastern Kamchatka, they shifted south. Similar trends in orientation, but not so pronounced, were noted in relation to distribution centers. The southern borders of the ranges showed slightly different dynamics. Whereas in the Sea of Okhotsk and Pacific waters the borders shifted northward, in the Bering Sea, they showed a weak but opposite trend.

*Catch dynamics* (Figure 4 and Supplementary Figure S28). In general, during the research period, the value of catches per unit effort of Okhotsk skate increased by almost 3 times. At the same time, an increase in catches was noted in the western Bering Sea and Pacific waters (especially sharp, starting from the fourth decade). In the Sea of Okhotsk, on the contrary, a decrease in the catches rates (especially at the beginning of the research period) was observed.

**Figure 2.** Decadal changes in the spatial distribution of the Okhotsk skate *Bathyraja violacea* in the Northwestern Pacific. Here and on Figures 2–9, red dots indicate stations where this species was not recorded.

**Figure 3.** Decadal changes in the position of the range of nine skate species in the Northwestern Pacific. The Bering Sea is in blue, Sea of Okhotsk in green, Sea of Japan in red, and Pacific Ocean in yellow, and all areas combined are in black. The northern boundaries of the ranges are shown by dash-dotted lines, southern boundaries by dashed lines, and middle of the ranges by solid lines with circles.

**Figure 4.** Decadal changes in the CPUE of nine skate species in the Northwestern Pacific. The Bering Sea is in blue, Sea of Okhotsk in green, Sea of Japan in red, and Pacific Ocean in yellow, and all areas combined in black. The solid lines with circles indicate the average CPUEs, and the dotted lines indicate the 95% confidence intervals.

**Figure 5.** Decadal changes in the spatial distribution of the Aleutian skate *Bathyraja aleutica* in the Northwestern Pacific. Designations as in Figure 1.

**Figure 6.** Decadal changes in the spatial distribution of the dusky-purple skate *Bathyraja matsubarai* in the Northwestern Pacific. Designations as in Figure 1.

**Figure 7.** Decadal changes in the spatial distribution of the whiteblotched skate *Bathyraja maculata* in the Northwestern Pacific. Designations as in Figure 1.

**Figure 8.** Decadal changes in the spatial distribution of the bottom skate *Bathyraja bergi* in the Northwestern Pacific. Designations as in Figure 1.

**Figure 9.** Decadal changes in the spatial distribution of the mud skate *Bathyraja taranetzi* in the Northwestern Pacific. Designations as in Figure 1.

#### *3.2. Aleutian Skate*

*Spatial distribution* (Figure 5). In the 1970s, the Aleutian skate was sporadically observed in catches off the southern and northern Kuril Islands, off southwestern Kamchatka and in the western Bering Sea. Its most frequent captures were recorded in the Gulf of Anadyr. The maximum catches were recorded in the central part of the Koryak coast and in the Pacific waters of the northern Kuril Islands. During the next decade, this species was most widely distributed within the study area. Its captures were recorded from the southern Kuril Islands in the south to the northern part of the Gulf of Anadyr in the north. In the Sea of Okhotsk, it was found almost everywhere except in the southwestern part of the sea. The maximum catches were in the Pacific waters of the northern Kuril Islands, off the central parts of the East Kamchatka and Koryak coasts and off the Navarin Cape. In the third decade, catches of the Aleutian skate were recorded mainly along the western and eastern coasts of Kamchatka, in the Pacific waters of the northern Kuril Islands and in the western Bering Sea. Several catches were recorded off the Eastern Sakhalin to the east of Terpeniya Bay. The maximum catches were recorded in the Olyutor Bay of the Bering Sea and in the Pacific waters of the northern Kuril Islands. In the 2000s, the distribution of the Aleutian skate changed significantly. Catches from Eastern Sakhalin have shifted far to the north. Off the western Kamchatka and northern Kuril Islands, their values have significantly decreased, and the captures themselves have shifted south to the Simushir Island. In the western Bering Sea, the distribution pattern has practically not changed, but the value of catches has increased significantly. In the next decade, the distribution of this species changed again. The occurrence has increased near northeastern Sakhalin, in the northern Sea of Okhotsk and off western Kamchatka. Individual catches have been recorded in the Pacific waters of the southern Kuril Islands. In the western Bering Sea, the catches in the Gulf of Anadyr have not only increased, but they have also been recorded almost within the entire Gulf of Anadyr. For several decades, the maximum catches were recorded in the Olyutor Bay, in the central part of the Koryak coast and off the Navarin Cape. At the beginning of the current decade, the Aleutian skate was found mainly in the western Bering Sea and also in a few catches in the Pacific waters of the northern Kuril Islands.

*Boundaries of the ranges* (Figure 3, Supplementary Figures S4–S6). In general, the research area was characterized by a shift of the distribution center in the north direction. At the same time, the northern border of the range has shifted significantly to the north in the Sea of Okhotsk. In the Bering Sea and Pacific waters, on the contrary, it has shifted to the south (in the latter area it is noticeably stronger). The center of the range has shifted to the north in the Sea of Okhotsk and Pacific waters and to a lesser extent to the south in the Bering Sea. The southern distribution boundaries in all three areas shifted to the north (most noticeably in the last two areas).

*Catch dynamics* (Figure 4, Supplementary Figure S29). During the research period, the value of the average catch per unit effort for all areas increased several times. At the same time, catches of the Aleutian skate in the western Bering Sea demonstrated the maximum growth. To a lesser extent, the increase in catches was typical for Pacific waters. In the Sea of Okhotsk, despite a sharp decrease in the number of catches per unit effort in recent years, the general trend showed a certain increase in them as a whole for the entire period of research.

#### *3.3. Dusky-Purple Skate*

*Spatial distribution* (Figure 6). In the first decade, there were no captures of the dusky-purple skate within the surveyed area. In the 1980s, it was most widely distributed, occurring from the southern Kuril Islands in the south to m. Navarin in the Bering Sea in the north. In the Sea of Okhotsk, it was found almost everywhere, with the exception of a deep-water basin in the southwestern part of the sea. The maximum catches in this period were typical for the waters of Eastern Sakhalin, central and northern Kuril Islands, eastern Kamchatka and the western Bering Sea. In the 1990s, the area of occurrence and number of catches of this species decreased significantly. The main catches during this period were observed in the waters of Western Kamchatka. The second most common area of occurrence was the western Bering Sea. Occasional catches were also recorded off Eastern Sakhalin and in the Pacific waters of the central Kuril Islands. In the next decade, the range of the dusky-purple skate decreased even more. It has practically disappeared from catches in the waters of Western Kamchatka and Eastern Sakhalin but has become more numerous in the western Bering Sea and Pacific waters of the Kuril Islands from Paramushir Island to Iturup Island. In the 2010s, the area of distribution of this species increased markedly. It reappeared in the waters of Eastern Sakhalin and Western Kamchatka. In the Pacific waters of the Kuril Islands, catches shifted northward; in the western Bering Sea, the occurrence expanded westward to the Karagin Bay. The maximum catches were recorded in the western Bering Sea. In the early years of the current decade, the dusky-purple skates were recorded in the western Bering Sea only and also caught once in the Pacific waters of southeastern Kamchatka. The maximum catches were recorded in the Olyutor Bay of the Bering Sea.

*Boundaries of the ranges* (Figure 3, Supplementary Figures S7–S9). In general, the period of research was characterized by a significant shift of the distribution center to the north. At the same time, the position of the northern border of the range in the Bering Sea has not changed, while in the Sea of Okhotsk and Pacific waters it has shifted to the south. The center of distribution in the Bering Sea also did not undergo displacement; while in the Sea of Okhotsk, it shifted somewhat southward; and in Pacific waters, on the contrary, to the north. The southern boundary of the range in all three areas shifted northward (to a lesser extent in the western Bering Sea).

*Catch dynamics* (Figure 4, Supplementary Figure S30). In general, for the research period, catches per unit effort increased many times over the course of six decades. At the same time, the maximum increase in catches was recorded in the Pacific waters and western Bering Sea. In the former area, after a sharp increase in the catch rates in the fifth decade, an equally sharp decline has been noted in recent years. In the Sea of Okhotsk during the research periods, there was an alternation of ups and downs of the catch rates, as a result of which the overall trend turned out to be weakly positive.

#### *3.4. Whiteblotched Skate*

*Spatial distribution* (Figure 7). In the 1970s, the whiteblotched skate in the research area was marked by a single capture in the waters of the southern Kuril Islands to the east of the Iturup Island. In the second decade, it was found in catches from southern Kuril Islands in the south to Navarin Cape in the Bering Sea in the north. In the Sea of Okhotsk, it was recorded near Eastern Sakhalin; in the northwestern part of the sea, off Western Kamchatka; and the northern Kuril Islands. The maximum catches were typical for the waters of the northern Kuril Islands and southeastern Kamchatka. In the 1990s, the range of the whiteblotched skate significantly decreased. It disappeared from catches in the waters of Eastern Sakhalin, southwestern Kamchatka and the Sea of Okhotsk waters off the northern Kuril Islands. The main area of occurrence and maximum catches remained the Pacific waters of the northern Kuril Islands. The second most important area was the western Bering Sea. The pattern of catch distribution of this species has not fundamentally changed in the next decade, except for individual catches off Eastern Sakhalin and in the Pacific waters of the southern and central Kuril Islands. The main area of occurrence and maximum catches were still the Pacific waters of the northern Kuril Islands. The value of catches in the waters of Western Kamchatka and the western Bering Sea has increased markedly in comparison with the previous decade. A similar distribution pattern of whiteblotched skate catches was typical for the 2010s. The only difference between these two adjacent decades was the appearance of this species in catches off Eastern Sakhalin and the northwestern Sea of Okhotsk. In addition to the Pacific waters of the northern Kuril Islands and the western Bering Sea, the maximum catches were also recorded in the waters of Western Kamchatka. In the current decade, the western Bering Sea has remained the main area of occurrence and maximum catches of whiteblotched skate. It was also recorded by single capture off the northwestern Kamchatka. Catches in the Pacific waters of the northern Kuril Islands have become occasional, although they have retained their high value.

*Boundaries of the ranges* (Figure 3, Supplementary Figures S10–S12). In general, in the research period, the center of the whiteblotched skate range showed a trend of displacement in the north direction. At the same time, the northern boundaries of its range in the Bering and Okhotsk Seas have hardly changed their position, and in the Pacific waters have noticeably shifted to the north (especially in the second decade). The center of the range in the Bering Sea had not undergone a significant shift, while in the other two areas it noticeably shifted northward. A similar pattern was observed with respect to the southern borders of the whiteblotched skate range.

*Catch dynamics* (Figure 4 and Supplementary Figure S31). In general, a noticeable increase in catches per unit effort (more than 2 times) can be noted for the study period. The catches of whiteblotched skates in Pacific waters demonstrate a particularly sharp increase, where after a considerable decrease in their value in the third and fourth decades, there was a sharp increase. In the Bering Sea during the entire period, catches showed a progressive increase. In the Sea of Okhotsk, despite the decrease in the value of catches in recent years, the general trends showed noticeable growth.

#### *3.5. Bottom Skate*

*Spatial distribution* (Figure 8). In the 1970s and 2020s, bottom skate was not recorded in catches within the study area. In the 1980s, the main area of its occurrence and maximum catches were the northwestern Sea of Japan. The second most important area was the waters around the southern Kuril Islands. Occasional catches were recorded in the southern part of the Tatar Strait (mainland) and Sea of Okhotsk to the southeast of the Aniva Bay. In the next decade, the northwestern Sea of Japan retained its importance as the main area of occurrence and maximum catches of bottom skate. The catches in the mainland part of the Tatar Strait began to be recorded somewhat northward in comparison with the previous decade. In the 2000s, the pattern of bottom skate distribution did not fundamentally change, but catches began to occur in the eastern part of the Tatar Strait off southwestern Sakhalin. In the next decade, the range of the bottom skate significantly expanded. It reappeared in the waters of the southern Kuril Islands, and also began to occur off northeastern Sakhalin. The maximum catches during this period were recorded in the waters of the Peter the Great Bay of the Sea of Japan.

*Boundaries of the ranges* (Figure 3, Supplementary Figures S13–S15). During decades two to five, the center of the range shifted slightly to the north on average. In the Sea of Japan, the northern border of the range has shifted noticeably to the north. The center of distribution and the southern border of the range, on the contrary, have shifted slightly southward. Data on other areas (the Sea of Okhotsk and Pacific waters) are not sufficient to judge the displacement of the boundaries of the range.

*Catch dynamics* (Figure 4, Supplementary Figure S32). The available data allow us to correctly judge only the dynamics of bottom skate catches in the Sea of Japan. During the period under review, from the second to the fifth decade, the catches per unit effort showed an increase with some decrease in the fourth decade. In general, the catches of this species were characterized by a trend of increasing catches both in the Sea of Japan and other areas.

#### *3.6. Mud Skate*

*Spatial distribution* (Figure 9). In the 1970s, the mud skate was not observed in the catches. During the second decade, this species was recorded in catches from the northern Kuril Islands in the south to the southern Gulf of Anadyr in the north. The maximum catches were recorded in the Pacific waters of the northern Kuril Islands and the bays of Eastern Kamchatka. In the 1990s, the mud skate was observed by occasional catches in the Pacific waters of the northern Kuril Islands and southeastern Kamchatka. In the next decade, the Pacific waters of the northern Kuril Islands again became the main area of occurrence and maximum catches. The second most important area was the western Bering Sea to the east of Navarin Cape. In the 2010s, the pattern of the distribution of catches in comparison with the previous period practically did not change. The Pacific waters of the northern Kuril Islands remained the main area of occurrence and maximum catches. However, in the western Bering Sea, large catches were recorded in the northern part of the Karagin Bay and Olyutor Bay. In recent years, the character of the spatial distribution has

hardly changed, but in the western Bering Sea, the main catches have shifted to the area of Navarin Cape.

*Boundaries of the ranges* (Figure 3, Supplementary Figures S16–S18). In general, the center of the range of the mud skate showed a southward shift from the second to the third decade, after which it shifted to the north. At the same time, the northern borders and distribution centers in the Bering Sea and Pacific waters showed trends of displacement to the south. Additionally, the trend of displacement of the southern boundary of the range in Pacific waters had a southward direction, while in the Bering Sea it was directed in the opposite direction.

*Catch dynamics* (Figure 4, Supplementary Figure S33). Catches per unit effort of mud skate both in the Bering Sea and Pacific waters showed similar dynamics, i.e., increasing in the second, fourth and fifth decades and a fall in the third and sixth decades. Nevertheless, trends showed an increase in the value of catches per unit effort during the entire research period, both for each area and for all areas combined.

#### *3.7. Whitebrow Skate*

*Spatial distribution* (Figure 10). In the first decade, catches of the whitebrow skate in the research area were not registered. In the 1980s, the area of maximum occurrence and catches was the western Bering Sea from the central part of the Koryak coast to Navarin Cape. The second most important area was the Kronotsk Bay. Single captures of this species have been recorded in the Olyutor Bay, off northwestern and southwestern Kamchatka, the central part of the Sea of Okhotsk and the Pacific waters of the northern Kuril Islands. In the 1990s, the whitebrow skate was recorded only in a single catch in the Pacific waters of the northern Kuril Islands to the east of Onekotan Island. In the 2020s, the pattern of its distribution was close to that of the 1980s, but the maximum occurrence and catch rates were characteristic of the Olyutor and Karagin Bays of the Bering Sea. In the 2010s, the range of the whitebrow skate has significantly expanded, and the number of catches has increased markedly. In addition to the western Bering Sea, which remained the main area of occurrence and maximum catches, this species began to be frequently observed along the entire Western Kamchatka and in the Pacific waters of the northern Kuril Islands. Also, its single capture was registered to the northeast of Sakhalin. In recent years, the whitebrow skate had been recorded only in the western Bering Sea from the northern part of the Karagin Bay to Navarin Cape.

*Boundaries of the ranges* (Figure 3, Supplementary Figures S19–S21). From the second to the third decade, as a whole, a significant southward shift of the center of the range of the whitebrow skate was observed, after which it consistently and significantly shifted to the north. The northern boundaries of the distribution in all areas showed a tendency to shift in a southerly direction (well expressed in Pacific waters and less noticeable in other areas). The center of the range shifted slightly to the south in the Bering Sea and more noticeably in Pacific waters, while in the Sea of Okhotsk it showed the opposite tendency. Similar trends were noted in relation to the southern borders of the range of the whitebrow skate.

*Catch dynamics* (Figure 4, Supplementary Figure S34). The dynamics of catches per unit effort of the whitebrow skate as a whole for all areas of its occurrence showed a series of their growth (2nd and 4th decades) and fall (3rd and 5th decades). Nevertheless, all areas combined were characterized by positive trends in the increase in catches per unit effort, i.e., the strongest in the Pacific waters and less pronounced in the Bering and Okhotsk Seas.

**Figure 10.** Decadal changes in the spatial distribution of the whitebrow skate *Bathyraja minispinosa* in the Northwestern Pacific. Designations as in Figure 1.

#### *3.8. Sandpaper Skate*

*Spatial distribution* (Figure 11). During the entire period of research, the sandpaper skate was observed in the western Bering Sea only. In the fourth and sixth decades, catches of this species in the research area were not recorded. In the 1970s, it was found in catches in the western Bering Sea from Goven Cape to Navarin Cape, and the maximum catches were recorded in the central part of the Koryak coast. In the 1980s, its number and distribution area in the western Bering Sea increased markedly. Its catches began to be observed south of Olyutor Bay in the waters of the Shirshov underwater ridge, as well as in the southern part of the Gulf of Anadyr. In the 1990s, the distribution pattern of the sandpaper skate in the western Bering Sea did not undergo significant changes, except for the extension of the range northward to the central part of the Gulf of Anadyr and absence of catches in the Olyutor Bay and waters of the western part of the Koryak coast. In the 2010s, the sandpaper skate was recorded in only a few catches off the Navarin Cape.

**Figure 11.** Decadal changes in the spatial distribution of the sandpaper skate *Bathyraja interrupta* in the Northwestern Pacific. Designations as in Figure 1.

*Boundaries of the ranges* (Figure 3, Supplementary Figures S22–S24). The fragmentary nature of the data and the short time series of observations do not allow us to judge with certainty about the dynamics of the position of sandpaper skate range boundaries. In general, we can only note a trend towards a southward shift of the northern border and a trend towards a northward shift of the center of the range and its southern border.

*Catch dynamics* (Figure 4 and Supplementary Figure S35). It is difficult to judge the dynamics of sandpaper skate catches in the western Bering Sea due to the lack of a continuous series of data. Nevertheless, the trend of catch per unit effort during the study period demonstrated their sharp decline.

#### *3.9. Challenger's Skate*

*Spatial distribution* (Figure 12). In the first decade, there were no captures of the Challenger's skate within the surveyed area. Occasional catches in the waters of southeastern Kamchatka and the Pacific waters of the northern Kuril Islands were recorded, respectively, in the 1990s and the current decade. This species was most widely distributed in the study area in the 1980s. During this period, it was most numerous in the Sea of Okhotsk, occurring mainly in the central part from southeastern and northeastern Sakhalin to southwestern Kamchatka. In Pacific waters, this species has been recorded in catches from Urup Island to Avacha Bay in Kamchatka. The maximum catches were recorded in the central part of the Eastern Sakhalin, southwestern Kamchatka and in the central part of the Sea of Okhotsk. The 2000s were characterized by the capture of the Challenger's skate exclusively in the Pacific waters of the Kuril Islands from Urup Island to Paramushir Island with maximum catches eastward to Onekotan Island. In the 2010s, catches of this species in the Pacific waters of the Kuril Islands were occasional, and its greatest occurrence was observed off Eastern Sakhalin and Western Kamchatka. The maximum catches were recorded in the northern Sea of Okhotsk off the coast of northwestern Kamchatka.

**Figure 12.** Decadal changes in the spatial distribution of the Challenger's skate *Bathyraja isotrachys* in the Northwestern Pacific. Designations as in Figure 1.

*Boundaries of the ranges* (Figure 3, Supplementary Figures S25–S27). In general, the center of the range of the Challenger's skate was characterized by a series of shifts to the north (2nd and 4th decades) and south (3rd and 5th decades). At the same time, the most complete data are available only for Pacific waters. The northern border of the range in this area showed a well-marked trend of southward displacement; the center of the range also shifted to the south (but less sharply), and the southern border of the range, on the contrary, shifted to the north.

*Catch dynamics* (Figure 4, Supplementary Figure S36). Generalized data for all areas showed that Challenger's skate catches per unit effort increased during the study period, with the exception of recent years. Nevertheless, in general, they showed a positive trend. In the Pacific waters and Sea of Okhotsk (despite the limited data in the latter area), catches showed very similar long-term dynamics in increasing catches during the entire study period.

#### **4. Discussion**

#### *4.1. Spatial Distribution*

An analysis of long-term changes in the spatial distribution of nine studied species showed that in the 1970s only three species were recorded in catches: Okhotsk, Aleutian and sandpaper skates. We see the reason for this situation as follows. Two species (whiteblotched and whitebrow) were first described only in the late 1970s [52] while bottom and mud skates were described several years later [53]. However, two species described much earlier, Challenger's skate at the end of the 19th century [54] and dusky-purple skate in the middle of the 20th century [55], were also not recorded in catches. The reason for this, in our opinion, was the lack of field guides to skates, which were developed and published only in the early 1980s [46]. This did not allow species identification during surveys in the 1970s, and therefore all rays caught were identified only to the genus level. Additionally, data for the first decade include survey results for only one year for the Sea of Japan (1978) and three years (1977–1979) for the remaining areas. Therefore, they should hardly be considered as representative ones.

The vast majority of species differed in the widest distribution within the surveyed water area in the 1980s. The most likely reason for this is the specifics of surveys. During this period, research funding was at its maximum. The research fleet consisted of dozens of vessels, including medium- and large-tonnage ones, capable of bottom trawling to a depth of about 2 km [37,56,57]. This period was characterized by the maximum number of stations, surveyed depths and areas. All this, in our opinion, led to the most complete account of skates within the study area.

In the 1990s, the surveyed area, number of trawl stations (more than 2 times) and range of surveyed depth significantly decreased. This was due to the collapse of the former USSR and economic crisis that followed when research funding was significantly reduced. In the 2000s, with the improvement of the economic situation in Russia, the amount of funding for scientific research increased slightly, which made it possible to significantly increase the number of stations and expand surveyed areas. However, the range of surveyed depths noticeably narrowed, which affected the occurrence of skates in catches. In the 2010s, the number of hauls compared with the previous period hardly changed. However, the range of surveyed depths became even narrower, which resulted in a large degree of underestimation of skates in the study area. Additionally, it should be noted that in the fourth and fifth decades, the total surveys were not carried out in the waters of Eastern Kamchatka, or in the fifth decade in the Pacific waters of the central Kuril Islands. The data for the last (sixth) decade cannot be considered representative since they were obtained for only two years (2020–2021).

Thus, the nature of the data presented by us on the spatial distribution of nine species of deep-sea skates within the Russian EEZ in the northwestern Pacific is largely due to some artifacts. However, there are no other such detailed and long-term observations. The value of the presented data lies in the fact that, despite their scarcity, they still make it possible to present the nature of the spatial distribution and its long-term changes in nine species of skates in the study area, to determine the areas of their main concentrations based on the catches and assess the fishing potential.

Until recently, the features of the spatial distribution of the Okhotsk skate were described in general terms only for the Pacific waters of the northern Kuril Islands and southeastern Kamchatka [17,29]. The distribution of this species throughout the North Pacific, including the waters of Russia, is reviewed in a recent publication [33]. Despite the fact that the same value (ind./km2) was used to characterize the distribution of catches in our and the mentioned work, the data on the long-term dynamics of the distribution turned out to be difficult to compare, since different time intervals were used and the visual representation of the features of the distribution of species varies greatly.

Data on the spatial distribution of the Aleutian, whiteblotched, whitebrow, mud, dusky-purple and Challenger's skates have so far been limited only to information on the occurrence and general features of the localization of catches and their magnitude in the Pacific waters of the northern Kuril Islands and southeastern Kamchatka, obtained in the 1990s [17,29–31]. Data on the distribution of the sandpaper skate in the northwestern Pacific have not yet been presented in the literature. Therefore, this article significantly expands the understanding of the occurrence, size of the catches, features of the spatial distribution and its long-term changes of the nine most common species of deep-sea skates within the Russian EEZ

#### *4.2. Boundaries of the Ranges*