Philippa Kaur

Massive drift of pumice along the northeastern coast of Okinawa IslandA large amount of pumice stones reached and was deposited along the northeastern coast of Okinawa Island, that were brought by strong seasonal northeasterly winds (Supplementary Video 1). The pumice was thought to be brought by the Kuroshio countercurrent from sites near the Ogasawara Archipelago 1300 km away. Because the Kuroshio countercurrent is composed of various medium-sized eddies in the ocean, the current does not always flow in one direction and as a continuous flow27,28. The pumice drift was more strongly controlled by the seasonal northwesterly winds to be transported to Okinawa across the Philippine Sea (Fig. 1a). The pumice raft reached the northern part of Okinawa approximately 2 months after the eruption (Figs. 2, 3 and 4). According to a very recent report, the pumice clasts were drifting ashore in Thailand (traveling 4000 km-long distance) across the South China Sea within half a year of this eruption29. Most pumice stones were gray, but some pumice was banded, and others were black reflecting some compositional variation25,29 (Fig. 2d,e). The Kuroshio Current is faster than the Kuroshio countercurrent27, so some pumice clasts have already reached the main island of Japan25. Tracking the dispersal of the pumice will allow a better forecasting model based on observed raft trajectories by considering exact wind effects in the Philippine Sea30.Figure 2An example of a natural beach on Okinawa Island where pumice has washed ashore. (a) Appearance of natural sandy beaches on the northern part of Okinawa Island (Ibu beach, Kunigami Village, 26°75′57.88″ N, 128°32′23.32″ E). Photo was taken on 24 October 2021. Pumice drifted onto the sandy beach and formed a striped pattern. The white-capped waves indicate on the place where the reef edge exist. The white arrow points to the mangrove river estuary corresponding to Fig. 9. (b) Estimation of the pumice sedimentation depth on the original sand beach surface. (c) The high tide zone of the natural sandy beach is covered with pumice pebbles and stones. Yellow arrows indicate black pumice stones. Scale bar: 10 cm. (d, e) Front and back of examples of relatively large pumice stones from the same beach. The left image is mostly light brown, whereas the right image is almost black. Scale bars: 5 cm.Full size imageFigure 3Short-term migration of pumice from beaches as revealed by stationary observations. These four photos were taken at two sites on northern Okinawa Island on two consecutive days, 23 and 24 October 2021. (a, b) A sandy beach along the Sate Coast (26°78′84.56″ N, 128°22′30.57″ E). It was windy on the first day, and pumice stones were washed up with the waves. Almost all the pumice stones were removed from the beach and transported offshore on the following day. The black arrow in photo (a) indicates Cape Hedo, the northernmost tip of Okinawa Island. (c, d) At this gravelly beach (26°80′83.25″ N, 128°23′38.56″ E), pumice fully covers the seawall on the first day, but all of the pumice stones washed away, leaving the original gravels, on the following day. The white arrow in each photo indicates an identical marker stone placed on the beach. Weather data of northern Okinawa (https://www.data.jma.go.jp/obd/stats/etrn/view/daily_a1.php?prec_no=91&block_no=0901&year=2021&month=10&day=23&view=g_wsp) and tidal data (Naha: 26°13′ N, 127°40′ E) (https://www.data.jma.go.jp/gmd/kaiyou/db/tide/genbo/genbo.php) are provided by Japan Meteorological Agency.Full size imageFigure 4Pumice stones settled by marine organisms. (a) Pumice collected from Ibu beach on 31 October 2021. Two marine benthos coexist close together on a pumice stone. Scale bar: 1 cm. (b) Enlarged image of the Lepas barnacle. Scale bar: 3 mm. (c) Enlarged image of the bryozoan. Scale bar: 3 mm. (d) Stereo microscopic image of pumice pebbles of a few millimeters in diameter collected from Ibu beach on 15 January 2022. The light brown coloration indicates some algal/cyanobacterial growth on the pumice. Scale bar 1 mm. (e) Red autofluorescence was detected from pumice pebbles. Image corresponds to (d). Autofluorescence from microalgae was confirmed by Supplementary Fig. 2. Scale bar 1 mm. (f) Enlarged image of the center of the figure of (e) shows red autofluorescent signals with a diameter of 10–30 µm. Scale bar: 200 µm.Full size imageChanges in the coastal landscape: natural beaches and estuariesMarine calcifiers, including corals, calcareous algae, and foraminifers, produce white sandy beaches on Okinawa Island. However, the gray pumice drifting ashore changed the white sand beach, especially along the northeastern coastline. We observed several lines of pumice aggregations, suggesting that pumice was brought ashore by wavefronts several times produced by a strong north wind at the tide lines (Supplementary Video 1; Fig. 2a). At the same sampling site, the thickest depth of beached pumice was more than 30 cm (Fig. 2b; Supplementary Video 2). Most of the pumice stones were from 0.5 cm to 3 cm in diameter, with a few black pumice stones included (Fig. 2c: yellow arrow). Pumice stones arrived at the estuaries of some brackish rivers (Fig. 8, Supplementary Fig. 1a) and mangrove forests in northwest Okinawa (Fig. 9).Pumice stones and pumice rafts show dynamic behavior in a short period. We captured photographs 24 h apart at two positions on the shore of Okinawa, which allowed us to compare the pumice dynamics during this period (Fig. 3). Within that time frame, there were two high tides, and the tide level changed by up to 170 cm. As seen in Fig. 3a, on the first day, the coast was covered with pumice, and floating pumice could be seen on the seafront. The north wind was strong that day, as shown by the relatively high waves near the shore as well as white‐crested waves near the reef edge. By the following day, most of the pumice had been moved offshore by tides and winds (Fig. 3b), indicating that newly beached pumice raft deposits were removed quickly from open beach areas. At another site on a gravelly beach, pumice fully covered the seawall on the first day, but almost all of the pumice stones were washed away, leaving the original gravels, on the following day (Fig. 3c,d). Japan Meteorological Agency (Oku station: 232 m above sea level, latitude 26°50.1, longitude 128°16.3′) reported that northerly winds were blowing (mean wind speed: 3.4 m/s) on 23rd October in northern Okinawa. The following day, the wind direction changed to the east-southeast; blowing offshore (mean wind speed: 2.9 m/s), resulting in the dramatic removal of pumice form the coast (Fig. 3). These observations indicate that surface winds rather than ocean currents had a strong influence on the raft trajectory and residence time on beaches, and are consistent with past research5. These observations lead us to expect that the pumice rafts will disappear from the coast of Okinawa fairly quickly, but in fact, there have been many cases where they have come back again in a few days. Although the overall amount of pumice drifting has been decreasing, a small amount of pumice has been drifting in coastal area of Okinawa in May, 202231. It is unlikely that large amounts of pumice will drift repeatedly throughout Okinawa Prefecture as reported in this report, but it should be noted that detached pumice material remains in beach and river runoff.Biofouling of sessile organisms on pumice arriving to OkinawaIt is noteworthy that the pumice rafts traveled over the deep Philippine Sea for over 2 months, and on arrival in Okinawa there was little to no biofouling of the pumice (Fig. 2). Some stranded pumices observed on Okinawa beaches had become habitats for sessile organisms (Fig. 4), as reported in previous studies1,2,3,4,5,6,29. Goose barnacles (Lepas sp.) without external damage to the shell were the most abundant species observed on the pumice (Fig. 4b). Lepas is a common biofouling taxon distributed globally and plays a role in biofouling as a foundation organism. The shell growth rate is more than 1 mm/day in some Lepas species32 suggesting that the Lepas had been growing on the pumice for about two weeks. Measurements of the shell size of Lepas attached to the pumice collections conducted in the same area (Supplementary Video 2) showed a bias toward larger sizes in the second collection (5.92 ± 3.86 mm (average ± S.D.), n = 75, 13 November 2021) than in the first one (3.43 ± 1.08 mm, n = 21, 31 October 2021), and significant differences were detected between the measurement periods (Mann–Whitney U test, p  More

Temporal differentiation of wild populations of T. hispidus and climatic parametersComparison of morphometric features of T. hispidus shells collected in different years in two geographic regions, i.e., Wrocław and Lubawka, showed significant differences depending on the year of collection. The largest number of differences was revealed in shells from Wrocław (Figs. 1 and 2A; Additional file 2: Table S1). Out of 210 comparisons (15 pairs of collection years × 14 features), 84 were statistically significant (Additional file 2: Table S2), e.g., shell diameter (D) was significantly different in 11 cases, shell height (H) and shell width (W) in 10 cases, body whorl height (bwH), the number of whorls (whl), umbilicus major (U) and minor (u) diameters in 9 cases and aperture height/width ratio (h/w) in 7 cases. Nine features obtained more than 10% difference between shells in at least one comparison of mean values, e.g., U 24%, u 19%, H 16% and D 15% (Additional file 2: Table S2). Umbilicus major (U) and minor (u) diameters showed the largest average percentage difference, i.e., 12% and 10%, respectively, in comparisons of all years.Figure 1Shells of Trochulus hispidus collected in different years in Wrocław.Full size imageFigure 2Changes in: mean values of selected morphometric features of shells collected in various years in Wrocław (A) as well as the mean temperature (B) and the relative humidity (C) recorded in four seasons in Wrocław in eight-year period. Abbreviations: D—shell diameter (in mm), H—shell height (in mm), h/w—aperture height/width ratio, whl—number of whorls. The summary statistics for A is included in Table S1 and original data in Table S10 in Additional file 2.Full size imageFor snails from Lubawka, out of 84 comparisons (6 pairs of collection years × 14 features) only 8 were statistically significant (Additional file 2: Table S3). The shells differed significantly in their aperture height (h) and width (w) in 3 comparisons. The h feature showed the percentage difference up to 9% (Additional file 2: Table S3) and the largest average difference was 4.5%.Besides the phenotypic variation, climatic parameters also showed high fluctuations in the studied period (Fig. 2B,C, Additional file 2: Table S4). The maximum difference reported between temperature parameters in some years prior to sample collection in Wrocław was up to 3.7 °C for the maximum winter temperature, while the maximum difference in the relative humidity was up to 11% for autumn. The maximum temperature difference in Jelenia Góra close to Lubawka was up to 3.5 °C for the minimum winter temperature, while the relative humidity differed at most by up to 8% in summer.Differences in shell morphometry under various climatic conditionsThe distinction between shells collected in individual years and changes in climatic parameters along the same period suggest that these differences can be associated with the climate. Therefore, we calculated the average value of a given climatic parameter for each season and studied region and next divided the collected shell data into two groups according to this value. The first group included the shells that developed in conditions above this average and the second below this average (Additional file 2: Table S5). The differences between these groups were statistically significant for 15 out of 16 considered climatic parameters for at least two shell features (Fig. 3). Similarly, each of 14 features significantly separated the groups based on at least two climatic conditions. The results demonstrated that the mean winter temperature substantially influenced nine morphometric shell features, whereas eight characters were changed due to the maximum winter temperature as well as the mean and minimum temperatures in spring, summer and autumn. Umbilicus major (U) and minor (u) diameters as well as umbilicus relative diameter (U/D) were significantly different in 14 pairs of groups characterized by various climatic parameters. In 11 pairs, the height/width ratio (H/W) was significantly different and shell height (H) in 10 pairs.Figure 3Mean percentage differences in morphometric features between shells that were grown in different conditions. The shells were divided into two groups according to the average value of a given climatic parameter for each season and studied region. The first group included the shells that developed in conditions above this average and the second below this average. Positive values indicate that the given feature was greater in the first group, whereas negative values indicate that this feature was greater in the second group. Dendrograms cluster the features and the parameters according to their similarity in the percentage differences. Values marked in bold indicate statistically significant differences between the compared groups of shells. Values at the dendrogram nodes indicate significance assessed according to approximately unbiased test (au) and bootstrap resampling (bp).Full size imageThe umbilicus diameters (u and U) as well as umbilicus relative diameter (U/D) clustered together in the dendrogram based on the mean percentage difference, which indicates that they similarly responded to climatic conditions (Fig. 3). The features u and U revealed the strongest average increase of all features, from 4.1 to 10.5% in shells developed in higher temperatures in all seasons. The largest percentage difference exceeding 10% was recorded for groups separated according to the mean summer and autumn temperatures as well as the maximum summer and minimum autumn temperatures. The U/D ratio was also significantly greater with the mean percentage difference of 2.8–7.6% in shells grown under high temperatures for all seasons and almost all temperature types. On the other hand, the u and U diameters as well as the U/D ratio were on average by 3.7–6.0% significantly smaller in shells developed under higher humidity in summer and winter.The height/width shell ratio (H/W) was grouped with H and bwH features in the dendrogram and was on average by up to 3.6% significantly smaller in shells grown under higher temperatures in all seasons for almost all types of parameters. The maximum winter temperature caused a significant increase, on average by ca. 3%, in shell height (H) and body whorl height (bwH), whereas higher temperatures in other seasons led to their decrease by up to 3.4% (Fig. 3).The shells that were grown in autumn with a relatively high maximum temperature were characterized by ca. 3% significantly smaller aperture height (h) and aperture height/width ratio (h/w), which were clustered together in the dendrogram (Fig. 3).Other four features, shell diameter (D), number of whorls (whl) as well as shell (W) and aperture width (w), formed an additional cluster in the dendrogram (Fig. 3). All of them were on average significantly greater in shells collected one year after winter that was characterized by relatively higher mean and maximum temperatures. The percentage difference was greater, with 3.6–3.9% for W and D.In the dendrogram, the climatic parameters were clustered in several groups indicating their similar influence on the morphometric features of shells (Fig. 3). There are separate clusters for temperature and humidity parameters with the exception of the autumn maximum temperature and autumn humidity, which are grouped together. The other temperature parameters for warmer seasons are separated from those for winter, which indicates that they differently influenced the shell morphometry.Correlations between morphometric shell features and climatic parametersThe influence of climatic conditions on the shells collected in individual years was also assessed using Spearman’s correlation coefficient between the morphometric features and climatic parameters (Fig. 4). Of 224 potential relationships 113 were statistically significant. The spring mean temperature was significantly correlated with 10 morphometric features. Summer humidity and six temperature parameters, i.e., the minimum temperatures as well as the spring and winter maximum temperatures, significantly correlated with eight shell features. Minor umbilicus diameter (u) and umbilicus relative diameter (U/D) were significantly correlated with almost all climatic parameters, i.e., 15, umbilicus major diameter (U) and height/width ratio (H/W) with 13 and the ratio of umbilicus minor to its major diameter (u/U) with 11.Figure 4Spearman’s correlation coefficients between morphometric features of shells with climatic parameters under which the snails were grown. Dendrograms cluster the features and the parameters according to their similarity in the coefficients. Values marked in bold are statistically significant. Values at the dendrogram nodes indicate significance assessed according to approximately unbiased test (au) and bootstrap resampling (bp).Full size imageAs in the case of percentage difference, we can also recognize groups of morphometric features that were similarly correlated with climatic parameters (Fig. 4). Features U/D, U and u were significantly positively correlated with all or almost all temperature parameters for four seasons with the coefficients up to 0.34, 0.30 and 0.36, respectively. On the other hand, the significant correlation coefficients between these features and the humidity in spring, summer and winter were negative and reached − 0.34.Another group of features included shell height/width ratio (H/W), shell height (H) and body whorl height (bwH) (Fig. 4). All of them showed significant negative correlations with all temperature parameters for spring and summer as well as the minimum autumn temperature, and H/W also with the mean and maximum autumn temperatures as well as the mean and minimum winter temperatures. The correlation coefficients reached − 0.28, − 0.27 and − 0.28, respectively. These three features significantly correlated with summer and spring humidity, at up to 0.23.The number of whorls (whl), shell width (W), shell diameter (D), demonstrated a similar correlation with climatic parameters (Fig. 4). They showed the largest and significant correlation coefficients with winter temperatures: up to 0.24, 0.22 and 0.22, respectively. The ratio of umbilicus minor to its major diameter (u/U) showed significant positive correlation up to 0.22 with temperature of warmer seasons.The climatic parameters were grouped into several clusters indicating their similar relationships with morphometric features (Fig. 4). Humidity parameters of warmer seasons formed a separate cluster and temperature parameters were grouped according to seasons. The winter parameters were connected with autumn humidity and separated from temperatures for warmer seasons.Modelling relationships between morphometric shell features and climatic parametersThe joint influence of many climatic parameters on morphometry of shells collected in individual years was studied using a linear mixed-effects (LME) model after exclusion of correlated parameters and a linear ridge regression (LRR) model including all climatic parameters. The latter allows for the inclusion of correlated variables. We separately investigated the seasonal maximum, mean and minimum temperature parameters in combination with seasonal humidity parameters (Additional file 2: Table S6) because they are obviously correlated.Umbilicus minor (u) and major (U) diameters as well as umbilicus relative diameter (U/D) turned out best explained by the climatic parameters (with R2  > 0.15) in two models (Additional file 2: Table S6). Moreover, u, U and U/D were described in LME models by the largest number of significant climatic parameters, i.e., 15. The features u and U had also the largest number of significant parameters in LRR models, i.e., 18 out of 24 possibilities. The largest average values of temperature coefficients for the LRR models were 0.66 for D, 0.58 for W, 0.32 for H, 0.26 for U and 0.22 for u. Thus, all the above-mentioned features were under the strongest influence of the climatic conditions.In the case of LRR models, the coefficients at the winter mean temperature were most often selected as significant, in 12 out of 14 possibilities (Additional file 2: Table S6). The humidity coefficients for autumn were significant in 30 cases of 42 possibilities. The highest average absolute values of coefficients in climatic variables were those for the summer (0.63), spring (0.31) and autumn (0.24) minimum temperatures as well as the summer mean temperature (0.31). Thus, the temperatures of warmer seasons were more important for developing shell morphology. Seasonal humidity coefficients showed similar values compared to each other.Comparison of shell morphometry of T. hispidus and T. sericeus kept under various conditionsIn order to verify the influence of different climatic parameters on Trochulus shell morphometry in selected conditions, we compared shells from three groups of T. hispidus, which represented several subsequent generations: (1) parental snails collected in the wild in Wrocław-Jarnołtów, (2) their offspring bred in the laboratory for two generations and (3) offspring of the second laboratory-bred generation transplanted again into a garden in Wrocław (Fig. 5A–C). The comparison of the group 2 and 1 was to verify if laboratory conditions with controlled temperature and humidity can influence the shell morphometry within only one generation, whereas including the group 3 in the comparison, we wanted to check if snails raised in wild garden conditions can recover the original phenotype. Furthermore, we transplanted into the same garden conditions T. sericeus, which was collected in the wild in Muszkowice (Fig. 5D,E). In this case, we verified if two originally different ecophenotypes T. hispidus and T. sericeus, develop the same shell morphometry under the same conditions.Figure 5Shells of two Trochulus ecophenotypes: parental T. hispidus from wild habitat in Wrocław (A), the first generation of T. hispidus raised in laboratory (B); T. hispidus reared in garden in Wrocław (C); T. sericeus from wild habitat in Muszkowice (D); T. sericeus reared in garden in Wrocław (E).Full size imageConditions in which these snails developed were different. According to WorldClim, the wild environment of T. hispidus in Wrocław was generally warmer than that of T. sericeus in Muszkowice (Additional file 2: Table S7). The largest difference was 1.4 °C for the maximum summer temperature. Relative humidity was lower in Wrocław by up to 2% for warmer seasons but was higher in winter by 1.6%. The difference between the wild and garden localities in Wrocław was much smaller and did not exceed 0.41 °C. The garden conditions were less humid, by up to 2%. However, data from WorldClim are generalized over a longer period and wider regions, so may not well reflect local conditions in the studied places. Actually, the Wrocław site was an open habitat covered with a nettle community like a garden patch, while the Muszkowice site was overgrown by a beech forest, which most likely maintained a higher humidity and a more stable temperature.Laboratory temperatures were substantially different from those in the field, especially for winter (by 18–19.7 °C) as well as for spring and autumn (by 8.2–12 °C). Laboratory humidity was by up to 4.5% lower compared to winter and 5.9–9.9% higher than in spring and summer.A discriminant function analysis (DFA) for the defined groups of snails provided their interesting grouping and separation (Fig. 6). The analysis identified three significant discriminant functions (p  More

FAO. The State of World Fisheries and Aquaculture 2020., https://doi.org/10.4060/ca9229en (FAO, 2020).Watson, J. W. & Kerstetter, D. W. Pelagic Longline Fishing Gear: A Brief History and Review of Research Efforts to Improve Selectivity. Mar. Technol. Soc. J. 40, 6–11 (2006).Article 

Google Scholar 
Hare, S. R. et al. The western an our d central Pacific tuna fishery: 2019 overview and status of stocks. (SPC, 2020).Wang, S.-P. Stock assessment of blue marlin in the Indian Ocean using Stock Synthesis. (IOTC, 2019).Ohshimo, S. et al. Horizontal distribution and habitat of Pacific bluefin tuna, Thunnus orientalis, larvae in the waters around Japan. Bull. Mar. Sci. 93, 769–787 (2017).ADS 
Article 

Google Scholar 
Margulies, D., Scholey, V. P., Wexler, J. B. & Stein, M. S. Chapter 5 – Research on the Reproductive Biology and Early Life History of Yellowfin Tuna Thunnus albacares in Panama. In Advances in Tuna Aquaculture: From Hatchery to Market (eds. Benetti, D. D., Partridge, G. J. & Buentello, A.) 77–114, https://doi.org/10.1016/B978-0-12-411459-3.00004-7 (Academic Press, 2016).Madigan, D. J. et al. Intrinsic tracers reveal recent foraging ecology of giant Pacific bluefin tuna at their primary spawning grounds. Mar. Ecol. Prog. Ser. 553, 253–266 (2016).ADS 
CAS 
Article 

Google Scholar 
Ward, T. M., Staunton-Smith, J., Hoyle, S. & Halliday, I. A. Spawning patterns of four species of predominantly temperate pelagic fishes in the sub-tropical waters of southern Queensland. Estuar. Coast. Shelf Sci. 56, 1125–1140 (2003).ADS 
Article 

Google Scholar 
Kiyofuji, H. & Ochi, D. Proposal of alternative spatial structure for skipjack stock assessment in the WCPO. (WCPFC, 2016).Reglero, P., Tittensor, D., Álvarez-Berastegui, D., Aparicio-González, A. & Worm, B. Worldwide distributions of tuna larvae: revisiting hypotheses on environmental requirements for spawning habitats. Mar. Ecol. Prog. Ser. 501, 207–224 (2014).ADS 
Article 

Google Scholar 
Schilling, H. T. et al. Multiple spawning events promote increased larval dispersal of a predatory fish in a western boundary current. Fish. Oceanogr. 29, 309–323 (2020).Article 

Google Scholar 
Richardson, D. E., Hare, J. A., Overholtz, W. J. & Johnson, D. L. Development of long-term larval indices for Atlantic herring (Clupea harengus) on the northeast US continental shelf. ICES J. Mar. Sci. 67, 617–627 (2010).Article 

Google Scholar 
Muhling, B. A. et al. Overlap between Atlantic bluefin tuna spawning grounds and observed Deepwater Horizon surface oil in the northern Gulf of Mexico. Mar. Pollut. Bull. 64, 679–687 (2012).CAS 
Article 

Google Scholar 
Brown, S. K., Buja, K. R., Jury, S. H., Monaco, M. E. & Banner, A. Habitat Suitability Index Models for Eight Fish and Invertebrate Species in Casco and Sheepscot Bays, Maine. North Am. J. Fish. Manag. 20, 408–435 (2000).Article 

Google Scholar 
Hernández, C. M. et al. Evidence and patterns of tuna spawning inside a large no-take Marine Protected Area. Sci. Rep. 9, 1–11 (2019).Article 

Google Scholar 
Janßen, H. et al. Integration of fisheries into marine spatial planning: Quo vadis? Estuar. Coast. Shelf Sci. 201, 105–113 (2018).ADS 
Article 

Google Scholar 
Richardson, A. J. et al. Residency and reproductive status of yellowfin tuna in a proposed large-scale pelagic marine protected area. Aquat. Conserv. Mar. Freshw. Ecosyst. 28, 1308–1316 (2018).Article 

Google Scholar 
CBD. Protected areas and other effective area-based conservation measures (Decision 14/8). https://www.cbd.int/doc/decisions/cop-14/cop-14-dec-08-en.pdf (2018).Gurney, G. G. et al. Biodiversity needs every tool in the box: use OECMs. Nature 595, 646–649 (2021).ADS 
CAS 
Article 

Google Scholar 
Hiltz, E., Fuller, S. & Mitchell, J. Disko Fan Conservation Area: a Canadian case study. Parks 24, 17–30 (2018).Article 

Google Scholar 
IUCN-WCPA Task Force on OECMs. Recognising and reporting other effective area-based conservation measures. https://portals.iucn.org/library/sites/library/files/documents/PATRS-003-En.pdf (2019).Berkeley, S. A., Hixon, M. A., Larson, R. J. & Love, M. S. Fisheries Sustainability via Protection of Age Structure and Spatial Distribution of Fish Populations. Fisheries 29, 23–32 (2004).Article 

Google Scholar 
Hall, S. In A Fishery manager’s guidebook – Second Edition (eds. Garcia, S. M. & Cochrane, K. L.) 196–219 Ch.8 – Area and time restrictions. (Wiley-Blackwell, 2009).Jonas, H. D., Barbuto, V., Jonas, H. C., Kothari, A. & Nelson, F. New steps of change: looking beyond protected areas to consider other effective area-based conservation measures. Parks 20, 111–128 (2014).Article 

Google Scholar 
Dunn, D. C., Maxwell, S. M., Boustany, A. M. & Halpin, P. N. Dynamic ocean management increases the efficiency and efficacy of fisheries management. Proc. Natl. Acad. Sci. 113, 668–673 (2016).ADS 
CAS 
Article 

Google Scholar 
Le Quesne, W. J. F. & Codling, E. A. Managing mobile species with MPAs: the effects of mobility, larval dispersal, and fishing mortality on closure size. ICES J. Mar. Sci. 66, 122–131 (2009).Article 

Google Scholar 
Richardson, D. E. et al. Discovery of a spawning ground reveals diverse migration strategies in Atlantic bluefin tuna (Thunnus thynnus). Proc. Natl. Acad. Sci. 113, 3299–3304 (2016).ADS 
CAS 
Article 

Google Scholar 
Dueri, S., Bopp, L. & Maury, O. Projecting the impacts of climate change on skipjack tuna abundance and spatial distribution. Glob. Change Biol. 20, 742–753 (2014).ADS 
Article 

Google Scholar 
Asch, R. G. Climate change and decadal shifts in the phenology of larval fishes in the California Current ecosystem. Proc. Natl. Acad. Sci. 112, E4065–E4074 (2015).ADS 
CAS 
Article 

Google Scholar 
Mountain, D. G. & Murawski, S. A. Variation in the distribution of fish stocks on the northeast continental shelf in relation to their environment, 1980–1989. ICES mar. Sci. Symp. 195, 424–432.Muhling, B. A. et al. Potential impact of climate change on the Intra-Americas Sea: Part 2. Implications for Atlantic bluefin tuna and skipjack tuna adult and larval habitats. J. Mar. Syst. 148, 1–13 (2015).Article 

Google Scholar 
Asch, R. G., Stock, C. A. & Sarmiento, J. L. Climate change impacts on mismatches between phytoplankton blooms and fish spawning phenology. Glob. Change Biol. 25, 2544–2559 (2019).ADS 
Article 

Google Scholar 
Pankhurst, N. W. & Munday, P. L. Effects of climate change on fish reproduction and early life history stages. Mar. Freshw. Res. 62, 1015–1026 (2011).CAS 
Article 

Google Scholar 
Lehodey, P., Senina, I., Nicol, S. & Hampton, J. Modelling the impact of climate change on South Pacific albacore tuna. Deep Sea Res. Part II Top. Stud. Oceanogr. 113, 246–259 (2015).ADS 
CAS 
Article 

Google Scholar 
Muhling, B. A. et al. Collection of Larval Bluefin Tuna (Thunnus Thynnus) Outside Documented Western Atlantic Spawning Grounds. Bull. Mar. Sci. 87, 687–694 (2011).Article 

Google Scholar 
Nishikawa, Y., Honma, M., Ueyanagi, S. & Kikawa, S. Average Distribution of Larvae of Oceanic Species of Scombroid Fishes, 1956–1981. (Far Seas Fisheries Research Laboratory, 1985).Maunder, M. N. & Punt, A. E. Standardizing catch and effort data: a review of recent approaches. Fish. Res. 70, 141–159 (2004).Article 

Google Scholar 
Collette, B. B. In Annotated Checklist of Fishes Vol. 19 Family Scombridae Rafinesque 1815 – mackerels, tunas, and bonitos. (California Academy of Sciences, 2003).Fricke, R., Eschmeyer, W. N. & Van der Laan, R. Eschmeyer’s catalog of fishes; Genera, species, references. Electronic version. (California Academy of Sciences, 2022).R Core Team. R: A language and environment for statistical computing. (R Foundation for Statistical Computing, 2022).Buenafe, KCV. tinbuenafe/DigitizingNishikawa: Digitizing Nishikawa v3.0, Zenodo, https://doi.org/10.5281/zenodo.6592148 (2022).Hijmans, R. J. et al. terra: Spatial Data Analysis. (2022).Pebesma, E. Simple Features for R: Standardized Support for Spatial Vector Data. R J. 10, 439–446 (2018).Article 

Google Scholar 
Blondel, E. & Billet, N. RFigisGeo: A R package to handle utilities for geospatial processing. (2022).FAO. FAO Major Fishing Areas. https://www.fao.org/fishery/en/area/search (2022).Richards, W. J. & Potthoff, T. Analysis of the Taxonomic Characters of Young Scombrid Fishes, Genus Thunnus. In The Early life history of Fish (ed. Blaxter, J. H. S.) 623–648, https://doi.org/10.1007/978-3-642-65852-5_50 (Springer, 1974).Luthy, S. A., Cowen, R. K., Serafy, J. E. & McDowell, J. R. Toward identification of larval sailfish (Istiophorus platypterus), white marlin (Tetrapturus albidus), and blue marlin (Makaira nigricans) in the western North Atlantic Ocean. Fish. Bull. 103 (2004). More

Berthold, P. Bird Migration: A General Survey. (Oxford University Press, 2001).Harrison, X. A., Blount, J. D., Inger, R., Norris, D. R. & Bearhop, S. Carry-over effects as drivers of fitness differences in animals. J. Anim. Ecol. 80, 4–18 (2011).PubMed 
Article 

Google Scholar 
Webster, M. S., Marra, P. P., Haig, S. M., Bensch, S. & Holmes, R. T. Links between worlds: Unraveling migratory connectivity. Trends Ecol. Evol. 17, 76–83 (2002).Article 

Google Scholar 
Saurola, P., Valkama, J. & Velmala, W. Suomen rengastusatlas Osa I/The Finnish Bird Ringing Atlas Vol. I. (Finnish Museum of Natural History and Ministry of Environment, 2013).Bergman, G. Allin ja mustalinnun muuttokannat keväällä 1960 (in Finnish). Suomen Riista 14, 69–74 (1961).
Google Scholar 
Skov, H. et al. Waterbird Populations and Pressures in the Baltic Sea. (TemaNord 550, 2011).Grenquist, P. Öljytuhoista Suomen aluevesillä v. 1948–1955. Suomen Riista 10, 105–116 (1956).Hario, M., Rintala, J. & Nordenswan, G. Dynamics of wintering long-tailed ducks in the Baltic Sea–the connection with lemming cycles, oil disasters, and hunting. Suomen Riista 55, 83–96 (2009).
Google Scholar 
Ellermaa, M. & Pettay, T. Põõsaspean niemen arktinen muutto syksyllä 2004. Linnut Vuosik. 2005, 99–112 (2005).
Google Scholar 
Delany, S. & Scott, D. Waterbird Population Estimates. (Wetlands International, 2006).Nolet, B. A. et al. Faltering lemming cycles reduce productivity and population size of a migratory Arctic goose species. J. Anim. Ecol. 82, 804–813 (2013).PubMed 
PubMed Central 
Article 

Google Scholar 
Sokolov, V., Vardeh, S. & Quillfeldt, P. Long-tailed Duck (Clangula hyemalis) ecology: Insights from the Russian literature. Part 1: Asian part of the Russian breeding range. Polar Biol. 42, 2259–2276 (2019).Article 

Google Scholar 
Summers, R. W. & Underhill, L. G. Factors related to breeding production of Brent Geese Branta b. bernicla and waders (Charadrii) on the Taimyr Peninsula. Bird Study 34(161), 171 (1987).
Google Scholar 
Summers, R. W., Underhill, L. G. & Syroechkovski, J. The breeding productivity of dark-bellied brent geese and curlew sandpipers in relation to changes in the numbers of arctic foxes and lemmings on the Taimyr Peninsula Siberia. Ecography 21, 573–580 (1998).Article 

Google Scholar 
Underhill, L. G. et al. Breeding of waders (Charadrii) and Brent Geese Branta bernicla bernicla at Pronchishcheva Lake, northeastern Taimyr, Russia, in a peak and a decreasing lemming year. Ibis 135, 277–292 (1993).Article 

Google Scholar 
Gauthier, G., Bëty, J., Giroux, J.-F. & Rochefort, L. Trophic interactions in a High Arctic snow goose colony. Integr. Comp. Biol. 44, 119–129 (2004).PubMed 
Article 

Google Scholar 
Elton, C. Voles, Mice and Lemmings: Problems in Population Dynamics. (Clarendon Press, 1942).Ehrich, D. et al. Documenting lemming population change in the Arctic: Can we detect trends?. Ambio https://doi.org/10.1007/s13280-019-01198-7 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
Kokorev, Y. I. & Kuksov, V. A. Population dynamics of lemmings, Lemmus sibirica and Dicrostonyx torquatus, and Arctic Fox Alopex lagopus on the Taimyr peninsula, Siberia, 1960–2001. Ornis Svecica 12, 139–145 (2002).
Google Scholar 
Angerbjörn, A., Tannerfeldt, M. & Erlinge, S. Predator-prey relationships: Arctic foxes and lemmings. J. Anim. Ecol. 68, 34–49 (1999).Article 

Google Scholar 
Fauteux, D., Gauthier, G. & Berteaux, D. Seasonal demography of a cyclic lemming population in the Canadian Arctic. J. Anim. Ecol. 84, 1412–1422 (2015).PubMed 
Article 

Google Scholar 
Gilg, O., Sittler, B. & Hanski, I. Climate change and cyclic predator–prey population dynamics in the high Arctic. Glob. Chang. Biol. 15, 2634–2652 (2009).ADS 
Article 

Google Scholar 
Berryman, A. A. The orgins and evolution of predator-prey theory. Ecology 73, 1530–1535 (1992).Article 

Google Scholar 
Framstad, E., Stenseth, N. C., Bjørnstad, O. N. & Falck, W. Limit cycles in Norwegian lemmings: Tensions between phase-dependence and density-dependence. Proc. R Soc. London. Ser. B Biol. Sci. 264, 31–38 (1997).ADS 
Article 

Google Scholar 
Hanski, I. & Korpimaki, E. Microtine rodent dynamics in northern Europe: Parameterized models for the predator-prey interaction. Ecology 76, 840–850 (1995).Article 

Google Scholar 
May, R. M. Limit cycles in predator-prey communities. Science 177, 900–902 (1972).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Gilg, O., Hanski, I. & Sittler, B. Cyclic dynamics in a simple vertebrate predator-prey community. Science 302, 866–868 (2003).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Juhasz, C. C., Shipley, B., Gauthier, G., Berteaux, D. & Lecomte, N. Direct and indirect effects of regional and local climatic factors on trophic interactions in the Arctic tundra. J. Anim. Ecol. 89, 704–715 (2020).PubMed 
Article 

Google Scholar 
McKinnon, L., Berteaux, D., Gauthier, G. & Bêty, J. Predator-mediated interactions between preferred, alternative and incidental prey in the arctic tundra. Oikos 122, 1042–1048 (2013).Article 

Google Scholar 
Angelstam, P., Lindström, E. & Widén, P. Role of predation in short-term population fluctuations of some birds and mammals in Fennoscandia. Oecologia 62, 199–208 (1984).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Ehrich, D. et al. Vole abundance and reindeer carcasses determine breeding activity of Arctic foxes in low Arctic Yamal Russia. BMC Ecol. 17, 1–13 (2017).Article 

Google Scholar 
Brook, R. W., Duncan, D. C., Hines, J. E., Carrière, S. & Clark, R. G. Effects of small mammal cycles on productivity of boreal ducks. Wildlife Biol. 11, 3–11 (2005).Article 

Google Scholar 
Guillemain, M. et al. Effects of climate change on European ducks: what do we know and what do we need to know?. Wildlife Biol. 19, 404–419 (2013).Article 

Google Scholar 
Pehrsson, O. Duckling production of the Oldsquaw in relation to spring weather and small-rodent fluctuations. Can. J. Zool. 64, 1835–1841 (1986).Article 

Google Scholar 
ACIA. Impacts of a Warming Arctic: Arctic Climate Impact Assessment. (Cambridge University Press, 2004).Høye, T. T., Post, E., Meltofte, H., Schmidt, N. M. & Forchhammer, M. C. Rapid advancement of spring in the High Arctic. Curr. Biol. 17, R449–R451 (2007).PubMed 
Article 
CAS 

Google Scholar 
Post, E. et al. Ecological dynamics across the Arctic associated with recent climate change. Science 325, 1355–1358 (2009).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Kausrud, K. L. et al. Linking climate change to lemming cycles. Nature 456, 93–97 (2008).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Berteaux, D. et al. Effects of changing permafrost and snow conditions on tundra wildlife: Critical places and times. Arct. Sci. 3, 65–90 (2017).Article 

Google Scholar 
Bilodeau, F., Gauthier, G. & Berteaux, D. The effect of snow cover on lemming population cycles in the Canadian High Arctic. Oecologia 172, 1007–1016 (2013).ADS 
PubMed 
Article 

Google Scholar 
Madsen, F. J. On the food habits of the diving ducks in Denmark. Danish Rev. Game Biol. 3, 2–83 (1954).
Google Scholar 
Nilsson, L. Habitat selection, food choice, and feeding habits of diving ducks in coastal waters of South Sweden during the non-breeding season. Ornis Scand. 3, 55–78 (1972).Article 

Google Scholar 
Žydelis, R. & Ruškytė, D. Winter foraging of long-tailed ducks (Clangula hyemalis) exploiting different benthic communities in the Baltic Sea. Wilson Bull. 117, 133–141 (2005).Article 

Google Scholar 
Skabeikis, A. et al. Effect of round goby (Neogobius melanostomus) invasion on blue mussel (Mytilus edulis trossulus) population and winter diet of the long-tailed duck (Clangula hyemalis). Biol. Invasions 21, 911–923 (2019).Article 

Google Scholar 
Laursen, K. & Møller, A. P. Long-Term changes in nutrients and mussel stocks are related to numbers of breeding eiders Somateria mollissima at a large Baltic colony. PLoS ONE 9, e95851 (2014).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Carstensen, J., Andersen, J. H., Gustafsson, B. G. & Conley, D. J. Deoxygenation of the baltic sea during the last century. Proc. Natl. Acad. Sci. USA 111, 5628–5633 (2014).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Savchuk, O. P. Large-scale nutrient dynamics in the Baltic Sea, 1970–2016. Front. Mar. Sci. 5, 95 (2018).Article 

Google Scholar 
Møller, A. P., Flensted-Jensen, E. & Mardal, W. Agriculture, fertilizers and life history of a coastal seabird. J. Anim. Ecol. 76, 515–525 (2007).PubMed 
Article 

Google Scholar 
Møller, A. P., Thorup, O. & Laursen, K. Predation and nutrients drive population declines in breeding waders. Ecol. Appl. 28, 1292–1301 (2018).PubMed 
Article 

Google Scholar 
Gelman, A., Carlin, J. B., Stern, H. S. & Rubin, D. B. Bayesian Data Analysis. (Chapman & Hall/CRC, 2004).Lebreton, J.-D. & Gimenez, O. Detecting and estimating density dependence in wildlife populations. J. Wildl. Manage. 77, 12–23 (2013).Article 

Google Scholar 
Bergman, G. The spring migration of the Long-tailed Duck and the Common Scoter in western Finland. Ornis Fenn. 51, 129–145 (1974).
Google Scholar 
Richardson, W. J. Timing and amount of bird migration in relation to weather: A Review. Oikos 30, 224–272 (1978).Article 

Google Scholar 
Alerstam, T. Bird flight and optimal migration. Trends Ecol. Evol. 6, 210–215 (1991).CAS 
PubMed 
Article 

Google Scholar 
Richardson, W. J. Wind and Orientation of Migrating Birds: A Review. in Orientation in Birds (ed. Berthold, P.) 226–249 (Birkhäuser, 1991). https://doi.org/10.1007/978-3-0348-7208-9_11.Christensen, T. K. & Fox, A. D. Changes in age and sex ratios amongst samples of hunter-shot wings from common duck species in Denmark 1982–2010. Eur. J. Wildl. Res. 60, 303–312 (2014).Article 

Google Scholar 
Fox, A. D., Clausen, K. K., Dalby, L., Christensen, T. K. & Sunde, P. Age-ratio bias among hunter-based surveys of Eurasian Wigeon Anas penelope based on wing vs. field samples. Ibis 157, 391–395 (2015).Article 

Google Scholar 
Møller, A. P., Flensted-Jensen, E., Laursen, K. & Mardal, W. Fertilizer leakage to the marine environment, ecosystem effects and population trends of waterbirds in Denmark. Ecosystems 18, 30–44 (2015).Article 
CAS 

Google Scholar 
Scott, D. A. & Rose, P. M. Atlas of Anatidae Populations in Africa and Western Eurasia. Wetlands International Publication 41 (Wetlands International, 1996).Fick, S. E. & Hijmans, R. J. WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017).Article 

Google Scholar 
Harris, I., Jones, P. D., Osborn, T. J. & Lister, D. H. Updated high-resolution grids of monthly climatic observations–the CRU TS3.10 Dataset. Int. J. Climatol. 34, 623–642 (2014).Article 

Google Scholar 
Hijmans, R. J. Introduction to the ’raster’ package (version 3.0–12). https://rspatial.org/raster/pkg/index.html (2020).National Center for Atmospheric Research Staff. The climate data guide: Hurrell North Atlantic Oscillation (NAO) index (PC-based). https://climatedataguide.ucar.edu/climate-data/hurrell-north-atlantic-oscillation-nao-index-pc-based (2019).Hurrell, J. W. Decadal trends in the north atlantic oscillation: Regional temperatures and precipitation. Science 269, 676–679 (1995).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Büttger, H., Nehls, G. & Stoddard, P. The history of intertidal blue mussel beds in the North Frisian Wadden Sea in the 20th century: Can we define reference conditions for conservation targets by analysing aerial photographs?. J. Sea Res. 87, 91–102 (2014).ADS 
Article 

Google Scholar 
Kristensen, P. S. & Borgstrøm, R. The Danish Wadden Sea: Fishery of mussels (Mytilus edulis L.) in a wildlife reserve? in Proceedings from the 11. Scientific Wadden Sea Symposium, Esbjerg, Denmark, 4.-8. April 2005. NERI technical report (ed. Laursen, K.) vol. 573 107–111 (National Environmental Research Institute. Department of Wildlife Ecology and Biodiversity, 2006).Baird, R. H. Measurement of condition in mussels and oysters. ICES J. Mar. Sci. 23, 249–257 (1958).Article 

Google Scholar 
Waldeck, P. & Larsson, K. Effects of winter water temperature on mass loss in Baltic blue mussels: Implications for foraging sea ducks. J. Exp. Mar. Bio. Ecol. 444, 24–30 (2013).Article 

Google Scholar 
Nehls, G. et al. Beds of blue mussels and Pacific oysters. Quality Status Report, Thematic Report; No. 11. Wadden Sea Ecosystem; No. 25 (2009).Laursen, K., Møller, A. P., Haugaard, L., Öst, M. & Vainio, J. Allocation of body reserves during winter in eider Somateria mollissima as preparation for spring migration and reproduction. J. Sea Res. 144, 49–56 (2019).ADS 
Article 

Google Scholar 
Morelli, F., Laursen, K., Svitok, M., Benedetti, Y. & Møller, A. P. Eiders, nutrients and eagles: Bottom-up and top-down population dynamics in a marine bird. J. Anim. Ecol. https://doi.org/10.1111/1365-2656.13498 (2021).Article 
PubMed 

Google Scholar 
Westerbom, M., Kilpi, M. & Mustonen, O. Blue mussels, Mytilus edulis, at the edge of the range: population structure, growth and biomass along a salinity gradient in the north-eastern Baltic Sea. Mar. Biol. 140, 991–999 (2002).Article 

Google Scholar 
Kery, M. & Schaub, M. Bayesian Population Analysis Using WinBUGS: A Hierarchical Perspective. (Elsevier, 2012).Kerman, J. Neutral noninformative and informative conjugate beta and gamma prior distributions. Electron. J. Stat. 5, 1450–1470 (2011).MathSciNet 
MATH 
Article 

Google Scholar 
Crainiceanu, C. M., Ruppert, D. & Wand, M. P. Bayesian analysis for penalized spline regression using WinBUGS. J. Stat. Softw. 14, (2005).Saha, K. & Paul, S. Bias-corrected maximum likelihood estimator of the negative binomial dispersion parameter. Biometrics 61, 179–185 (2005).MathSciNet 
PubMed 
MATH 
Article 

Google Scholar 
Mutshinda, C. M., O’Hara, R. B. & Woiwod, I. P. A multispecies perspective on ecological impacts of climatic forcing. J. Anim. Ecol. 80, 101–107 (2011).PubMed 
Article 

Google Scholar 
Pöysä, H. et al. Environmental variability and population dynamics: Do European and North American ducks play by the same rules?. Ecol. Evol. 6, 7004–7014 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
Almaraz, P., Green, A. J., Aguilera, E., Rendón, M. A. & Bustamante, J. Estimating partial observability and nonlinear climate effects on stochastic community dynamics of migratory waterfowl. J. Anim. Ecol. https://doi.org/10.1111/j.1365-2656.2012.01972.x (2012).Article 
PubMed 

Google Scholar 
Schmidt, N. M. et al. Response of an arctic predator guild to collapsing lemming cycles. Proc. R. Soc. B Biol. Sci. 279, 4417–4422 (2012).Article 

Google Scholar 
Ebbinge, B. S., Heesterbeek, H. J. A. P., Ens, B. J. & Goedhart, P. W. Density dependent population limitation in dark-bellied brent geese Branta b. bernicla. Avian Sci. 2, 63–75 (2002).
Google Scholar 
Domine, F. et al. Snow physical properties may be a significant determinant of lemming population dynamics in the high Arctic. Arct. Sci. 4, 813–826 (2018).Article 

Google Scholar 
Ims, R. A., Henden, J.-A. & Killengreen, S. T. Collapsing population cycles. Trends Ecol. Evol. 23, 79–86 (2008).PubMed 
Article 

Google Scholar 
Korslund, L. & Steen, H. Small rodent winter survival: Snow conditions limit access to food resources. J. Anim. Ecol. 75, 156–166 (2006).PubMed 
Article 

Google Scholar 
Callaghan, T. V. et al. The changing face of Arctic snow cover: A synthesis of observed and projected changes. Ambio 40, 17–31 (2011).Article 

Google Scholar 
Machín, P. et al. The role of ecological and environmental conditions on the nesting success of waders in sub-Arctic Sweden. Polar Biol. 42, 1571–1579 (2019).Article 

Google Scholar 
Koneff, M. D. et al. Evaluation of harvest and information needs for North American sea ducks. PLoS ONE 12, e0175411 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Benton, T. G. & Grant, A. Elasticity analysis as an important tool in evolutionary and population ecology. Trends Ecol. Evol. 14, 467–471 (1999).CAS 
PubMed 
Article 

Google Scholar 
Heppell, S. S., Caswell, H. & Crowder, L. B. Life histories and elasticity patterns: Perturbation analysis for species with minimal demographic data. Ecology 81, 654–665 (2000).Article 

Google Scholar 
Sæther, B.-E. & Bakke, O. Avian life history variation and contribution of demographic traits to the population growth rate. Ecology 81, 642–653 (2000).Article 

Google Scholar 
Öst, M., Ramula, S., Lindén, A., Karell, P. & Kilpi, M. Small-scale spatial and temporal variation in the demographic processes underlying the large-scale decline of eiders in the Baltic Sea. Popul. Ecol. 58, 121–133 (2016).Article 

Google Scholar 
Holopainen, S. & Fox, A. D. Associations between duck harvest, hunting wing ratios and measures of reproductive output in Northern Europe. Eur. J. Wildl. Res. 64, (2018).Conley, D. J., Humborg, C., Rahm, L., Savchuk, O. P. & Wulff, F. Hypoxia in the Baltic Sea and basin-scale changes in phosphorus biogeochemistry. Environ. Sci. Technol. 36, 5315–5320 (2002).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Carstensen, J. et al. Hypoxia in the Baltic Sea: Biogeochemical cycles, benthic fauna, and management. Ambio 43, 26–36 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Conley, D. J. et al. Hypoxia-related processes in the Baltic Sea. Environ. Sci. Technol. 43, 3412–3420 (2009).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Conley, D. J. et al. Long-term changes and impacts of hypoxia in Danish coastal waters. Ecol. Appl. 17, 165–184 (2007).Article 

Google Scholar 
Diaz, R. J. & Rosenberg, R. Spreading dead zones and consequences for marine ecosystems. Science 321, 926–929 (2008).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Fox, A. D. et al. Current and potential threats to Nordic duck populations–a horizon scanning exercise. Ann. Zool. Fennici 52, 193–220 (2015).Article 

Google Scholar 
Møller, A. P. Biological consequences of global change for birds. Integr. Zool. 8, 136–144 (2013).PubMed 
Article 

Google Scholar  More

All figures were produced using the R package ggplot2 v3.3.534.eBird and covariate dataeBird data are structured as follows. Birders submit observations as species checklists with counts of each species they identify. They report associated metadata, such as location, date and time, duration of the observation period, number of observers, and sampling protocol25,26,31. The birder indicates whether their checklist is “complete”; complete checklists yield inferred zeroes for all species not reported on a checklist.We retrieved the eBird Basic Dataset containing all eBird observations and sampling metadata. We extracted all complete checklists that occurred within the U.S. state of California between April 1 and June 30, 2019. Four survey-level covariates were retrieved from eBird checklist metadata as detection covariates: number of observers, checklist duration, date of year, and time of day; any checklist that failed to report one or more of these variables was dropped. Corresponding to best practices for use of eBird data, we filtered the data for quality according to the following criteria: we discarded checklists other than those following the “Stationary” survey protocol (observations made at a single spatial location) with duration shorter than 4 hours and at most 10 observers in the group31,35.We selected twenty circular regions of high sampling intensity with 10 km radii across California (Supplemental Fig. 1). These spanned the state’s many habitats including coastal, agricultural, wetland, and mountain areas, and contained active birding areas such as parks and human population centers. In each subregion, we selected 10 species with the highest reporting rate (proportion of checklists including that species) and 10 representing an intermediate reporting rate. An additional 10 species were selected that were detected in many regions to enable cross-region comparisons, yielding 407 species-subregion (SSR) datasets (with overlaps between the two species selection protocols; see Supplemental Section 2 for the full algorithm). Across 20 subregions, we accepted 6094 eBird checklists for analysis, each with an associated count (potentially zero) for each species. Observations were aggregated to sampling sites defined by a 50 m spatial grid. The 50 m grid was chosen to conservatively identify related surveys and was not motivated by biological processes, nor does it represent the sampling area of each survey. In this context, the concept of “closure” in the latent state is already suspect due to the fact that eBird checklist sampling areas are inconsistent. Data were processed in R using the ‘auk’ package36,37.An elevation surface for the state of California was retrieved from WorldClim at (8.3 times 10^{-3}) decimal degrees resolution using the R package raster38,39. This commonly used covariate was included as a baseline spatial covariate to enable comparison of estimation properties across sites, but its biological relevance to abundance is not crucial to our analysis31. Land cover data were retrieved from the LandFire GIS database’s Existing Vegetation Type layer40. For each unique survey location, a 500 m buffer was calculated around the reported location, and the percent of the buffer which was water, tree cover, agriculture or other vegetation (shrub or grassland) was calculated. We used the following five site-level covariates: elevation, and percent of the landscape within a 500 m buffer of the site that was water, trees, agricultural land, or other vegetation. We included six checklist-level covariates: duration, number of observers, time of day, time of day squared, Julian date, and Julian date squared. Covariates were dropped in datasets where only a single unique value was observed for that covariate.Model implementation and selectionWe considered four variants of the N-mixture model and two variants of the GLMM comprising a total of 6 distinct models, defined by the distributions used in the model or sub-model.The GLMM for count data that we considered is defined as$$begin{array}{*{20}l} {y_{{ij}} sim D(mu _{{ij}} ,[theta ])} hfill \ {log (mu _{{ij}} ) = beta _{0} + {mathbf{x}}_{{ij}}^{T} user2{beta } + alpha _{i} } hfill \ {alpha _{i} sim {mathcal{N}}(0,sigma _{alpha } )} hfill \ end{array}$$where (y_{ij}) is th jth observation at site i, D is a probability distribution (which may contain an extra parameter (theta) to account for overdispersion), (mu _{ij}) represents the mean expected count and is a logit-linear combination of observed site- and observation-level covariates (x_{ij}), (beta) are coefficients representing the effect of those covariates, (beta _0) is a log-scale intercept corresponding to the expected log count at the mean site (i.e. with all centered covariates set to 0), and (alpha _i) is the random effect of site i following a normal distribution. Due to the right skew of (exp (y_{ij})), by log-normal distribution theory the log of the expected count at the mean site is (beta _0 + 0.5 sigma _{alpha }^2). We considered two forms of this model, where D was either a Poisson or a negative binomial distribution, in the latter case with the extra parameter (theta).The N-mixture model is defined as$$begin{array}{*{20}l} {y_{{ij}} sim D_{w} (N_{i} ,p_{{ij}} ,[theta _{w} ])} hfill \ {N_{i} sim D_{b} (lambda _{i} ,[theta _{b} ])} hfill \ {{text{logit}}(p_{{ij}} ) = {text{}}{text{logit}}(p_{0} ) + {mathbf{x}}_{{ij(w)}} {mathbf{beta }}_{w} } hfill \ {log (lambda _{i} ) = log (lambda _{0} ) + {mathbf{x}}_{{i(b)}} {mathbf{beta }}_{b} } hfill \ {p_{0} = e^{{frac{{phi _{1} + phi _{2} }}{2}}} } hfill \ {lambda _{0} = e^{{frac{{phi _{1} – phi _{2} }}{2}}} } hfill \ end{array}$$where (D_b) and (D_w) are probability distributions representing between- and within-site variation, respectively; (N_i) is a site-level latent variable normally representing the “true” abundance at site i; (p_{ij}) is the detection probability of each individual on the jth observation event at site i; (lambda _i) is the mean abundance at site i; and (x_{(w)}) and (x_{(b)}) are covariate vectors for detection and abundance, respectively, with corresponding coefficients (beta _w) and (beta _b). For reasons described below, we reparameterize the intercept parameters of the N-mixture submodels, (log (lambda _0)) and (text{ logit }(p_0)), in terms of two orthogonal parameters (phi _1 = log (lambda _0 p_0)) and (phi _2=log (p_0 / lambda _0)). Now (phi _1) and (phi _2) represent the expected log count and the contrast between detection and abundance, respectively, at the mean site. This parameterization allows us to investigate stability of parameter estimation. The log-scale expected count of the N-mixture model is (phi _1 = log (lambda _0 p_0)), analogous to (beta _0 + 0.5 sigma _{alpha }^2) in the GLMM (see Supplemental Section 6). Each submodel distribution D could include or not include an overdispersion parameter ((theta _w) and (theta _b)), yielding four possible N-mixture model variants: binomial-Poisson (B-P), binomial-negative binomial (B-NB), beta-binomial-Poisson (BB-P), and beta-binomial-negative binomial (BB-NB)8,11.We chose to fit models with maximum likelihood estimation (MLE) for computational feasibility and because key diagnostic tools, such as AIC and methods for checking goodness-of-fit and autocorrelation, were best suited to MLE estimation15. We fit N-mixture models with the nimble and nimbleEcology R packages starting with a conservatively large choice of K, the truncation value of the infinite sum in the N-mixture likelihood calculation33,41 (see Supplemental Section 4 for a discussion of maximum likelihood estimation with NIMBLE). We fit GLMMs with the R package glmmTMB42. We applied forward AIC selection to choose the best covariates for each model with each dataset (illustrated in Fig. S1). One spatial covariate (elevation) and two checklist metadata covariates (duration and number of observers) were treated as a priori important and were included in all models. In the N-mixture model, checklist-specific sampling metadata were only allowed in the detection submodel, while land cover covariates and the interactions between them were allowed in both the detection and abundance submodels. Interactions were dropped in datasets when interaction values showed a correlation of > 0.8 with one of their first-order terms. In N-mixture models, additions to both submodels were considered simultaneously during forward AIC selection.For comparisons between models, we selected a heuristic threshold of (Delta text {AIC} > 2) to say that one model is supported over another30.Fit, estimation, and computationGoodness-of-fitWe used the Kolmogorov-Smirnov (KS) test, a p-value based metric, to evaluate goodness-of-fit on each selected model. For GLMMs, residuals were obtained using the DHARMa R package’s ‘simulateResiduals’ and the KS test was applied using the ‘testUniformity’ function43. For N-mixture models, we considered the site-sum randomized quantile (SSRQ) residuals described by Knape et al.15, computing these for each N-mixture model and running a KS test against the normal CDF. We assumed that covariate effects did not vary by space within subregions and chose not to use spatially explicit models31,44. To test this assumption, we applied Moran’s I test to the SSRQ or DHARMa-generated residuals for each site or observation.Parameter estimationWe compared two abundance parameters of interest across models: coefficients for elevation and log expected count at a standard site (in the GLMM, (beta _0 + 0.5 sigma _alpha ^2); in the N-mixture model, (log (lambda _0 p_0))). We examined absolute differences in point estimates and the log-scale ratios between their standard errors.Stability of estimated parametersAttempting to decompose the expected value of observed data into within- and between-site components can lead to ridged likelihood surfaces with difficult-to-estimate optima. Kéry found that instability of model estimates with increasing K occurred when there was a likelihood tradeoff between detection and abundance, resulting in a tendency in abundance toward positive infinity restrained only by K10. Dennis et al. showed that N-mixture models could in fact yield estimates of absolute abundance at infinity18. We interpreted this as a case of a boundary parameter estimate rather than non-identifiability and explored it by reparametrizing as follows. We estimated the intercepts for detection and abundance with two orthogonal parameters (rotated in log space) (phi _1 = log (lambda _0 p_0)) and (phi _2 = log (p_0 / lambda _0)), where (lambda _0) and (p_0) are real-scale abundance and detection probability at the mean site. We hypothesized that in unstable cases, (phi _1), log expected count, is well-informed by the data, but (phi _2), the contrast between abundance and detection, is not well-informed, corresponding to a likelihood ridge as (phi _2 rightarrow -infty) due to detection probability approaching 0 and abundance approaching infinity. This reparameterization isolates the likelihood ridge to one parameter direction, similar to a boundary estimate as (exp (phi _2) rightarrow 0). Boundary estimates occur in many models and are distinct from non-identifiability in that they result from particular datasets. Confidence regions extending from a boundary estimate may include reasonable parameters, reflecting that there is information in the data. We defined a practical lower bound for (phi _2). When (phi _2) was estimated very near that bound, we conditioned on that boundary for (phi _2) when estimating confidence regions for other parameters.In the N-mixture case, diagnosing a boundary estimate for (phi _2) is made more difficult by the need to increase K for large negative (phi _2) to calculate the likelihood accurately. We used an approach like that of Dennis et al.18 to numerically diagnose unstable cases. For each N-mixture variant in each SSR, the final model was refitted twice, using values of K 2000 and 4000 greater than the initial choice. Estimates were considered unstable if the absolute value of the difference in AIC between these two large-K refits was above a tolerance of 0.1. We monitored whether MLE estimates of (phi _1) and (phi _2) also varied with increasing K.Evaluating the fast N-mixture calculationWe extended previous work by Meehan et al. to drastically improve the efficiency of N-mixture models using negative binomial or beta-binomial distributions in submodels45 (see Supplemental Section 3).We ran benchmarks of this likelihood calculation for a single site against the traditional algorithm, which involves iterating over values of N to compute a truncated infinite sum. We calculated the N-mixture likelihood at 5,000 sites and compared the computation time between the two methods for all four N-mixture model variations. We ran benchmarks along gradients of (text {length}(y_i)) (number of replicate observations at the simulated site) and K (the upper bound of the truncated infinite sum) for each variant. More

Nash, R. D. M., Valencia, A. H. & Geffen, A. J. The origin of Fulton’s condition factor—setting the record straight. Fisheries 31(5), 236–238 (2006).
Google Scholar 
Tarkan, A. S., Gaygusuz, Ö., Acıpınar, H., Gürsoy, Ç. & Özuluğ, M. Length–weight relationships of fishes from the Marmara region (NW-Turkey). J. Appl. Ichthyol. 22(4), 271–273 (2006).Article 

Google Scholar 
Al Nahdi, A., de Leaniz, C. G. & King, A. J. Spatio-temporal variation in length-weight relationships and condition of ribbonfish Trichiurus lepturus (Linnaeus, 1758): Implications for fisheries. PLoS One 11(8), e0161989 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Froese, R., Tsikliras, A. C. & Stergiou, K. I. Editorial note on weight–length relations of fishes. Acta Ichthyol. Piscat. 41(4), 261–263 (2011).Article 

Google Scholar 
Torres, M. A. et al. Length–weight relationships for 22 crustecans and cephalopods from the Gulf of Cadiz (SW Spain). Aquat. Liv. Resour. 30, 12 (2017).Article 

Google Scholar 
Rocha, F., Guerra, A. & Gonzalez, A. F. A review of reproductive strategies in cephalopods. Biol. Rev. 76, 291–304 (2001).CAS 
PubMed 
Article 

Google Scholar 
Laptikhovsky, V. & Salman, A. On reproductive strategies of the epipelagic octopods of the superfamily Argonautoidea (Cephalopoda: Octopoda). Mar. Biol. 142, 321–326 (2003).Article 

Google Scholar 
Forsythe, J. W. & van Heukelem, W. F. Growth. In Cephalopod Life Cycles (ed. Boyle, P. R.) 135–156 (Academic Press, 1987).
Google Scholar 
Jereb, P., et al. (eds) 2015. Cephalopod biology and fisheries in Europe: II. Species Accounts. ICES Cooperative Research Report No. 325, p. 360.Salman, A. Cephalopod research in the eastern Mediterranean (East of 23°E): A review. Boll. Malacol. 45, 47–59 (2009).
Google Scholar 
Salman, A. & İzmirli, C. Ege Üniversitesi Su Ürünleri Fakültesi Müzesi (ESFM)’nin cephalopod envanteri. EgeJFAS 37(4), 357–361. https://doi.org/10.12714/egejfas.37.4.06) (2020) (in Turkish).Article 

Google Scholar 
Önsoy, B. & Salman, A. Reproductive biology of the common cuttlefish Sepia officinalis L. (Sepiida: Cephalpoda) in the Aegean Sea. Turk. J. Vet. Anim. Sci. 29, 613–619 (2005).
Google Scholar 
Salman, A. Reproductive biology of the elegant cuttlefish (Sepia elegans) in the Eastern Mediterranean. Turk. J. Fish. Aquat. Sci. 15(2), 265–272 (2015).Article 

Google Scholar 
Dursun, D., Eronat, E. G. T., Akalın, M. & Salman, M. A. Reproductive biology of pink cuttlefish Sepia orbignyana in the Aegean Sea (eastern Mediterranean). Turk. J. Zool. 37, 576–581 (2013).Article 

Google Scholar 
Salman, A. Reproductive biology of Sepietta oweniana (Pfeffer, 1908) (Sepiolidae: Cephalopoda) in the Aegean Sea. Sci. Mar. 62(4), 379–383 (1998).Article 

Google Scholar 
Salman, A. & Önsoy, B. Reproductive biology of the bobtail squid Rossia macrosoma (Cephalopoda: Sepiolidea) from the eastern Mediterranean. Turk. J. Fish. Aquat. Sci. 10, 81–86 (2010).Article 

Google Scholar 
Salman, A. Fecundity and spawning strategy of shortfin squid Illex coindetii (Oegopsida: Ommastrephidae), in the eastern Mediterranean. Turk. J. Fish. Aquat. Sci. 17, 841–849 (2017).
Google Scholar 
Mangold-Wirz, K. Biologie des céphalopodes benthiques et nectoniques de la Mer Catalane. Vie Millieu suppl. 13, 1–285 (1963).
Google Scholar 
Salman, A. Fecundity, spawning strategy and oocyte development of shortfin squid Alloteuthis media (Myopsida: Loliginidae) in the eastern Mediterranean. Cah. Biol. Mar. 55, 163–171 (2014).
Google Scholar 
Önsoy, B. & Salman, A. Reproduction patterns of the Mediterranean endemic, Eledone moschata (Lamarck, 1798) (Octopoda: Cephalopoda) in the eastern Mediterranean. (In Turkish) 1st National Malacology Congress, 1–3 September 2004, Izmir-Turkey (Bilal Öztürk & Alp Salman, eds). Turk. J. Aquat. Life 2(2), 55–60 (2004).
Google Scholar 
Tesch, F. W. Age and growth. In Methods for Assessment of Fish Production in Fresh Waters (ed. Ricker, W. E.) 99–130 (Blackwell Scientific Publications, 1971).
Google Scholar 
Merella, P., Quetglas, A., Alemany, F. & Carbonell, A. Length–weight relationship of fishes and cephalopods from the Balearic Islands (western Mediterranean). Naga ICLARM Q. 20(3–4), 66–68 (1997).
Google Scholar 
Manfrin Piccinetti, G. & Giovanardi, O. Données sur la biologie de Sepia officinalis L. dans l’Adriatique obtenues lors de expéditions pipeta. FAO Fish. Rep. 290, 135–138 (1984).
Google Scholar 
Bello, G. Length–weight relationship in males and females of Sepia orbignyana and Sepia elegans (Cephalopoda: Sepiidae). Rapp. Comm. Int. Mer. Médit. 31(2), 254 (1988).
Google Scholar 
Ragonese, S. & Jereb, P. Length-weight relationship and growth of the pink and elegant cuttlefish Sepia orbignyana and Sepia elegans in the Sicilian Channel. In Acta of the 1st International Symposium on the Cuttlefish (ed. Boucaud-Camou, E.) 31–47 (SEPIA. Centre de Publications de l’Universite de Caen, 1991).
Google Scholar 
Akyol, O. & Metin, G. An investigation on determination of some morphological characteristics of Cephalopods in Izmir Bay (Aegean Sea). EU J. Fish. Aquat. Sci. 18(3–4), 357–365 (2001).
Google Scholar 
Lefkaditou, E., Verriopoulos, G. & Valavanis, V. VII9. Research on Cephalopod resources in Hellas. In State of Hellenic Fisheries (eds Papaconstantinou, C. et al.) 440–451 (HCMR Publications, 2007).
Google Scholar 
Duysak, Ö., Sendão, J., Borges, T., Türeli, C. & Erdem, Ü. Cephalopod distribution in Iskenderun bay (eastern Mediterranean–Turkey). J. Fish. Sci. 2, 118–125 (2008).
Google Scholar 
Giordano, D. et al. Distribution and biology of Sepietta oweniana (Pfeffer, 1908) (Cephalopoda: Sepiolidae) in the southern Tyrrhenian Sea (central Mediterranean Sea). Cah. Biol. Mar. 50, 1–10 (2009).
Google Scholar 
Andriguetto, J. M. Jr. & Haimovici, M. Effects of fixation and preservation methods on the morphology of a Loliginid squid (Cephalopoda: Myopsida). Am. Malac Bull. 6(2), 213–217 (1988).
Google Scholar 
Sanchez, P. Donnés preliminaires sur la biologie de trois species de cephalopods de la Mer Catalan. Rapp. Comm. Int. Mer. Médit. 30(2), 247 (1986).
Google Scholar 
Belcari, P., Sartor, P., Nannini, N. & De Ranieri, S. Length-weight relationship of Toda- ropsis eblanae (Cephalopoda: Ommastrephidae) of the northern Tyrrhenian Sea in relation to sexual maturation. Biol. Mar. Mediter. 6, 524–528 (1999).
Google Scholar 
Belcari, P. Length–weight relationship in relation to sexual maturation of Illex coindetii (Cephalopoda: Ommastrephidae) in the northern Tyrrhenian Sea (western Mediterranean). Sci. Mar. 60, 379–384 (1996).
Google Scholar 
Petric, M., Ferri, J., Skeljo, F. & Krstulovic Sifner, S. Body and beak measures of Illex coindetii (Cephalopoda: Ommastrephidae) and their relation to growth and maturity. Cah. Biol. Mar. 51, 275–287 (2010).
Google Scholar 
Ceriola, L., Ungaro, N. & Toteda, F. Some information on the biology of Illex coindetii Verany, 1839 (Cephalopoda, Ommastrephidae) in the south-western Adriatic Sea (central Mediterranean). Fish. Res. 82, 41–49 (2006).Article 

Google Scholar 
Arvanitidis, C. et al. A comparison of the fishery biology of three Illex coindetii Verany, 1839 (Cephalopoda: Ommastrephidae) populations from the European Atlantic and Mediterranean Waters. Bull. Mar. Sci. 71, 129–146 (2002).
Google Scholar 
Quetglas, A., Alemany, F., Carbonell, A., Merella, P. & Sanchez, P. Some aspects of the biology of Todarodes sagittatus (Cephalopoda: Ommastrephidae) from the Balearic Sea (western Mediterranean). Sci. Mar. 62, 73–82 (1998).Article 

Google Scholar 
Krstulovic Sifner, S. K. & Vrgoc, N. Population structure, maturation and reproduction of the European squid, Loligo vulgaris, in the central Adriatic Sea. Fish. Res. 69, 239–249 (2004).Article 

Google Scholar 
Moreno, A. et al. Biological variation of Loligo vulgaris (Cephalopoda: Loliginidae) in the eastern Atlantic and Mediterranean. Bull. Mar. Sci. 71(1), 515–534 (2002).
Google Scholar 
Guerra, A. & Manriquez, M. Parametros biometricos de Octopus vulgaris. Invest. Pesq. 44, 177–198 (1980).
Google Scholar 
Quetglas, A., Alemany, F., Carbonell, A., Merella, P. & Sanchez, P. Biology and fishery of Octopus vulgaris Cuvier, 1797, caught by trawlers in Mallorca (Balearic Sea, western Mediterranean). Fish. Res. 36, 237–249 (1998).Article 

Google Scholar 
Sanchez, P., & Obarti, R. 1993. The biology and fishery of Octopus vulgaris caught with clay pots on the Spanish Mediterranean coast. In: Jereb, P., Allcock, A. L., Lefkaditou, E., Piatkowski, U., Hastie, L. C., Pierce, G. J. (Eds.) 2015. Cephalopod biology and fisheries in Europe: II. Species Accounts. ICES Cooperative Research Report No. 325, p 360.Gonzalez, M., Barcala, E., Perez-Gil, J. L., Carrasco, M. N. & Garcia-Martinez, M. C. Fisheries and reproductive biology of Octopus vulgaris (Mollusca: Cephalopoda) in the Gulf of Alicante (Northwestern Mediterranean). Medit. Mar. Sci. 12, 369–389 (2011).Article 

Google Scholar 
Jabeur, C., Nouira, T., Khoufi, W., Mosbahi, D. S. & Ezzedine-Najai, S. Age and growth of Octopus vulgaris Cuvier, 1797 along the east coast of Tunisia. J. Shellf. Res. 31, 119–124 (2012).Article 

Google Scholar 
Quetglas, A., Ordines, F., Gonzalez, M. & Franco, I. Life history of the bathyal octopus Pteroctopus tetracirrhus (Mollusca, Cephalopoda) in the Mediterranean Sea. Deep Sea Res. Part I 56, 1379–1390 (2009).Article 

Google Scholar 
Quetglas, A., Gonzalez, M. & Franco, I. Biology of the upper-slope cephalopod Octopus salutii from the western Mediterranean Sea. Mar. Biol. 146, 1131–1138 (2005).Article 

Google Scholar 
Moriyasu, M. Etude biometrique de la croissance d’E. cirrhosa [LAM. 1798 (Cephalopoda, Octopoda)] du Golfe du Lion. Oceanol. Acta 6, 35–41 (1983).
Google Scholar 
Massi, D. Effetti del congelamento sull’accuratezza delle misure in Eledone cirrhosa (Lamarck, 1798). Biol. Mar. Suppl. al Notiziario S.I.B.M. 1, 379–380 (1993).
Google Scholar 
Agnesi, S., Belluscio, A. & Ardizzone, G. D. Biologia e dinamica di populazione di Eledone cirrhosa (Cephalopoda: Octopoda) nel Tirreno Centrale. Biol. Mar. Mediterr. 5, 336–348 (1998).
Google Scholar 
Giordano, D. et al. Population dynamics and distribution of Eledone cirrhosa (Lamarck, 1798) in the Southern Tyrrhenian Sea (Central Mediterranean). Cah. Biol. Mar. 51, 213–227 (2010).
Google Scholar 
Krstulovic Sifner, S. K. & Vrgoc, N. Reproductive cycle and sexual maturation of the musky octopus Eledone moschata (Cephalopoda: Octopodidae) in the northern and central Adriatic Sea. Sci. Mar. 73, 439–447 (2009).Article 

Google Scholar 
Ikica, Z., Krstulovic Sifner, S. & Joksimovic, A. Some preliminary data on biological aspects of the musky octopus, Eledone moschata (Lamarck, 1798) (Cephalopoda: Octopodidae) in Montenegrin waters. Stud. Mar. 25, 21–36 (2011).
Google Scholar 
Akyol, O., Şen, H. & Kinacigil, H. T. Reproductive biology of Eledone moschata (Cephalopoda: Octopodidae) in the Aegean Sea (Izmir Bay, Turkey). J. Mar. Biol. Assoc. UK 87, 967–970 (2007).Article 

Google Scholar 
Quetglas, A., Gonzalez, M., Carbonell, A. & Sanchez, P. Biology of the deep-sea octopus Bathypolypus sponsalis (Cephalopoda: Octopodidae) from the western Mediterranean Sea. Mar. Biol. 138, 785–792 (2001).Article 

Google Scholar  More

For a schematic overview of the experimental design, see Fig. 2.Experimental evolutionProtocolThe stock population (Stock population below) was allowed to expand for one generation and from this we established eight replicate experimental evolution populations, four selected for fighter morphs (F-lines) and four selected for scrambler morphs (S-lines). Each population was founded by 1,000 recently eclosed adults (500 random females and 500 random males of the desired morph). The classification of the morphs was based on visual inspection using a stereoscopic microscope and was unambiguous due to the discontinuous distribution of the phenotypes (Classifying male morphs below). Adults were allowed to interact freely for 6 days, all surviving adults (with previously laid eggs discarded) were transferred to a new container for 24 h of egg laying, after which adults were removed. The resulting offspring were allowed to mature over 13 days and 1,000 individuals from the newly eclosed adults selected for founding the following generation, again 500 random females and 500 random males of the desired morph, with this protocol repeated every generation (Extended Data Fig. 1). The isolation of nymphs to use virgins was unfeasible with our experimental design and population sizes. However, the period of 6 days after selecting the founders of the next generation and collecting eggs for the next generation was probably enough to displace most sperm stored by females mated with any unselected males due to the high number of remating that will be occurring over this duration (females on average remate after 80 min, ref. 88) and last male sperm precedence89. The timing of generation was chosen to reflect maturation rates from our stock population to avoid indirect selection on this trait. Moreover, a previous study90 showed there was no difference between male morphs in maturation rates and that over similar lengths of time to the protocol here the fertility of both morphs remains similar. Therefore, our protocol was not likely to impose strong differential selection on morph life histories.Tracking morph proportionWe assayed the proportion of male morph in each population every 6–7 generations, by isolating 200 larvae (ten per vial) from the container, allowing maturation within vials and recording the morph of all males that eclosed (mean n = 86 per population, per generation, range 71–109). Our selection protocol was highly effective in driving an increase in the frequency of the desired male morph to >90% after 20 generations in both treatments, with this effect considerably faster within F-lines indicted by a significant two-way interaction between proportion of the desired morph and generation (χ2 = 39.9, d.f. = 6, P 90% is probably a consequence of a longer interaction period (3 versus 6 days) in which the stored sperm of males before selection was able to be displaced and/or because selection was acting more efficiently in our larger populations. The difference between rate of changes in morph proportion between F- and S-lines in the current study, and also found by Plesnar-Bielak et al.39, may be associated with the genetic architecture of morph expression. Alternatively, selection could be less effective in scrambler lines if they are less efficient than fighters in displacing sperm of previous females’ partners, but this is unlikely as R. robini male morphs have previously been demonstrated to not differ in their sperm competiveness89.Stock populationWe established a stock population by mixing three laboratory populations that were collected from three sites in Poland (Krakόw, collected in 1998 and 2008, Kwiejce, collected in 2017 and Mosina, collected in 2017; Extended Data Fig. 1), where the line derived for material used in creating the reference genome (below) was also established from the same collections at Mosina in 2017. All populations were maintained in cultures with several hundred individuals per generation before mixing and establishment of the stock population. The mixing of distinct populations increased the genetic variance in the stock population, which otherwise would probably have been limited due to founder events and the limited population size of each of the contributing populations73, thus decreasing our power to detect the effects of SSTs on genetic variation. The newly mixed stock population was maintained with several hundred individuals per generation for roughly 12 generations before the onset of this experiment. This time period is probably enough to break linkage disequilibria that could have arisen due to mixing (for unlinked loci, linkage disequilibrium should decay by half each generation91).One generation before establishing experimental evolution populations the proportion of male morphs was determined from 176 random males, indicating a roughly equal morph ratio (95 fighters, 81 scramblers) of the stock population (Extended Data Fig. 3).General housing and husbandryThe stock population and experimental evolution populations were maintained in plastic containers (approximate length, 9.5 cm; width, 7 cm; height, 4.5 cm), filled with roughly 1 cm of plaster-of-Paris. The same containers were used when sampling mites for sequencing for the reference genome or resequencing from experimental evolution populations, but either replaced the plaster-of-Paris with 5% agarose gel or added a thin layer of 5% agarose gel above the plaster-of-Paris, respectively. The agarose gel was used to reduce the number of contaminates within our samples and on the basis of preliminary extractions that indicated that small pieces of plaster-of-Paris may reduce the quality of DNA during extractions. Individuals, pairs and small groups of ten mites were housed in glass vials (approximate height, 2 cm; diameter, 0.8 cm) and large groups of 60 or 150 mites in plastic containers (approximate height, 1.5 cm; diameter, 2 cm diameter or height, 1.5 cm; diameter, 3.5 cm diameter, respectively) all with an approximate 1 cm base of plaster-of-Paris. All plaster-of-Paris bases were completely soaked in water before mites were transferred into them. All mites were reared at a constant 23 °C, at high humidity ( >90%) and were provided an excess of powdered yeast ad libitum.Classifying male morphsTo illustrate the discontinuous distribution of the weapon and to demonstrate that this classification based on visual inspection is non-subjective, we performed phenotypic measurements from male mites from a population collected near Krakόw, Poland, that had previously been fixed onto microscope slides for a separate study66. The measurements taken were idiosoma (body without mouthparts) length and width of third proximal segment of the third right leg (genu). Measurements were preformed using Lecia DM5500B microscope and Lecia Application Suite v.4.6.1. We then performed an analysis to, first, determine whether the allometric relationship between idiosoma length and width of third pair of legs is best described as discontinuous and, second, to verify that classification by simple visual inspection matches the same classification from allometric analysis. One researcher performed all the measurements and classified each male as a fighter (n = 50) or scrambler (n = 50), a separate researcher was then given the measurements but not the classification of the male morph.Broadly, guidelines for the analysis of non-linear allometries92 were followed. The log–log scatterplots of idiosoma length against leg width were visualized, which showed there was clear evidence for non-linear scaling relationships. Next histograms of idiosoma length, leg width and relative leg width (leg width/idiosoma length) were visualized (Extended Data Fig. 2a–c). Where a normal distribution of idiosoma length, and a binomial distribution in leg width and relative leg width are further indications of a discontinuous relationship. On the basis of the lowest point between the two peaks of the density plot of relative leg width (Extended Data Fig. 2c) males were classified as scramblers (relative leg width 0.125). Replotting the log–log scatterplot of idiosoma length and leg width, and using the classification of morph described above clearly demonstrates the discontinuous allometric relationship of idiosoma length and leg width in R. robini (Extended Data Fig. 2d). Moreover, on the basis of the Akaike information criterion (AIC), the discontinuous model where males were assigned a morph (AIC = 646.5) clearly has a substantially better fit than a simple linear and quadratic models (AIC = 918.5 and 920.2, respectively). Further models were omitted from comparison (for example, breakpoint or sigmoidal) due to the clear discontinuous allometry observed. Finally, all 100 males were assigned the same morph by visual inspection and blind allometric analysis, demonstrating that the former is effective and accurate in classifying male morph.Phenotypic assaysFecundity assays were performed using experimental evolution females at F20 and F32. Eggs laid by females between days 4–8 were counted, encompassing the window of time of most evolutionary relevance for female fitness during maintenance of selection lines (that is, egg laying period in selection lines was between days 6–7) and also likely to capture variation in lifetime fecundity that remains largely consistent throughout the first 3 weeks of life93. Nymphs were individually isolated to gain virgin females, which on maturation females from each experimental evolution population (n = 30) were paired with a male from the stock population (15 with fighters and 15 with scramblers). Pairs were transferred to a new vial on day 4, with the pair being removed from the second vial after a further 4 days and all eggs in the second vial counted. If the male had died in the first vial, they were replaced with a stock male of the same morph. Any female deaths in the first or second vials were recorded.Longevity assays were also performed at F20 and F32. At F20, females used in fecundity assays, including the stock male they were paired with (replaced if dead), were transferred to a new vial at day 8. After this point, vials were then checked every 2 days for female deaths and pairs were moved to new vials every 4 days. Males were replaced with stock males of the same morph if found dead. Similarly, at F20, on maturation males from experimental evolution populations (n = 30) were paired with stock females, vials were checked every 2 days and changed every 4 days, with females being replaced if dead. At F32, only female longevity was determined and was performed in groups; 30 experimental evolution females and 30 stock males (15 of each morph) were placed in plastic containers, two per experimental population. This logistically easier estimate of longevity was done due to local restrictions during the SARS-CoV-2 pandemic and the imposed limitations on people working closely together. Groups were checked for dead females every other day and all remaining live mites transferred to a new container every 4 days. When mites were transferred to a new container the sex and morph ratio were balanced to that of the remaining females, by either removing or adding males of the desired morph from the stock population.To determine whether the survival of mites differed between F- and S-lines when competition between males was allowed, at F45 we created small colonies from each population and survival of males and females recorded over 6 days, the same period as used between selecting founders of the next generation and subsequent egg laying period. Colonies were at a 50:50 sex ratio, established with 150 newly eclosed mites placed into small plastic containers. This was approximately the same density after selection of the next generations founders during the maintenance of experimental evolution populations (150 mites in roughly 9.5 cm2 = 16 mites per 1 cm2; 1,000 mites in roughly 67 cm2 = 15 mites per 1 cm2). After 3 days, all colonies were checked and any dead mites identified by sex. After another 3 days, again dead mites were recorded and all surviving mites sexed and counted.Additionally, at F45 we performed further fecundity assays to obtain estimates of inbreeding depression within experimental evolution populations. To establish family groups, larvae were isolated and on maturation F0 males and females (n = 16) from within the same experimental evolution population were paired together. Pairs were allowed to produce eggs for 48 h, after which adults were removed from vials. After hatching from each pair, 12 F1 larvae were isolated into new vials. On their maturation, these F1 mites were either paired with a full sibling, that is, from the same family, or with an individual from a different family but from the same experimental evolution population. When possible, we made two inbred and two outbred pairs with same family lines used. Again, pairs were allowed to produce eggs for 48 h before their removal for the vial. After a further 5 days, vials were checked for larvae, if larvae were present in the first vial six were individually isolated and the second vial discarded, if no larvae were present in the first vial the second vial was checked for larvae and, if present, they were isolated. This protocol therefore produced inbred and outbred individuals from within the same experimental evolution population. Which, as above, on maturation F2 inbred and outbred females were paired with stock males (fighter males only) and number of eggs laid between days 4 and 8 counted. Only a single female from each unique inbred or outbred family was used. Either due to pairs failing to produce offspring or there being no F2 females, samples sizes were not exactly equal. In total, 59 outbred and 55 inbred females from F-lines, and 56 outbred and 54 inbred females from S-lines were paired with stock males.Phenotypic assay statistical analysesAll phenotypic analysis was conducted using R statistical software94 (v.3.5.2) and data were visualized using ggplot2 (ref. 95).Analysis of male morph proportion was performed using a generalized linear mixed model with binomial error structure, fitted using lme4 (ref. 96). Where the proportion of desired morph was compared in model with morph selection and generation (as a factor) including their two-way interaction as explanatory variables, and population included as a random effect.All fecundity data were analysed using generalized linear mixed models with Poisson error structures, fitted using lme4. Due to the differences in stock population males used between F45 and earlier generations, and slightly different rearing conditions between females in the fecundity assays from generations F20 and F32, they were analysed separately from data collected in F45. However, we noted that the fecundity of females in Fig. 5a was comparable to the outbred females in Fig. 5b. Explanatory variables fitted to fecundity data from F20 and F32 were, morph selection treatment, generation, including their two-way interaction term, and stock male morph. The explanatory variables fitted to fecundity data from inbreeding depression data were, morph selection treatment and status of female (that is, inbred or outbred), including their two-way interaction term. In both analyses, we included population as a random effect and an observation level random effect to account for overdispersion, we omitted fitting random slopes due to issues with increasing the complexity of random effects close to reaching a singular fit. Females that died before the end of the fecundity assay and those that laid zero eggs were removed from analysis. This excluded five females from F20 (three F-line and two S-line), 20 from F32 (13 F-line and seven S-line) and 16 from F45 (three inbred and three outbred F-line, and nine inbred and one outbred S-line).Longevities of females at F20 and F32, and males at F20, were analysed separately using mixed effects Cox models, fitted using coxme97. In all analyses, we included a random effect of population, with morph selection treatment as an explanatory variable and extra variable of male morph included in female longevity analysis at F20. Survival of mites over 6 days at F45 was analysed using a GLM with counts of dead and surviving mites fitted with a quasibinomial error structure, the model included morph selection treatment and sex, including their interaction term, as explanatory variables. If individuals were lost due to handling error (that is, killed or escaped) they were right-censored during analysis.Genome assemblySample originA line of R. robini originated from a wild-collected population from the Mosina region (Wielkopolska, Poland). In October 2017, onions were collected from the field and approximately 200 individuals of R. robini were identified under dissecting microscope. The line used for DNA isolation in the genome sequencing project was developed from full sib × sib mating for 14 generations (to maximize homozygosity) following and continuing the protocol described in ref. 67.DNA extractionFor DNA extraction we used only mite eggs, that were laid by 500 females, collected in a container (see above for a description) Females were kept in this container for 3 days. After that time, they were removed, and eggs were filtered using fine sieves and washed for 1 min in 0.3% sodium hypochlorite solution and in Milli-Q water for 2 × 2 min to remove any potential foreign DNA contamination. These eggs were collected in 1.5 ml Eppendorf tube and after short centrifugation, the remains of the water (supernatant) removed with a pipette. The sample was immediately transferred to ice and prepared for DNA extraction. DNA was extracted using Bionano Prep Animal Tissue DNA Kit for HMW DNA isolation according to the manufacturer’s instructions. Briefly, eggs were smashed with a sterile pellet pestle on ice in 500 μl homogenization buffer; the sample was fixed with 500 μl cold ethanol and incubated 60 min on ice, after that time the sample was centrifuged at 1500g for 5 min at 4 °C and the supernatant was discarded. Next, after resuspension in a homogenization buffer pellet, this was cast in four agarose plugs as described in the original protocol. Agarose plugs were incubated with Proteinase K and Lysis buffer solution for 2 h with intermittent mixing. After that time, the digestion solution was replaced with a freshly made one and incubated overnight with intermittent mixing. According to the original protocol, after RNase A digestion and plug washing, DNA was recovered by incubation of the plugs in TE buffer, followed by plug melting and addition of agarase. Recovered DNA was dialysed and homogenized on a membrane for 45 min at room temperature and transferred to a clean tube with a wide bore tip.SequencingSequencing was done using Oxford Nanopore Technologies (ONT, MinION). Isolated DNA purified using AMPure XP beads and resuspended in H2O before library preparation. Two separate libraries were prepared using ligation sequencing kit, SQK-LSK109 and Rapid Sequencing Kit SQK-RAD004, respectively, according to the manufacturer’s protocols and were sequenced on a FLO-MIN106 R9.4.1 SpotON flow cell on a MinION Mk 1B sequencer (ONT). The total yields from sequencing were 484,700 reads (2,417,068,187 nt) with a read-N50 of 10,044 nt (ranging from 216,403 to 100). Base calling of the raw reads was done using Guppy (v.3.3) resulting in a total sum of the reads 7,979,616,172, equivalent to 26× coverage aiming for a genome of 300 megabases (Mb). The reads N50/N90 were estimated at 7,958/1,719.Assembling reference genomeReads aligning with the Mitochondrion genome were identified using BLASTN and filtered from the raw reads before assembling the genome. The remaining ONT reads were assembled using the Flye software (v.2.6), with –min-overlap 3,000 to increase stringency at the initial overlay step, and default parameters including five rounds of polishing through consensus, contigs were additionally polished two times with Medaka (v.0.11.2). Illumina paired-end RNA dataset is assembled using CLC Assembler (CLC Assembly Cell). Both RNA assemblies and paired-end 10X genomic dataset (unpublished data) were mapped onto the contigs using minimap2 (v.2.16) and BWA mapper (v.0.7.17), respectively, and the assembly was further polished using PILON (v.1.20) to error correct potential low-quality regions. The resulting assembly yielded a genome of 307 Mb, assembled into 1,533 contigs ranging from 10,840,357 to 100 basepairs (bp) and an assembly-N50 of 1.670 Mb. Moreover, the BUSCO completeness analysis using the Arachnida (odb10) reference set confirmed our assembly represents the complete genome C:94.8%(S:89.1%,D:5.7%),F:0.9%,M:4.3%,n:2934 (=arachnida_odb10), only missing 126 genes from the whole reference set. Knowing that BUSCO only gives a rough estimation, we remain confident that this assembly represents well the bulb mite genome.Flow-cytometryWhole individual R. robini were homogenized in 500 μl of ice-cold LB01 detergent buffer along with the head of a male Drosophila melanogaster (1 C = 0.18 pg) as an internal standard. The homogenized tissue was filtered through a 30-μm nylon filter. Then 12 μl of propidium iodide with 2 μl of RNase was added, and stained for 1 h on ice in the dark. All samples were run on an FC500 flow cytometer (Beckman-Coulter) using a 488-nm blue laser, providing output as single-parameter histograms showing relative fluorescence between the standard nuclei and the R. robini nuclei. Six replicate samples were run to account for variation in fluorescence outputs. The genome size of R. robini was estimated at 0.30 pg, or about 293 Mb, and consistent with estimation of size from the genome assembly described above.Mitochondrial genomeONT reads aligned with R. robini mitochondrion genome were de novo assembled with Flye (v.2.6) assembler and polished with Racon. Mitochondrion genome is assembled in one single contig with a size of 15,335 bases.Gene predictionOn the polished final genome, protein coding genes have been predicted. For this, AUGUSTUS was used including hints coming from R. robini RNA-sequencing (RNA-seq) (samples SRR3934324, SRR3934325, SRR3934326, SRR3934327, SRR3934328, SRR3934329, SRR3934330, SRR3934331, SRR3934332, SRR3934333, SRR3934335, SRR3934337, SRR3934338 and SRR3934339 from the PRJNA330592 BioProject deposited at the National Center for Biotechnology Information (NCBI) Short Read Archive) and proteins coming from highly curated Tetranychus urticae (v.2020-03-20) as well as proteins from the previous version of the unpublished, Illumina-sequenced R. robini genome (https://public-docs.crg.eu/rguigo/Data/fcamara/bulbmite.v4a/). The PE RNA-seq reads were mapped on the genome using HISAT2 (-k 1 —no-unal) and further processed with Regtools to extract junction hints and filtered for junctions with a minimum coverage of 10. All the RNA-seq reads were also assembled with CLC Assembly Cell (v.5.2.0) software, setting the word size for the Bruijn graph at 50 and maximum bubble size at 31. The reads were assembled into 689,563 contigs (ranging from 10,675 to 180 bp), which were later mapped on the genome with GenomeTheader to generate complementary DNA hints. Protein hints were generated by using with Exonerate (v.2.2) with Protein2Genome model. To reduce the amount of overprediction due to repeated elements (transposable elements, simple sequence repeats) we de novo predicted high abundant repeats using RepeatModeler. The accompanying parameter file for extrinsic data for AUGUSTUS was adapted to include these hints as well as the softmasking of the genomic sequence. The resulting gene predictions from AUGUSTUS were further curated with EvidenceModeler using the same extrinsic data. The BUSCO analysis confirmed that our gene prediction indeed captured the expected genes well (C:94.6%(S:86.3%,D:8.3%),F:0.4%,M:5.0%,n:2934 (=arachnida_odb10)). The final predicted gene set was subsequently processed to be uploaded into ORCAE (https://bioinformatics.psb.ugent.be/orcae/overview/Rhrob)98.ResequencingGenomic sampling and mappingFor genomic analyses we sampled material from each of the morph selection lines (n = 8) at F1, F12 and F29. Following the experimental evolution protocols, after the first 24 h of egg laying all adults were transferred to a new container (described above) for a second 24 h to lay eggs and from these second dishes genomic material was sampled. On maturation, adults were transferred to and kept for 3 days in containers. Adults were then randomly selected and placed into Digestion Solution for MagJET gDNA Kit (F1&12) or ATL buffer (F29) before freezing at −20 °C. From each population two samples were collected consisting of 100 individuals (1 × 100 females and 1 × 100 males of random morph), the two samples separated by sex were used as technical replicates. The tissue from the 100 individuals within each sample was homogenized and DNA was extracted by Proteinase K digestion (24 h) followed by standard procedures using MagJET Genomic DNA Kit (ThermoScientific, F1&12) or DNeasy Blood and Tissue (Qiagen, F29). DNA concentration was controlled with the Qubit double-stranded DNA HS Assay Kit and DNA quality was assessed on agarose gels. The library preparation was performed using NEBNext Ultra II FS DNA Library Prep kit for Illumina.Whole-genome resequencing was carried out by National Genomics Infrastructure (Uppsala, Sweden) using the Illumina Nova-Seq 6000 platform with S4 flow cell to produce 2 × 150 bp reads (average 160.7 × 106: range 130.7 × 106 − 189.9 × 106). Adaptors were trimmed from reads using Trimmomatic99 software (v.0.39) and unpaired reads discarded. Fastq files were mapped to the assembled genome with bwa mem100 (v.0.7.17-r1188) using default settings. Sam files were converted to bam files, sorted, duplicates marked and ambiguously mapped reads removed using samtools101 (v.1.9). On average, 90% (range, 86–93%) of the reads from each sample were mapped successfully, of which an average of 17% (range, 15–19%) were marked as duplicates. This left us with an average of 117.7 × 106 pair end reads per sample, ranging between 99.6 × 106 and 145.9 × 106 (Supplementary Table 1).Genomic analysisFile preparation and filteringPreparation of files used in genomic analysis was done as follows: bam files were converted to a pile-up file using samtools, following which indels and surrounding windows (5 bp either side) were filtered, using identify-genomic-indel-regions.pl and filter-pileup-by-gtf.pl in PoPoolation102 (v.1.2.2) to avoid false SNPs, with the resulting filtered pile-up file converted to a sync file using mpileup2sync.pl in PoPoolation2 (ref. 103) (v.1.201). Using custom python scripts, the distribution of coverage from each sample (single sex) was determined by recording the coverage of positions every 10 kb across the genome from the sync file to give information on expected coverage (Supplementary Fig. 1). On the basis of this, we filtered the sync and pile-up files to contain only regions within a range of informative coverage, where the mean coverage of all samples at every position was between 50% of the expected coverage and 200% of the expected coverage (56×, range 23−112×). The pile-up and sync files containing individual male and female samples (48 in total) were then merged by sex to give files containing allele frequencies from 24 samples (eight populations across three generations), each consisting of allele frequencies of 200 individuals (100 males and 100 females, above) and used in all subsequent analysis (unless stated otherwise). Similarly, we drew coverage of a position every 10 kb from each sample in the sex-merged sync file to determine a distribution from which we decided to subsample to (Supplementary Fig. 1). We putatively identified X-linked contigs (below) and excluded them autosomal analysis. A similar, but, separate analysis on genes and SNPs from X-linked contigs was performed by using different parameters (below).Estimating nucleotide diversityUsing PoPoolation we determined various estimates of genetic diversity per sample (that is, 24 sex-merged samples). The pile-up file from each sample was subsampled using subsample-pileup.pl to a coverage of 63× (max coverage, 252×) to standardize estimations of genetic diversity across the genome, between populations and across generations. First, nucleotide diversity (Tajima’s Pi, π) and number of segregating sites (Watterson’s theta, ϴ) were estimated within genes. We performed analysis of exons using Syn-nonsyn-at-position.pl, in which genetic diversity of synonymous and non-synonymous positions were determined. Further analysis of overall genetic diversity within exons and introns were performed using Variance-at-position.pl, Tajima’s D (D) also estimated in the former. We used a minimum count of three (equal to a minor allele frequency of roughly 5%) for a SNP to be called, and a phred score >30 and a pool size of 400. Further analysis using 10 kb sliding windows (step size 10 kb) across the genome were performed using Variance-sliding.pl, and also included estimation of D. Estimates of D require the minimum count to be 2, but otherwise all the same parameters were used.We filtered genes to be included in our analysis (and all subsequent analysis) on the basis of a number of criteria. On the basis of extensive RNA-seq data from both males and females (Plesnar-Bielak, unpublished data with NCBI accession number PRJNA796800), we only included genes in our analyses that were expressed at a mean level of fragments per kilobase of transcript per million mapped reads >1 across 72 samples originated from both sexes and both morphs rearing in three different temperatures (18, 23 and 28 °C). A further filtering step was performed to remove genes with inconsistent mapping between samples, only genes with >60% exons mapped to (calculated from positions used to calculate parameters in the Syn-nonsyn-at-position.pl π outputs), with 63−252× coverage, in all 24 samples were included in the analysis. The final dataset contained 13,389 autosomal genes and subsequently used to filter other datasets to retain this set of genes only (see Supplementary Table 8 for a list of genes). Similarly, windows were discarded from outputs if 60% of genes being mapped to in all 24 samples) and reducing the final X-linked dataset to contain fewer than 200 genes. We therefore opted to reduce the target coverage further to 40×, in an attempt to retain more genes. This slight reduction of target coverage increased the number of genes in the final dataset substantially to 587 genes. We therefore opted to use a minimum coverage of 40× in all analysis of X-linked SNPs, genes and windows.Diverging SNPsTo determine divergent SNPs between F- and S-lines, we extracted the allele frequencies of all samples from the sex combined sync file. Samples from F29 were then used to filter the entire dataset to only contain SNPs on the basis of a number of criteria. First, positions within all samples were required to have a coverage >63× and 5% (that is, the average of all samples but not necessarily above >5% in all samples). Thus, our dataset contained only positions with the target coverage in all F29 samples and in which polymorphisms were unlikely to be a consequence of sequencing errors. After this filtering we were left with roughly 6 million SNPs used in further analysis. We performed a GLM, at each position by comparing the count of the major allele against counts of minor alleles at F29, to determine consistent allele frequency changes between treatments70. If any population had minor or major allele count of 0, +1 was added to minor and major alleles from all samples. To correct for multiple testing, we converted P values to q values using the qvalues R package (v.2.14.1)104 and applied a FDR with a q 900,000), GLMs were performed (identical to above) on the simulated major and minor allele counts. Using a FDR with a q  More

Martin WF, Bryant DA, Beatty JT. A physiological perspective on the origin and evolution of photosynthesis. FEMS Microbiol Rev. 2018;42:205–31.CAS 
PubMed 
Article 

Google Scholar 
Partensky F, Blanchot J, Vaulot D. Differential distribution and ecology of Prochlorococcus and Synechococcus in oceanic waters: a review. Bull Oceanogr Monaco, no Spec. 1999;19:457–76.
Google Scholar 
Zwirglmaier K, Jardillier L, Ostrowski M, Mazard S, Garczarek L, Vaulot D, et al. Global phylogeography of marine Synechococcus and Prochlorococcus reveals a distinct partitioning of lineages among oceanic biomes. Environ Microbiol. 2008;10:147–61.PubMed 

Google Scholar 
Callieri C. Picophytoplankton in freshwater ecosystems: the importance of small-sized phototrophs. Freshw Rev. 2008;1:1–28.Article 

Google Scholar 
Stal LJ. Physiological ecology of cyanobacteria in microbial mats and other communities. New Phytol. 1995;131:1–32.CAS 
PubMed 
Article 

Google Scholar 
Rikkinen J. Cyanobacteria in terrestrial symbiotic systems. In: Hallenbeck PC editor. Modern topics in the phototrophic prokaryotes. Switzerland: Springer; 2017. p. 243–94.Badger MR, Price GD, Long BM, Woodger FJ. The environmental plasticity and ecological genomics of the cyanobacterial CO2 concentrating mechanism. J Exp Bot. 2006;57:249–65.CAS 
PubMed 
Article 

Google Scholar 
Rae BD, Long BM, Badger MR, Price GD. Functions, compositions, and evolution of the two types of carboxysomes: polyhedral microcompartments that facilitate CO2 fixation in cyanobacteria and some proteobacteria. Microbiol Mol Biol Rev. 2013;77:357–79.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Bar-On YM, Phillips R, Milo R. The biomass distribution on Earth. Proc Natl Acad Sci USA. 2018;115:6506–11.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Buitenhuis ET, Li WKW, Vaulot D, Lomas MW, Landry MR, Partensky F, et al. Picophytoplankton biomass distribution in the global ocean. Earth Syst Sci Data. 2012;4:37–46.Article 

Google Scholar 
Flombaum P, Gallegos JL, Gordillo RA, Rincón J, Zabala LL, Jiao N, et al. Present and future global distributions of the marine Cyanobacteria Prochlorococcus and Synechococcus. Proc Natl Acad Sci USA. 2013;110:9824–9.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Field CB, Behrenfeld MJ, Randerson JT, Falkowski P. Primary production of the biosphere: integrating terrestrial and oceanic components. Science. 1998;281:237–40.CAS 
PubMed 
Article 

Google Scholar 
Garcia-Pichel F, Belnap J, Neuer S, Schanz F. Estimates of global cyanobacterial biomass and its distribution. Arch Hydrobiol Suppl Algol Stud. 2003;109:213.
Google Scholar 
Scanlan DJ, Ostrowski M, Mazard S, Dufresne A, Garczarek L, Hess WR, et al. Ecological genomics of marine picocyanobacteria. Microbiol Mol Biol Rev. 2009;73:249–99.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Doré H, Farrant GK, Guyet U, Haguait J, Humily F, Ratin M, et al. Evolutionary mechanisms of long-term genome diversification associated with niche partitioning in marine picocyanobacteria. Front Microbiol. 2020;11:2129.Article 

Google Scholar 
Dufresne A, Ostrowski M, Scanlan DJ, Garczarek L, Mazard S, Palenik BP, et al. Unraveling the genomic mosaic of a ubiquitous genus of marine cyanobacteria. Genome Biol. 2008;9:R90.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Badger MR, Hanson D, Price GD. Evolution and diversity of CO2 concentrating mechanisms in cyanobacteria. Funct Plant Biol. 2002;29:161–73.CAS 
PubMed 
Article 

Google Scholar 
Whitehead L, Long BM, Price GD, Badger MR. Comparing the in vivo function of α-carboxysomes and β-carboxysomes in two model cyanobacteria. Plant Physiol. 2014;165:398–411.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Castenholz RW, Wilmotte A, Herdman M, Rippka R, Waterbury JB, Iteman I, et al. Phylum BX. cyanobacteria. In: Boone DR, Castenholz RW, Garrity GM editors. Bergey’s manual of systematic bacteriology. New York, NY: Springer; 2001. p. 473–599.Cabello‐Yeves PJ, Picazo A, Camacho A, Callieri C, Rosselli R, Roda-Garcia JJ, et al. Ecological and genomic features of two widespread freshwater picocyanobacteria. Environ Microbiol. 2018;20:3757–71.PubMed 
Article 
CAS 

Google Scholar 
Di Cesare A, Cabello-Yeves PJ, Chrismas NAM, Sánchez-Baracaldo P, Salcher MM, Callieri C, et al. Genome analysis of the freshwater planktonic Vulcanococcus limneticus sp. nov. reveals horizontal transfer of nitrogenase operon and alternative pathways of nitrogen utilization. BMC Genomics. 2018;19:259.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Sánchez-Baracaldo P, Bianchini G, Di Cesare A, Callieri C, Chrismas NAM. Insights into the evolution of picocyanobacteria and phycoerythrin genes (mpeBA and cpeBA). Front Microbiol. 2019;10:a45.Article 

Google Scholar 
Callieri C, Mandolini E, Bertoni R, Lauceri R, Picazo A, Camacho A, et al. Atlas of picocyanobacteria monoclonal strains from the collection of CNR-IRSA, Italy. J Limnol. 2021;80:2002.Article 

Google Scholar 
Herdman M, Castenholz RW, Iteman I, Waterbury JB, Rippka R. Subsection I (Formerly Chroococcales Wettstein 1924, emend. Rippka, Deruelles, Waterbury, Herdman and Stanier 1979). In: Boone DR, Castenholz RW, Garrity GM, editors. Bergey’s manual of systematic bacteriology, Vol 1, 2nd ed., The archaea and the deeply branching and phototrophic bacteria. New York: Springer; 2001. p. 493–514.Cabello-Yeves PJ, Haro-Moreno JM, Martin-Cuadrado A, Ghai R, Picazo A, Camacho A, et al. Novel Synechococcus genomes reconstructed from freshwater reservoirs. Front Microbiol. 2017;8:1151.PubMed 
PubMed Central 
Article 

Google Scholar 
Badger MR, Price GD. CO2 concentrating mechanisms in cyanobacteria: molecular components, their diversity and evolution. J Exp Bot. 2003;54:609–22.CAS 
PubMed 
Article 

Google Scholar 
Wheatley NM, Sundberg CD, Gidaniyan SD, Cascio D, Yeates TO. Structure and identification of a pterin dehydratase-like protein as a ribulose-bisphosphate carboxylase/oxygenase (RuBisCO) assembly factor in the α-carboxysome. J Biol Chem. 2014;289:7973–81.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Huang F, Kong WW, Sun Y, Chen T, Dykes GF, Jiang Y, et al. Rubisco accumulation factor 1 (Raf1) plays essential roles in mediating Rubisco assembly and carboxysome biogenesis. Proc Natl Acad Sci USA. 2020;117:17418–28.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Kerfeld CA, Melnicki MR. Assembly, function and evolution of cyanobacterial carboxysomes. Curr Opin Plant Biol. 2016;31:66–75.CAS 
PubMed 
Article 

Google Scholar 
Kupriyanova E, Pronina N, Los D. Carbonic anhydrase—a universal enzyme of the carbon-based life. Photosynthetica. 2017;55:3–19.CAS 
Article 

Google Scholar 
DiMario RJ, Machingura MC, Waldrop GL, Moroney JV. The many types of carbonic anhydrases in photosynthetic organisms. Plant Sci. 2018;268:11–17.Heinhorst S, Cannon GC. A novel evolutionary lineage of carbonic anhydrase (epsilon class) is a component of the carboxysome shell. J Bacteriol. 2004;186:623–30.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Sawaya MR, Cannon GC, Heinhorst S, Tanaka S, Williams EB, Yeates TO, et al. The structure of beta-carbonic anhydrase from the carboxysomal shell reveals a distinct subclass with one active site for the price of two. J Biol Chem. 2006;281:7546–55.CAS 
PubMed 
Article 

Google Scholar 
Heinhorst S, Williams EB, Cai F, Murin CD, Shively JM, Cannon GC. Characterization of the carboxysomal carbonic anhydrase CsoSCA from Halothiobacillus neapolitanus. J Bacteriol. 2006;188:8087–94.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Omata T, Takahashi Y, Yamaguchi O, Nishimura T. Structure, function and regulation of the cyanobacterial high-affinity bicarbonate transporter, BCT1. Funct Plant Biol. 2002;29:151–9.CAS 
PubMed 
Article 

Google Scholar 
Omata T, Price GD, Badger MR, Okamura M, Gohta S, Ogawa T. Identification of an ATP-binding cassette transporter involved in bicarbonate uptake in the cyanobacterium Synechococcus sp. strain PCC 7942. Proc Natl Acad Sci USA. 1999;96:13571–6.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Shelden MC, Howitt SM, Price GD. Membrane topology of the cyanobacterial bicarbonate transporter, BicA, a member of the SulP (SLC26A) family. Mol Membr Biol. 2010;27:12–22.PubMed 
Article 
CAS 

Google Scholar 
Price GD, Howitt SM. The cyanobacterial bicarbonate transporter BicA: its physiological role and the implications of structural similarities with human SLC26 transporters. Biochem Cell Biol. 2011;89:178–88.CAS 
PubMed 
Article 

Google Scholar 
Price GD, Woodger FJ, Badger MR, Howitt SM, Tucker L. Identification of a SulP-type bicarbonate transporter in marine cyanobacteria. Proc Natl Acad Sci USA. 2004;101:18228–33.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Bonfil DJ, Ronen-Tarazia M, Sültemeyer D, Lieman-Hurwitza J, Schatz D, Kaplan A. A putative HCO−3 transporter in the cyanobacterium Synechococcus sp. strain PCC 7942. FEBS Lett. 1998;430:236–40.CAS 
PubMed 
Article 

Google Scholar 
Price GD, Shelden MC, Howitt SM. Membrane topology of the cyanobacterial bicarbonate transporter, SbtA, and identification of potential regulatory loops. Mol Membr Biol. 2011;28:265–75.CAS 
PubMed 
Article 

Google Scholar 
Shibata M, Katoh H, Sonoda M, Ohkawa H, Shimoyama M, Fukuzawa H. Genes essential to sodium-dependent bicarbonate transport in cyanobacteria: function and phylogenetic analysis. J Biol Chem. 2002;277:18658–64.CAS 
PubMed 
Article 

Google Scholar 
Zhang P, Battchikova N, Jansen T, Appel J, Ogawa T, Aro EM. Expression and functional roles of the two distinct NDH-1 complexes and the carbon acquisition complex NdhD3/NdhF3/CupA/Sll1735 in Synechocystis sp PCC 6803. Plant Cell. 2004;16:3326–40.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Price GD, Badger MR, Woodger FJ, Long BM. Advances in understanding the cyanobacterial CO2-concentrating-mechanism (CCM): functional components, Ci transporters, diversity, genetic regulation and prospects for engineering into plants. J Exp Bot. 2008;59:1441–61.CAS 
PubMed 
Article 

Google Scholar 
Battchikova N, Eisenhut M, Aro E-M. Cyanobacterial NDH-1 complexes: novel insights and remaining puzzles. Biochim Biophys Acta Bioenerg. 2011;1807:935–44.CAS 
Article 

Google Scholar 
Koester RP, Pignon CP, Kesler DC, Willison RS, Kang M, Shen Y. Transgenic insertion of the cyanobacterial membrane protein ictB increases grain yield in Zea mays through increased photosynthesis and carbohydrate production. PLoS ONE. 2021;16:e0246359.Johnson ZI, Zinser ER, Coe A, McNulty NP, Woodward EMS, Chisholm SW. Niche partitioning among Prochlorococcus ecotypes along ocean-scale environmental gradients. Science. 2006;311:1737–40.CAS 
PubMed 
Article 

Google Scholar 
Callieri C, Coci M, Corno G, Macek M, Modenutti B, Balseiro E. Phylogenetic diversity of nonmarine picocyanobacteria. FEMS Microbiol Ecol. 2013;85:293–301.CAS 
PubMed 
Article 

Google Scholar 
Schallenberg LA, Pearman JK, Burns CW, Wood SA. Spatial abundance and distribution of picocyanobacterial communities in two contrasting lakes revealed using environmental DNA metabarcoding. FEMS Microbiol Ecol. 2021;97:fiab075.CAS 
PubMed 
Article 

Google Scholar 
Mózes A, Présing M, Vörös L. Seasonal dynamics of picocyanobacteria and picoeukaryotes in a large shallow lake (Lake Balaton, Hungary). Int Rev Hydrobiol. 2006;91:38–50.Article 
CAS 

Google Scholar 
Vörös L, Callieri C, Balogh KV, Bertoni R. Freshwater picocyanobacteria along a trophic gradient and light quality range. Hydrobiologia. 1998;369/370:117–25.Watanabe MF, Harada K, Carmichael WW, Fujiki H. Toxic microcystis. Boca Raton, FL: CRC Press; 1995.Stockner J, Callieri C, Cronberg G. Picoplankton and other non-bloom-forming cyanobacteria in lakes. In: Whitton BA, Potts M editors. The ecology of cyanobacteria. The Netherlands: Springer; 2000. p. 195–231.Flamholz AI, Prywes N, Moran U, Davidi D, Bar-On YM, Oltrogge LM. Revisiting trade-offs between Rubisco kinetic parameters. Biochemistry. 2019;58:3365–76.CAS 
PubMed 
Article 

Google Scholar 
Filazzola A, Mahdiyan O, Shuvo A, Ewins C, Moslenko L, Sadid T. A database of chlorophyll and water chemistry in freshwater lakes. Sci Data. 2020;7:1–10.Article 

Google Scholar 
Cabello‐Yeves PJ, Zemskaya TI, Zakharenko AS, Sakirko MV, Ivanov VG, Ghai R, et al. Microbiome of the deep Lake Baikal, a unique oxic bathypelagic habitat. Limnol Oceanogr. 2019.Rodrigo MA, Miracle MR, Vicente E. The meromictic Lake La Cruz (Central Spain). Patterns of stratification. Aquat Sci. 2001;63:406–16.Article 

Google Scholar 
Camacho A, Picazo A, Miracle MR, Vicente E. Spatial distribution and temporal dynamics of picocyanobacteria in a meromictic karstic lake. Arch Hydrobiol Suppl Algol Stud. 2003;109:171–84.
Google Scholar 
Vicente E, Camacho A, Rodrigo MA. Morphometry and physico-chemistry of the crenogenic meromictic Lake El Tobar (Spain). Int Ver für Theor und Angew Limnol Verhandlungen. 1993;25:698–704.CAS 

Google Scholar 
Camacho A, Miracle MR, Vicente E. Which factors determine the abundance and distribution of picocyanobacteria in inland waters? A comparison among different types of lakes and ponds. Arch für Hydrobiol. 2003;157:321–38.Article 

Google Scholar 
Kaźmierczak J, Kempe S, Kremer B, López-García P, Moreira D, Tavera R. Hydrochemistry and microbialites of the alkaline crater lake Alchichica, Mexico. Facies. 2011;57:543–70.Article 

Google Scholar 
Ghai R, Mizuno CM, Picazo A, Camacho A, Rodriguez‐Valera F. Key roles for freshwater Actinobacteria revealed by deep metagenomic sequencing. Mol Ecol. 2014;23:6073–90.CAS 
PubMed 
Article 

Google Scholar 
Cabello-Yeves PJ, Ghai R, Mehrshad M, Picazo A, Camacho A, Rodriguez-Valera F. Reconstruction of diverse Verrucomicrobial genomes from metagenome datasets of freshwater reservoirs. Front Microbiol. 2017;8:2131.PubMed 
PubMed Central 
Article 

Google Scholar 
de Hoyos C, Negro AI, Aldasoro JJ. Cyanobacteria distribution and abundance in the Spanish water reservoirs during thermal stratification. Limnetica. 2004;23:119–32.Article 

Google Scholar 
Raven J, Caldeira K, Eldefield H, Hoegh-Guldberg O, Liss P, Riebesell U, et al. Ocean acidification due to increasing atmospheric carbon dioxide. London: The Royal Society; 2005.Mangan NM, Flamholz A, Hood RD, Milo R, Savage DF. pH determines the energetic efficiency of the cyanobacterial CO2 concentrating mechanism. Proc Natl Acad Sci USA. 2016;113:E5354–62.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Tadesse I, Green FB, Puhakka JA. Seasonal and diurnal variations of temperature, pH and dissolved oxygen in advanced integrated wastewater pond system® treating tannery effluent. Water Res. 2004;38:645–54.CAS 
PubMed 
Article 

Google Scholar 
Gao Y, Zhang Z, Liu X, Yi N, Zhang L, Song W, et al. Seasonal and diurnal dynamics of physicochemical parameters and gas production in vertical water column of a eutrophic pond. Ecol Eng. 2016;87:313–23.Article 

Google Scholar 
Schindler DW. Recent advances in the understanding and management of eutrophication. Limnol Oceanogr. 2006;51:356–63.Article 

Google Scholar 
Bosak T, Bush JWM, Flynn MR, Liang B, Ono S, Petroff AP, et al. Formation and stability of oxygen‐rich bubbles that shape photosynthetic mats. Geobiology. 2010;8:45–55.CAS 
PubMed 
Article 

Google Scholar 
Zeebe RE, Wolf-Gladrow D. CO2 in seawater: equilibrium, kinetics, isotopes. Amsterdam: Elsevier Science B.V.; 2001.Rae BD, Förster B, Badger MR, Price GD. The CO2-concentrating mechanism of Synechococcus WH5701 is composed of native and horizontally-acquired components. Photosynth Res. 2011;109:59–72.CAS 
PubMed 
Article 

Google Scholar 
Rippka R, Deruelles J, Waterbury JB, Herdman M, Stanier RY. Generic assignments, strain histories and properties of pure cultures of cyanobacteria. Microbiology. 1979;111:1–61.Article 

Google Scholar 
Martín-Cuadrado A-B, López-García P, Alba J-C, Moreira D, Monticelli L, Strittmatter A, et al. Metagenomics of the deep Mediterranean, a warm bathypelagic habitat. PLoS ONE. 2007;2:e914.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Bankevich A, Nurk S, Antipov D, Gurevich AA, Dvorkin M, Kulikov AS, et al. SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. J Comput Biol. 2012;19:455–77.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Hyatt D, Chen G-L, LoCascio PF, Land ML, Larimer FW, Hauser LJ, et al. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinforma. 2010;11:1.Article 
CAS 

Google Scholar 
Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2015;12:59–60.CAS 
PubMed 
Article 

Google Scholar 
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28:27–30.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Overbeek R, Olson R, Pusch GD, Olsen GJ, Davis JJ, Disz T, et al. The SEED and the Rapid Annotation of microbial genomes using Subsystems Technology (RAST). Nucleic Acids Res. 2013;42:D206–14.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Tatusov RL, Natale DA, Garkavtsev IV, Tatusova TA, Shankavaram UT, Rao BS, et al. The COG database: new developments in phylogenetic classification of proteins from complete genomes. Nucleic Acids Res. 2001;29:22–8.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Haft DH, Loftus BJ, Richardson DL, Yang F, Eisen JA, Paulsen IT, et al. TIGRFAMs: a protein family resource for the functional identification of proteins. Nucleic Acids Res. 2001;29:41–3.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Kang D, Li F, Kirton E, Thomas A, Egan R, An H, et al. MetaBAT 2: an adaptive binning algorithm for robust and efficient genome reconstruction from metagenome assemblies. PeerJ Prepr. 2019;7:e27522v1.
Google Scholar 
Parks DH, Imelfort M, Skennerton CT, Hugenholtz P, Tyson GW. CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res. 2015;25:1043–55.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Parks DH, Chuvochina M, Waite DW, Rinke C, Skarshewski A, Chaumeil P-A, et al. A standardized bacterial taxonomy based on genome phylogeny substantially revises the tree of life. Nat Biotechnol. 2018;36:996–1004.CAS 
PubMed 
Article 

Google Scholar 
Asnicar F, Thomas AM, Beghini F, Mengoni C, Manara S, Manghi P, et al. Precise phylogenetic analysis of microbial isolates and genomes from metagenomes using PhyloPhlAn 3.0. Nat Commun. 2020;11:1–10.Article 
CAS 

Google Scholar 
Garczarek L, Guyet U, Doré H, Farrant GK, Hoebeke M, Brillet-Guéguen L, et al. Cyanorak v2.1: a scalable information system dedicated to the visualization and expert curation of marine and brackish picocyanobacteria genomes. Nucleic Acids Res. 2021;49:D667–76.CAS 
PubMed 
Article 

Google Scholar 
Erwin PM, Thacker RW. Cryptic diversity of the symbiotic cyanobacterium Synechococcus spongiarum among sponge hosts. Mol Ecol. 2008;17:2937–47.CAS 
PubMed 
Article 

Google Scholar 
Usher KM, Toze S, Fromont J, Ku J, Sutton DC. A new species of cyanobacterial symbiont from the marine sponge Chondrilla nucula. Symbiosis. 2004.Holtman CK, Chen Y, Sandoval P, Gonzales A, Nalty MS, Thomas TL, et al. High-throughput functional analysis of the Synechococcus elongatus PCC 7942 genome. DNA Res. 2005;12:103–15.CAS 
PubMed 
Article 

Google Scholar 
Chen M-Y, Teng W-K, Zhao L, Hu C-X, Zhou Y-K, Han B-P, et al. Comparative genomics reveals insights into cyanobacterial evolution and habitat adaptation. ISME J. 2021;15:211–27.PubMed 
Article 

Google Scholar 
Katoh K, Misawa K, Kuma K, Miyata T. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 2002;30:3059–66.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Price MN, Dehal PS, Arkin AP. FastTree 2–approximately maximum-likelihood trees for large alignments. PLoS ONE. 2010;5:e9490.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Rozen S, Skaletsky H. Primer3 on the WWW for general users and for biologist programmers. Methods Mol. Biol. 2000;132:365–86.CAS 
PubMed 

Google Scholar  More