Alexander Winterl

Sebastian Richter

Aymeric Houstin

Téo Barracho

Matthieu Boureau

Clément Cornec

Douglas Couet

Robin Cristofari

Claire Eiselt

Ben Fabry

Adélie Krellenstein

Christoph Mark

Astrid Mainka

Delphine Ménard

Jennifer Morinay

Susie Pottier

Elodie Schloesing

Céline Le Bohec

Daniel P Zitterbart

Remote sensing of emperor penguin abundance and breeding success

Emperor penguins (Aptenodytes forsteri) are under increasing environmental pressure. Monitoring colony size and trends of this Antarctic seabird relies primarily on satellite imagery recorded near the end of the breeding season, when illumination levels are sufficient to capture images, but colony occupancy is highly variable. To correct population estimates for this variability, we develop a phenological model that accurately predicts the number of breeding pairs and fledging chicks, as well as key phenological events such as arrival, hatching and foraging times, from as few as six data points from a single season. The ability to extrapolate occupancy from sparse data makes the model particularly useful for monitoring remotely sensed animal colonies where ground-based population estimates are very rare or unavailable.

INTRODUCTION

Emperor penguins (Aptenodytes forsteri) breed on land-fast sea ice, making them particularly vulnerable to the impact of global warming [START_REF] Barbraud | Emperor penguins and climate change[END_REF][START_REF] Forcada | Penguin responses to climate change in the Southern Ocean[END_REF][START_REF] Jenouvrier | The Paris Agreement objectives will likely halt future declines of emperor penguins[END_REF][START_REF] Jenouvrier | Projected continent-wide declines of the emperor penguin under climate change[END_REF][START_REF] Fretwell | Discovery of new colonies by Sentinel2 reveals good and bad news for emperor penguins[END_REF][START_REF] Trathan | The emperor penguin -Vulnerable to projected rates of warming and sea ice loss[END_REF]. For breeding, stable land-fast sea ice (or fast ice), which is sea ice anchored to land, ice shelves, or grounded icebergs, is required [START_REF] Fretwell | Discovery of new colonies by Sentinel2 reveals good and bad news for emperor penguins[END_REF][START_REF] Ainley | Impacts of cetaceans on the structure of Southern Ocean food webs[END_REF]. The extent of the Antarctic fast ice is predicted to rapidly decline in the coming decades (8), and the species is predicted to become extinct by the end of this century. A significant decline of sea ice could already be observed in 2023 [START_REF]Springing into summer | Arctic Sea Ice News and Analysis[END_REF]. Despite this, the majority of emperor penguin colonies remain insufficiently studied due to the remoteness and harsh environmental conditions of their habitat [START_REF] Fretwell | Discovery of new colonies by Sentinel2 reveals good and bad news for emperor penguins[END_REF][START_REF] Bohec | Polar Monitoring: Seabirds as Sentinels of Marine Ecosystems[END_REF]. New strategies to better understand this species and their responses to changing environmental conditions are urgently needed.

Of the 61 currently known emperor penguin breeding colonies [START_REF] Fretwell | Discovery of new colonies by Sentinel2 reveals good and bad news for emperor penguins[END_REF], ground truth population counts conducted during the winter, at regular, frequent intervals, are only available for the colonies at Pointe Géologie in Adélie Land (average population size of 3900 breeding pairs) [START_REF] Weimerskirch | Banding Recoveries and the Dispersal of Seabirds Breeding in French Austral and Antarctic Territories[END_REF][START_REF] Jenouvrier | Effects of climate change on an emperor penguin population: analysis of coupled demographic and climate models[END_REF] and Atka Bay in Dronning Maud Land (average population size of 8600 breeding pairs) [START_REF] Richter | A remote-controlled observatory for behavioural and ecological research: A case study on emperor penguins[END_REF]. Only from these closely-monitored colonies, precise parameters of the species' reproductive cycle and life cycle have been recorded.

Emperor penguins return to their breeding colony at the onset of the lightless Antarctic winter, between late March and early May, to begin their annual breeding cycle. After a period of courtship, copulation, and egg-laying lasting 6 to 10 weeks, females leave the colony to forage and replenish their body reserves at sea, while the males remain to incubate the single egg for an average of 64 days [START_REF] Stonehouse | The emperor penguin (Aptenodytes forsteri, Gray): I. Breeding behaviour and development[END_REF][START_REF] Prevost | Ecologie des Manchots Antarctiques[END_REF]. The birds limit their body heat loss during incubation by forming tight groups, so called huddles [START_REF] Stonehouse | The emperor penguin (Aptenodytes forsteri, Gray): I. Breeding behaviour and development[END_REF][START_REF] Prevost | Ecologie des Manchots Antarctiques[END_REF]. Chicks hatch in winter, i.e. between July and August, and are not thermally independent during the first 6 to 7 weeks of life, therefore one parent always stays with the chick while the other forages. After this brooding period, chicks are left alone at the colony, forming crèches with the other chicks, so that both parents can forage simultaneously to satisfy the chick's growing demands. Parents feed their chicks between 7 and 12 times until fledging between November and January [START_REF] Stonehouse | The emperor penguin (Aptenodytes forsteri, Gray): I. Breeding behaviour and development[END_REF][START_REF] Prevost | Ecologie des Manchots Antarctiques[END_REF], just before the land-fast sea ice begins to break up. The duration of these foraging trips declines from 15-29 days after hatching to less than 10 days before fledging [START_REF] Ancel | Foraging behaviour of emperor penguins as a resource detector in winter and summer[END_REF][START_REF] Kooyman | Diving Behavior of Emperor Penguins Nurturing Chicks at Coulman Island, Antarctica[END_REF][START_REF] Kirkwood | Seasonal change in the foraging ecology of emperor penguins on the Mawson Coast, Antarctica[END_REF][START_REF] Zimmer | Foraging movements of emperor penguins at Pointe Géologie, Antarctica[END_REF][START_REF] Watanabe | Activity Time Budget during Foraging Trips of Emperor Penguins[END_REF][START_REF] Houstin | Biologging of emperor penguins-Attachment techniques and associated deployment performance[END_REF].

Since the majority of colonies are not surveyed by ground-based observations, satellite-based surveys provide the bulk of available datasets for population size estimation. As the resolution of satellite imagery has improved over time, satellite-based surveys hold the greatest potential for estimating global populations and detecting trends [START_REF] Fretwell | Discovery of new colonies by Sentinel2 reveals good and bad news for emperor penguins[END_REF][START_REF] Labrousse | Quantifying the causes and consequences of variation in satellitederived population indices: a case study of emperor penguins[END_REF][START_REF] Fretwell | An Emperor Penguin Population Estimate: The First Global, Synoptic Survey of a Species from Space[END_REF]. However, even at the highest resolution currently feasible, which is on the order of 0.5 to 1 m/pixel, satellite-based surveys cannot yet provide measurements at the individual level, but rather estimates of the number of animals based on the area occupied by the colony.

The viability of using colony area to estimate abundance suffers from uncertainties introduced by the imaging process (satellite off-nadir, sun azimuth, and sun elevation [START_REF] Labrousse | Quantifying the causes and consequences of variation in satellitederived population indices: a case study of emperor penguins[END_REF]), the conversion of colony areas to numbers of individuals, and large fluctuations of colony occupancy due to the species' phenology. The conversion of colony area to numbers of individuals requires knowledge of the average area number density (number of animals per square meter), which is subject to hourly fluctuations, as the animals regulate their body temperature by huddling together or loosening up, depending on weather conditions such as temperature, wind speed, solar radiation, and humidity [START_REF] Labrousse | Quantifying the causes and consequences of variation in satellitederived population indices: a case study of emperor penguins[END_REF][START_REF] Gilbert | One for all and all for one: the energetic benefits of huddling in endotherms[END_REF][START_REF] Richter | Phase transitions in huddling emperor penguins[END_REF]. Moreover, the occupancy fluctuations due to the annual phenology pattern is modulated by many factors such as sea ice extent and prey availability, which both influences the duration of foraging trips and the foraging success [START_REF] Zimmer | Foraging movements of emperor penguins at Pointe Géologie, Antarctica[END_REF][START_REF] Massom | Fast ice distribution in Adélie Land, East Antarctica: interannual variability and implications for emperor penguins Aptenodytes forsteri[END_REF].

Due to polar night, satellite images are not available during the incubation stage (June to July), when the least variation in the numbers of individuals is expected. Instead, usable satellite images can only be obtained between September and January, when chicks and only a fraction of the adults are present at the colony. The low resolution of the images does not allow for distinguishing between adults and chicks. For these reasons, satellite-based surveys suffer from large uncertainties when estimating the number of breeding pairs. Moreover, the breeding success of a colony based solely on the number of surviving chicks can not be determined [START_REF] Labrousse | Quantifying the causes and consequences of variation in satellitederived population indices: a case study of emperor penguins[END_REF][START_REF] Fretwell | An Emperor Penguin Population Estimate: The First Global, Synoptic Survey of a Species from Space[END_REF]. Currently, satellite-based surveys may be sufficient to roughly gauge population size and long-term trends, but not to assess short-term selective forces impacting breeding success, unless the colony has completely disappeared [START_REF] Fretwell | Emperors on thin ice: three years of breeding failure at Halley Bay[END_REF].

The aim of this study is to develop a method to compensate for the uncertainties of satellitebased surveys and to provide an estimate of the annual number of breeding pairs as well as the annual breeding success of a colony, based on the colony area measured during the austral spring and summer (September to December). To achieve this, we combine data from ground-based manual counts and images with a mechanistic phenological model for emperor penguin population counts. The model parametrizes the phenological annual breeding pattern and, when applied to data, infers key phenological and breeding parameters such as arrival time at the colony for breeding, breeding cycle duration, number of breeding pairs, number of hatched chicks, and number of fledgling chicks. Furthermore, we demonstrate that the model remains accurate even when the model parameters are estimated from sparse population counts, and even when these counts are obtained over only a short time window. We then apply our model to satellite-based data to estimate animal numbers and breeding success.

RESULTS

Phenological model

We count the number of emperor penguin adults and chicks on a weekly basis at Atka Bay (AB, 70° 40'S, 8° 16'W) over 3 breeding seasons (2018 to 2020), and at Pointe Géologie (PG, 66° 40'S, 140° 01'E) over 10 breeding (2012 to 2021). The data reveal a characteristic pattern of animal counts depending on annual cycles of incubation, guarding, foraging at sea, and feeding chicks at the colony (Fig. 1).

Based on known patterns of the species' breeding activity [START_REF] Ancel | Foraging behaviour of emperor penguins as a resource detector in winter and summer[END_REF][START_REF] Kooyman | Diving Behavior of Emperor Penguins Nurturing Chicks at Coulman Island, Antarctica[END_REF][START_REF] Kirkwood | Seasonal change in the foraging ecology of emperor penguins on the Mawson Coast, Antarctica[END_REF][START_REF] Zimmer | Foraging movements of emperor penguins at Pointe Géologie, Antarctica[END_REF][START_REF] Watanabe | Activity Time Budget during Foraging Trips of Emperor Penguins[END_REF][START_REF] Houstin | Biologging of emperor penguins-Attachment techniques and associated deployment performance[END_REF], we develop a mechanistic phenological model that describes the observed colony abundance during a typical breeding season (from March 1 to February 28 of the next year). This model describes the presence/absence pattern for a breeding pair and its chick based on the following set of parameters: time of arrival at the colony site, courtship duration, first female absence duration after laying, foraging trip duration, and period at the colony to feed the chick. The model then calculates the total number of females, males, and chicks within the colony at any given time point from the presence/absence patterns. We fit the model parameters to the observed number of adults, using a Markov chain Monte Carlo approach, which iteratively optimizes the parameters and returns a distribution for each model parameter as well as a distribution for the number of chicks and adults for each day of the breeding season.

To validate the model, we compare the observed numbers of individuals to the model predictions (see Fig. 2). Because the count data usually vary on a logarithmic scale, we report not arithmetic, but geometric errors: for example, a ±25% geometric error corresponds to an overestimate of 25% (multiply by 1.25) or an underestimate of 20% (divide by 1.25). We find an average geometric error of 16% for all data points (AB: 23%, PG: 15%). 81% of the observations fall within one standard deviation around the model estimate (AB: 83%, PG 81%). The coefficient of determination (R 2 value) between model and data is 0.91 for all data points (AB: 0.73, PG: 0.95). The model also provides date estimates for phenological events such as arrival, first departure and first return of the females after laying, chick emancipation (when the chick is thermally independent to join the other chicks in crèches), and fledging (when the chick has molted and leaves the colony for the first time). We compare these estimates with ground-based observations for Pointe Géologie (Fig. 3C, Table S11). The timing of events as estimated by the model shows good agreement with the observed timing of events (R 2 =0.98, average absolute error of 10 days) when pooled over all years. For individual years, we find good agreement of the model prediction for all events. Note that a difference of 14 days arises between model prediction and ground truth times for female first departure and chick emancipation, because these events are recorded in the field as the time of the first observation of the respective behavior, whereas the model predicts the central moment of the event.

Predictive capacity of the phenological model

A key feature of our model is its ability to predict the annual breeding success expressed as the number of fledging chicks relative to the annual number of breeding pairs. Because we only fit the model parameters based on the number of adults, we can assess the predictive power of the model by comparing the model estimate of the number of chicks with the actual number of chicks observed. We find an overall R² value of 0.45, with R²=0.64 for Pointe Géologie and R²=0.12 for Atka Bay (Fig. 2). The average geometric error is 41% for all data, with 51% for AB data and 39% for PG data. The large relative error for Atka Bay is attributable to the unusually (compared to other AB and PG seasons) high number of chicks that were still present late in the 2018 season (see Fig. 2K), which the model fails to predict. Further, the average number of sample counts per season for Atka Bay (26 counts) is lower than for Pointe Géologie (44 counts). In addition to the number of alive adults and chicks per day, the total number of fledging chicks, dead chicks and lost eggs are available from Pointe Géologie but not from Atka Bay. For Pointe Géologie, the model predicts the number of fledging chicks and dead chicks with R² values of 0.74 and 0.32, and average geometric errors of 25% and 56%, respectively (Fig. 3A,B). The model predicts the number of dead eggs with a mean absolute error of +/-190 animals.

Phenological parameters

The model describes each season and each colony separately. Therefore, we extract 13 sets (from 13 breeding seasons) of 14 parameters (BP, H, F, NB, t0, Δt0, m, b, Δb, cmax, cmin, smax, sfem, smin) that provide insight into the interannual and inter-colony variation of the phenology (Fig. S2 in the Supplement).

Most prominently, the arrival time (t0) at Pointe Géologie (April 7 +/-3 d) is significantly (P<0.01) earlier than at Atka Bay (April 27 +/-7 d). The number of breeding pairs (BP) also shows significant (P<0.01) differences between those two colonies (AB: 8600 breeding pairs on average, PG: 3900 breeding pairs on average). In addition, we find statistically significant (see S8 in the supplement) interannual variance for the number of breeding pairs (BP, AB: +/-620), fledging success (F, PG: +/-0.66), time of arrival (t0, AB: +/-12.0 d, PG: +/-4.6 d), width of arrival date distribution (Δt0, PG: +/-3.8 d), width of female return date distribution (Δb, PG: +/-7.0 d), minimum time at sea (smin, PG: +/-6.4 d), and maximum time at sea (smax, +/-6.7 d). Note, however, that the standard deviation for the duration of each event is smaller than the sampling interval of 7 days for the ground truth data. Consequently, the biological significance of these interannual variations cannot be reliably verified.

We further investigate the correlation between model parameters versus breeding success, as estimated by the ground truth ratio of fledging chicks to breeding pairs. We find significant (P<0.02, n=13 seasons) correlations between breeding success and the foraging trip duration during the crèching period with the following R² values: 0.54 (smax), 0.52 (smin), 0.46 (cmin) (see Table 1, Fig. S5 in the Supplement). We then sum the duration of all foraging/feeding trips to obtain the total time spent at sea or at the colony during the crèching phase (Fig. 4). We find a significant (P<0.03) correlation between breeding success and total time spent at sea (R²=0.55) or time spent at the colony (R²=0.39). Noteworthy, a longer total time spent at the colony correlates with a higher breeding success, while a longer total time spent at sea correlates with a lower breeding success. . We find a significant (n=13, P<0.03) positive correlation between breeding success and time at colony, and a significant (n=13, P<0.01) negative correlation between breeding success and time at sea.

Density estimation from ground-based imagery

To estimate penguin abundance and breeding success based on ground-based imagery, we need to multiply the measured colony area (ground area covered by penguins) with an estimate of the colony density (number of animals per area) to obtain the number of individuals. To estimate the colony density, we use a so-called windchill model that predicts colony density from locally measured meteorological variables, specifically air temperature, wind speed, solar radiation, and relative humidity, as described in [START_REF] Richter | Phase transitions in huddling emperor penguins[END_REF]. All of these meteorological variables affect the animals' heat balance and contribute to an apparent or perceived temperature, analogous to the windchill-corrected temperature reported in US-American weather forecasts. If the perceived temperature falls below a critical temperature, the animals tend to form huddles. Accordingly, the model requires a critical temperature as a further model parameter. Specifically, the critical temperature corresponds to the apparent temperature at which the average density of the animals within a colony reaches half of the maximum packing density in a dense huddle of 12.8 animals per m 2 .

To obtain the windchill model parameters, we manually select the colony covered area from ground-based images of Pointe Géologie and Atka Bay (13, 28) using the image annotation software Clickpoints [START_REF] Gerum | ClickPoints : an expandable toolbox for scientific image annotation and analysis[END_REF], project the resulting polygons to the top view projection [START_REF] Gerum | CameraTransform: A Python package for perspective corrections and image mapping[END_REF], and calculate the total occupied area in m² (Fig. 5A &5B). From this, we compile a dataset of 538 measurements of colony area, number of individuals (chicks and adults), and meteorological variables (temperature, wind speed, solar radiation, and humidity) acquired at each colony between September and December of three (AB) and four (PG) seasons. From the correlation of colony density fluctuations (colony area divided by the number of individuals) with fluctuations in meteorological variables, we obtain the windchill model parameters (Fig. 5C-F, Fig. S17).

The windchill model predicts the measured colony density with an R² value of 0.32 and an average geometric error of 39%. This error may seem large, but when we multiply the density predicted from the model with the measured area in order to obtain a predicted number of individuals, this accurately matches the actual animal counts with an R² value of 0.93 and an average geometric error of 11%. Moreover, we find that 60% of the individual count data fall within the predicted 1-sigma interval (Fig. 6F). Fig. 5H and5I show an example of the measured colony area and predicted number of individuals at Pointe Géologie in 2014.

The windchill model allows us to correct for weather-induced fluctuations in colony density, which is important when computing the number of animals from colony area. Fig. 5H shows fluctuations in colony area for Pointe Géologie between September 1 and December 31. Accordingly, we find large fluctuations of colony density over the course of a day (e.g. from 0.20 to 2.41 animals per m² on 2014-10-08, with a total range between 0.03 and 7.00 animals per m² during the whole observation period; Fig. 5G). These large density fluctuations explain the uncertainties of satellite-based counts when a constant density of 0.93 animals per m 2 (23) is assumed.

Application to satellite imagery

To showcase an application, we use previously published data of colony area from the population of Coulman Island (CI), Atka Bay (AB), and Stancomb-Wills Glacier (SW) in 2011 [START_REF] Labrousse | Quantifying the causes and consequences of variation in satellitederived population indices: a case study of emperor penguins[END_REF]. We obtain the local meteorological data at the time of the satellite image acquisition from [START_REF] Hersbach | The ERA5 global reanalysis[END_REF] and, based on the windchill model, predict a colony density between 0.6 animals/m² to 1.8 animals/m². After multiplying the colony area with the colony density, we obtain animal counts, which are the sum of chicks and adults present at the colony. Finally, we fit the phenological model to the counts to obtain the number of breeding pairs and the number of fledging chicks for each year. The results are summarized in Table S13. For these three colonies, low quality ground truth data for the number of breeding pairs exist, albeit not from 2011 but from earlier years [START_REF] Fretwell | An Emperor Penguin Population Estimate: The First Global, Synoptic Survey of a Species from Space[END_REF]. Nonetheless we find an average geometric error of 28% and an R² score of 0.88.

We also benchmark this approach using ground-based images from AB and PG. The procedure is analogous: we estimate the colony area from the images, use the windchill model to estimate colony density based on meteorological data from station observatories, multiply both numbers to estimate the number of animals, fit the phenological model to the counts, and arrive at an estimate for the number of breeding pairs and fledging chicks for each year that we can compare with ground truth counts. We find good agreement, with an average geometric error of 21%, an R² value of 0.92, and 74% of the ground truth data being within one standard deviation around the model estimate (Fig. 6A-F). Ground truth values for Pointe Géologie and Atka Bay (Ground) are from manual counts, all others are from [START_REF] Fretwell | An Emperor Penguin Population Estimate: The First Global, Synoptic Survey of a Species from Space[END_REF]. Note that for satellite data, there is no ground truth for the number of fledging chicks.

DISCUSSION

Phenological model performance

We present a mechanistic model that estimates phenological information from incomplete time series of animal counts in emperor penguin colonies. The phenological information predicted by the model include time of first arrival at the colony site to breed, number of breeding pairs, number of hatched chicks, number of fledged chicks, timing of foraging trips, duration of courtship, and female absence (during incubation). Our phenological model is able to recapitulate the temporal fluctuations of the weekly individual counts at two colonies (Atka Bay and Pointe Géologie) over several seasons, and it accurately predicts the number of chicks per day and the number of fledging chicks solely based on the counts of adult animals.

The model accurately predicts the number of breeding pairs and number of fledging chicks for both colonies. Predicted and true counts of chicks and adults per sampling day align less well at Atka Bay than at Pointe Géologie, probably due to a lower number of sample counts for Atka Bay than for Pointe Géologie (26 vs 44 counts per season). Our model is therefore robust in detecting the main features of the phenological pattern, while struggling with very fine predictions like weekly counts, when the underlying data are sparse. We attribute this to the probabilistic inference method that blurs the parameter space when confronted with noisy and sparse data.

The model has difficulty in years with highly unusual presence/absence patterns at specific times in the breeding cycle. For example, in 2017 at Pointe Géologie, the number of adults present during the chick rearing period was much lower than the number of chicks. The model (which is informed only about the number of adults but not the number of chicks) therefore estimated a survival rate of less than 50%, while in reality, significantly more chicks survived (71%). We hypothesize that, in 2017, adults spent extremely long periods outside the colony, reflecting an extended sea ice cover at the end of chick rearing [START_REF] Labrousse | Landfast ice: a major driver of reproductive success in a polar seabird[END_REF], while food resources were likely available and abundant near the sea ice edge, enabling them to feed their chick sufficiently until they fledge. Moreover, the model is based on a very limited knowledge and data on the individual trip durations during the chick-rearing phase (i.e. the study of [START_REF] Kirkwood | Seasonal change in the foraging ecology of emperor penguins on the Mawson Coast, Antarctica[END_REF], which was carried out over a single year on the Auster and Taylor Glacier and Auster colonies). To refine the model, it would be important to collect foraging information on this critical period of the breeding cycle, over several years and on several colonies facing contrasting environmental conditions.

Overall, the model performs well at predicting the number of fledging chicks and breeding pairs, although it gives large numbers of dead chicks and lost eggs compared to manual counts (Fig. 3A). A part of the discrepancy arises because lost eggs and dead chicks are counted only if they can be found, but are often covered by snow. Therefore, the model may in fact provide more reliable estimates for the number of dead chicks and lost eggs in years with typical phenology.

Timing in phenological parameters

The model performs well for predicting phenological events on average. We could not quantify how well the model can predict annual variations in phenological events at the same colony, because most events vary less than the resolution of our observations (7 days, see Fig. 3C). Nonetheless, some interesting model predictions emerge. For example, the duration of courtship (m) and the females' first absence from the colony after laying (b) (the incubation period for the males), do not show significant variability between years and colonies, which is to be expected because such processes have likely evolved to be tightly regulated at a physiological level specific to the species. We also found a negative correlation of the total time spent at sea with breeding success: long (>10 days) foraging trips at the end of the chick-rearing phase (between September to December) lead to lower (<50%) breeding success, likely because the frequency of feeding events decreases, exposing the chick to a higher risk of starvation in between feeding events. We suggest that increased foraging trip durations are linked to longer distances to the land-fast ice edge, which have been shown to negatively impact fledging success [START_REF] Labrousse | Landfast ice: a major driver of reproductive success in a polar seabird[END_REF].

The model-predicted as well as measured arrival times (t0) and first hatchings (tH=t0+m+b) between Pointe Géologie and Atka Bay differ by approximately 30 days. At Pointe Géologie, the more northerly colony (66.39°S), first hatchings are observed in the first week of July, while at Atka Bay (70.60°S), hatchings are estimated (by our model) to occur in the first two weeks of August. In both colonies, the timing of hatchings correlates with the first sunrise after mid winter (June 29 at Pointe Géologie and July 28 at Atka Bay, Fig. S18). A positive correlation between arrival times to breed and latitudes has been previously reported [START_REF] Williams | The Penguins: Spheniscidae[END_REF]. A likely explanation of this correlation is that penguins synchronize the chick rearing phase and therefore their whole annual cycle with the increasing abundance of prey. Prey abundance, in turn, is triggered by a rise in primary production due to an increase in solar radiation and/or a decrease in sea ice cover near the Antarctic continent. A similar synchronization between prey abundance and breeding cycle has been observed in numerous species, including other seabirds [START_REF] Ramírez | Sea ice phenology and primary productivity pulses shape breeding success in Arctic seabirds[END_REF], and is in line with the match/mismatch hypothesis [START_REF] Durant | Climate and the match or mismatch between predator requirements and resource availability[END_REF][START_REF] Cushing | Plankton Production and Year-class Strength in Fish Populations: an Update of the Match/Mismatch Hypothesis[END_REF][START_REF] Cushing | The Regularity of the Spawning Season of Some Fishes[END_REF] that states that the reproductive success of a species depends on its ability to match its phenology with that of its prey species.

Windchill model

We find a relationship between colony density measured from ground-based images and meteorological conditions. This relationship is described by a previously reported model (which for simplicity we refer to as the windchill model) [START_REF] Richter | Phase transitions in huddling emperor penguins[END_REF]. We re-trained the windchill model with data from Atka Bay and Pointe Géologie obtained between September and December, as the previously reported model parameters were trained with data obtained during the high-density (up to 15 animals per m²) courtship and incubation phase of the austral winter, in April and May.

The re-trained model can only partially explain the fluctuations in colony density with a correlation of R²=0.32 and 60% of data within the predicted 1 sigma interval, which is slightly higher than the percentage of the original model (50%, ( 13)). This is to be expected because of the mild weather conditions of the austral spring, and because a mixed colony of adults and chicks rarely adopts a configuration where the majority of the animals are found in dense huddles.

Estimating abundance and breeding success from remote sensing

Satellite images can only be taken during periods when occupancy is most erratic and therefore the most challenging for inferring population size and success. One of the main objectives of our phenological model is to provide more reliable estimates of the number of breeding pairs within a colony from those images. Currently, it is common practice to multiply the colony area A (in units of m 2 ) as measured from satellite images with a conversion factor CF of 0.93 breeding pairs per m 2 [START_REF] Fretwell | An Emperor Penguin Population Estimate: The First Global, Synoptic Survey of a Species from Space[END_REF]. However, depending on the weather conditions and the time of year, the colony area and animal density fluctuate, potentially contributing to large uncertainties [START_REF] Labrousse | Quantifying the causes and consequences of variation in satellitederived population indices: a case study of emperor penguins[END_REF][START_REF] Fretwell | An Emperor Penguin Population Estimate: The First Global, Synoptic Survey of a Species from Space[END_REF]. To reduce these uncertainties, we suggest the following approach for a reliable estimate of the number of breeding pairs: (1) Collect several (5 to 10) satellite images from multiple time points t (over a time period of 2 months, from October to December) and measure the colony area A; (2) Acquire local weather conditions at the time of satellite image capture and estimate colony density 𝜚 using the windchill model;

(3) Compute the total number N of animals in the colony (N = 𝜚 A); ( 4) Fit the phenological model to N(t) and extract the number of breeding pairs BP and fledged chicks F.

To benchmark this method, we estimate the number of breeding pairs from ground-based images of the colony (instead of satellite images) for our two colonies (AB and PG) over a time period of 3 and 4 years, and find an excellent correlation (R²=0.93). Note that if multiple images are recorded on a single day, the predicted counts are averaged in order to reduce sampling bias. Furthermore, we compute a conversion factor CF from the number of breeding pairs and the colony area as CF = BP/A (see S15 in the Supplement). We obtain a value of CF = 1.09 breeding pairs per m², which is close to the conversion factor of 0.93 breeding pairs per m² as reported in [START_REF] Fretwell | An Emperor Penguin Population Estimate: The First Global, Synoptic Survey of a Species from Space[END_REF]. Note that when we compute the average CF, we exclude an outlier from the 2014 season at Pointe Géologie with an unusually low (<5%) breeding success, a very low number of animals in the colony (approximately 500 instead of the usual 1800 animals on average), and a very high density of 5.98 breeding pairs per m². However, this outlier is not excluded when we compute the correlation of the model estimates with manual counts. Taken together, we are convinced that previous estimates of the total number of breeding pairs from satellite images are accurate on average, but individual counts from any given image can be subject to large errors, as was already pointed out in [START_REF] Labrousse | Quantifying the causes and consequences of variation in satellitederived population indices: a case study of emperor penguins[END_REF][START_REF] Fretwell | An Emperor Penguin Population Estimate: The First Global, Synoptic Survey of a Species from Space[END_REF]. The number of surviving chicks that we estimate with our model following the approach outlined above also agrees well with manual counts (R²=0.81 and average geometric error of 26% based on data from two colonies over 3 or 4 years).

We find a rather striking difference in the number of breeding pairs at Atka Bay between the periods 2008-2011 and 2018-2020. In 2008-2011, the estimates range from 7300 (satellite images, Table S13) to 9657 (satellite images, ( 23)), while in 2018-2020, the estimates are around 10,000 (both from ground truth counts and ground-based images). This increase by more than 2,000 breeding pairs is likely due to an immigration of individuals from the neighboring Halley colony that experienced three consecutive years of breeding failure (from 2014 to 2016) and had completely abandoned the colony site in 2016 [START_REF] Fretwell | An Emperor Penguin Population Estimate: The First Global, Synoptic Survey of a Species from Space[END_REF].

While the number of breeding pairs provides valuable information about the current status and health of the species, the number of surviving chicks during a breeding period is a more sensitive predictor of future population trends, e.g. due to global change or changes in food supply. As of now, we have not detected a declining trend in the number of surviving chicks at AB or PG over the 2012-2021 period. However, given the recent development of a dramatically decreased sea ice extent around Antarctica (9), chances for the long term survival of the species stay dire.

The next milestone will be to apply our method on a long-term circum-Antarctic scale, enabling us to use the emperor penguin breeding success as an early warning indicator. This indicator will alert the scientific community on an annual level about potential breeding failures, allowing us to better correlate such breeding failures with local oceanographic events. Ultimately we aim to provide stakeholders and governments with information on the health of populations and their ecosystems to rapidly implement conservation measures.

MATERIALS AND METHODS

Individual counts

The individual counts were conducted at the Pointe Géologie (PG) emperor penguin colony in vicinity to the Dumont d'Urville French research station (66° 40'S, 140° 01'E) and at the Atka Bay (AB) colony in the vicinity to the Neumayer III German research station (70° 40'S, 8° 16'W). The data from Pointe Géologie were collected over 10 years (2012 to 2021), Atka Bay data were collected over 3 years (2018 to 2020) (Fig. 1). Weekly counts were acquired by imaging the colonies with either hand-held (Pointe Géologie) or remote-controlled cameras (Atka Bay, ( 13)) from an elevated position, resulting in one or several panoramic images per day that have sufficient resolution to count individual penguins. For manual counting, we used Adobe Photoshop or Clickpoints [START_REF] Gerum | ClickPoints : an expandable toolbox for scientific image annotation and analysis[END_REF].

Manual observations of phenological events such as time of arrival, female departure, or the beginning of fledging were recorded only at Pointe Géologie. However, some events (e.g. beginning of fledging, see Table S11 in the supplement) could not be recorded in some years for various reasons such as inaccessibility of the colony, occlusion of the colony by geographic landmarks, or a low number of breeding pairs.

The ground truth numbers of breedings pairs at AB and PG were taken from individual counts during the incubation period when only the incubating males are at the colony. Lost eggs/chicks were counted daily at PG. The number of fledged chicks was not counted. Therefore, we used the last known number of chicks before the onset of fledging as an estimate of the number of surviving chicks.

Colony area and colony density

Colony area was measured from images after perspective correction. The Pointe Géologie colony was monitored between 2014 and 2017 using automatic time lapse cameras [START_REF] Winterl | micrObs -A customizable time-lapse camera for ecological studies[END_REF]. For Atka Bay, we used the panoramic images from the seasons 2018, 2019, and 2020. We manually marked the colony boundaries on 538 images (509 of Pointe Géologie, 29 of Atka Bay) acquired between 1st of September and 31st of December of each season. An example of an annotated image is shown in Fig. 5A. From the colony boundaries we computed the surface area covered by the colony (in m²) by perspective correction and projection using the intrinsic and extrinsic camera parameters such as focal length, elevation, tilt, and roll [START_REF] Gerum | CameraTransform: A Python package for perspective corrections and image mapping[END_REF]. A projected top view is shown in Fig. 5B. The extracted areas of the Pointe Géologie colony between September 1 and December 31, 2014 are shown in Fig. 5H. Colony densities (average number of individuals (adults plus chicks) per area) for Atka Bay and Pointe Géologie are calculated from individual counts (weekly resolution) divided by colony area for each of the images. For time points where we had area measurements but not a corresponding individual count, we linearly interpolated the value from the two nearest available counts. We extended our dataset of colony area with satellite image based measurements of colony area from Coulman Island (12 days), Atka Bay (19 days), and Stancomb-Wills Glacier (18 days), between September 10 and December 11, 2011 published in [START_REF] Labrousse | Quantifying the causes and consequences of variation in satellitederived population indices: a case study of emperor penguins[END_REF].

Meteorological data

The meteorological data for Pointe Géologie and Atka Bay (seasons 2018-2020) were recorded every minute by the meteorological observatory at Dumont d'Urville station (operated by Météo France) and the meteorological observatory at Neumayer station. The meteorological data for Coulman Island, Atka Bay (season 2011), and Stancomb-Wills Glacier stem from the ECMWF model "The ERA5 global reanalysis", which provides 1h temporal and 1km² spatial resolution [START_REF] Hersbach | The ERA5 global reanalysis[END_REF] and is available online (https://www.ecmwf.int/en/forecasts/dataset/ecmwf-reanalysis-v5).

Phenological model

We developed a mechanistic phenological model to describe the temporary fluctuations of the number penguins at the colony over the course of a breeding season.

The most important model parameter is the number breeding males and females. This number is assumed to remain constant during a breeding cycle but can change between years. Additionally, the colony consists of a number of non-breeding animals. The model then describes the temporary fluctuations of the probability that male and female breeding animals or non-breeding animals are present at the colony. These probabilities are computed from the following phenological events: -arrival of the animals (breeders and non-breeders) at the colony site at the beginning of the breeding season, -departure of the females after courtship and laying an egg, and departure of the nonbreeders after unsuccessful courtship, -return of the females for chick feeding and departure of the males for foraging at sea after hatching of the chicks, -subsequent (2x, according to (18)) parental switching between feeding and foraging, departure of both parents after their chick becomes thermally independent, -subsequent (7x, according to ( 18)) feeding and foraging trips (females and males) until fledging of the chicks. Note, that the number of foraging trips (2x during guarding phase and 7x during crèching phase) represent the average over the whole colony. We could not introduce the trip number as a free parameter of the model without losing numeric stability. In the model, these events are described as probability densities with Gaussian shape that sets the mean time point and the standard deviation (temporal spreading) of the event (see Fig. 7A). Each event is specific for male, female, or non-breeding animals. The amplitude of the probability density for each event is chosen so that its integral corresponds to the total number of all animals of a specific sex (or of the total number of all non-breeding animals) that participate in the event. For example, as chicks die, the number of parents returning to the colony for feeding decreases, which we model by a smaller amplitude of the probability density for events describing the feeding trips. The number of male, female, or non-breeding animals is computed from the integration of the event probabilities over time (see Fig. 7B &7C). The number of chicks is computed from the integration of the female event probabilities over time, starting with their first return after incubation.

The model parameters (time point, standard deviation, and amplitude of each event) are determined by a least squares fit of the model predictions to data. Since observers usually cannot distinguish between male and female animals, the sum of all adult animals is considered. For fitting the model, the number of chicks is not considered, however if this number is available from observations, it can be used for benchmarking the predictive power of the model.

The time points of the two parental switching (between feeding at the colony and foraging at sea) events for each sex, and the 7 consecutive feeding and foraging trips, are not free parameters. Rather, they are determined by the duration of a foraging trip at sea and the duration of the stay at the colony for feeding the chick. Furthermore, we assume that both the foraging and feeding durations decrease over time, in line with field observations [START_REF] Ancel | Foraging behaviour of emperor penguins as a resource detector in winter and summer[END_REF][START_REF] Kooyman | Diving Behavior of Emperor Penguins Nurturing Chicks at Coulman Island, Antarctica[END_REF][START_REF] Kirkwood | Seasonal change in the foraging ecology of emperor penguins on the Mawson Coast, Antarctica[END_REF][START_REF] Zimmer | Foraging movements of emperor penguins at Pointe Géologie, Antarctica[END_REF][START_REF] Watanabe | Activity Time Budget during Foraging Trips of Emperor Penguins[END_REF][START_REF] Houstin | Biologging of emperor penguins-Attachment techniques and associated deployment performance[END_REF], which we model by two linear functions (see Fig. 8). Furthermore, the duration of the first foraging trip of the females is assumed to be shorter than the first foraging trip of the males, also in line with field observations [START_REF] Kirkwood | Seasonal change in the foraging ecology of emperor penguins on the Mawson Coast, Antarctica[END_REF][START_REF] Williams | The Penguins: Spheniscidae[END_REF].

Similarly, the standard deviations of all events are not free parameters but are constrained as follows: the standard deviation of arrival of all adult penguins, and the standard deviation of the first departure of females and non-breeders are the same, assuming that courtship is equally long for all breeders. The standard deviation of the first return of females and of all subsequent events are also equal and must be greater or equal than the standard deviation of arrival.

Taken together, the model has 14 free parameters (Table 1). We use Markov chain Monte Carlo (MCMC, see S1 in the Supplement) sampling to fit the model parameters to data (ground truth counts from manual observations). The best estimates of the model parameter for each colony and each season are listed in Table S3 in the Supplement. 18)). Red regression lines (A,B) show the decreasing trend in both time at sea and at the colony, after the initial brooding phase. (C) and (D) show the predicted time at sea and time at colony as predicted by the phenological model for 10 seasons (2012-2021) for Pointe Géologie and 3 seasons (2018-2020) for Atka Bay. Red, blue and gray dots denote the sex of the observed individuals (female/male/unknown). Gray lines indicate corresponding predictions for each season. Note that the model predicts large seasonal variations for the time at sea, but not for the time at colony.

Windchill model

While they incubate their egg during the harsh Antarctic winter, emperor penguins conserve energy by forming tight groups, the so-called huddles [START_REF] Gilbert | One for all and all for one: the energetic benefits of huddling in endotherms[END_REF][START_REF] Maho | The Emperor Penguin: A Strategy to Live and Breed in the Cold: Morphology, physiology, ecology, and behavior distinguish the polar emperor penguin from other penguin species, particularly from its close relative, the king penguin[END_REF]. The fraction of individuals of a colony that are currently in a huddle changes depending on an apparent (i.e. subjectively perceived) temperature. The apparent temperature depends on four environmental variables: ambient temperature (T), windspeed (W), solar radiation (R), and humidity (H). These variables linearly contribute to the apparent temperature Ta , with linear factors cW (windchill factor), cR (solar radiation factor), cH (humidity factor), and a factor of unity for ambient temperature (Eq. 1) [START_REF] Richter | Phase transitions in huddling emperor penguins[END_REF].

The model then describes the colony density as a sigmoidal function of the apparent temperature (Eq. 2). For very low apparent temperatures, the colony density tends to a maximum value of 12.8 animals / m 2 . This maximum density corresponds to a hexagonal packing of cylinders with a diameter of 30 cm. At very high apparent temperatures, the colony density tends to zero. The sigmoidal function has two free parameters, the critical temperature Tc at which the density is half of its maximum, and the steepness value b0 that describes how sensitively the animals respond to temperature changes.

Ta = T + cWW + cRR + cHH

(1)

N/A = ⍴ = 12.8 / (1 + exp(( Ta -Tc)/b0))

(

We fit the model parameters b0, cW, cH, cR, and Tc to the density data recorded at Pointe Géologie and Atka Bay between 2014 and 2020 for the months of September to December using Markov chain Monte Carlo sampling (see S1 in the Supplement). We can then predict the density of any colony at any given time point between September and December from weather data alone.

Applying the phenological model to data from satellite images

When we fit the phenological model to count estimates from satellite images, problems may arise due to the small number of usable images for each season (typically less than 10 images, available only during the months from September to December). The fit of 14 free model parameters to such a small number of data points can lead to numerical instability and large confidence intervals. Therefore, we fix all except three model parameters to the average values that we obtained from the model fit to ground truth counts at Atka Bay and Pointe Géologie (see Table S3 in the Supplement). The time of arrival is chosen so that the date of female return aligns with the first sunrise after mid winter + 27.4 d. The only two free fit parameters are the number of breeding pairs and fledging success. This approach is justified by the fact that the phenology of the two colonies we have studied here (AB and PG) are similar except for the arrival time, as seen by the similar values of the model parameters and the absence of large systematic fluctuations during different years.

of Dumont d'Urville. We are deeply grateful to all the wintering and summering members of projects IPEV 137, AWI-SPOT, AWI-MARE, and we also sincerely thank the IPEV and AWI logistics teams for their important and continued support in the field.

Fig. 1 .

 1 Fig. 1. Manual counts of the number of individuals. Number of individuals at Pointe Géologie (A,B) and Atka Bay (C,D) colonies. Circles indicate counts of individual adults (A,C), crosses indicate counts of individual chicks (B,D). Colors denote different breeding seasons. Gray bars show the time range of key phenological events extracted from manual observations at Pointe Géologie (A,B). Note that the number of chicks increases between August and November as the chicks become thermally independent and therefore less occluded.

Fig. 2 .

 2 Fig. 2. Phenological model fit to manual counts. Overlay of data and model output for all observed seasons at Pointe Géologie (10, A-J) and Atka Bay (3, K-M) colonies. Individual panel titles indicate colony (PG = Pointe Géologie, AB = Atka Bay) and season. Y-axes show the number of individual adults (circles) and chicks (crosses). Solid and dashed lines show the model prediction for the adult and chick counts, respectively. The shaded areas

Fig. 3 .

 3 Fig. 3. Model estimates of key phenological parameters. (A-B) Observed and predicted breeding success by year and colony. The solid black and dashed black bars show the true and predicted total number of breeding pairs. The colored areas split the bars into fledged chicks (green), dead chicks (red), and lost eggs (orange). The predicted numbers add up by model definition while the ground truth values do not, due to counting inaccuracy. The solid line separates between Pointe Géologie and Atka Bay data, and the left and right Y-axis. (C) Time of predicted phenological events vs observed phenological events. Each point

Fig. 4 .

 4 Fig.4. Correlation between foraging pattern and breeding success. The x-axis shows the total number of days spent at the colony (A) or at sea (B) during the crèching phase as predicted by the model, averaged for females and males. The y-axis shows the ratio of fledging chicks to breeding pairs (breeding success) from manual observations. Every dot (Pointe Géologie) and square (Atka Bay) denotes one season. The black lines show the regression lines for each point cloud (Pointe Géologie and Atka Bay merged). We find a significant (n=13, P<0.03) positive correlation between breeding success and time at colony, and a significant (n=13, P<0.01) negative correlation between breeding success and time at sea.

Fig. 5 .

 5 Fig. 5. Description of the windchill model. (A) Example of an image recorded by the automatic time lapse camera on 2017-08-22 04:00:00 UTC at Pointe Géologie. The manually-annotated outline of the colony is highlighted with a red polygon. (B) Projected top view of the image shown in (A). The number indicates the area covered by the colony.

Fig. 6 .

 6 Fig. 6. Application to satellite data.(A-E) Fit of the phenological model to individual counts inferred from measured colony area and the windchill model predictions (dots). Lines show the phenological model predictions for adults (dashed), chicks (dotted), and total counts (solid). In (A), adult and chick counts are omitted for readability. Those data are provided in Supplement S7. Colors denote seasons. (F) Comparison of ground truth counts and counts inferred from colony area and windchill model predictions for Pointe Géologie and Atka Bay over 7 seasons. Colors denote seasons. Note: Pointe Géologie data covers 2012-2017, Atka Bay 2018-2020, therefore colors are also colony specific. The back line shows identity. (G,H) Comparison of ground truth and phenological model predictions for the number of breeding pairs (G) and fledging chicks (H). Ground truth values for Pointe Géologie and Atka Bay (Ground) are from manual counts, all others are from[START_REF] Fretwell | An Emperor Penguin Population Estimate: The First Global, Synoptic Survey of a Species from Space[END_REF]. Note that for satellite data, there is no ground truth for the number of fledging chicks.

Fig. 7 .

 7 Fig. 7. Illustration of the breeding cycle and the phenological model (based on data from 2012 at Pointe Géologie). The model describes the number of individuals present at the

Fig. 8 .

 8 Fig. 8. Individual trip durations during chick rearing phase. Time (in days, mean ± 1 standard deviation) a breeding emperor penguin spends at sea (A) or at the colony (B) at different stages of the breeding cycle (data from (18)). Red regression lines (A,B) show the

Table 1 .

 1 Phenological Model Parameters. Parameter names, descriptions, and numerical ranges.

	Symb	Name	Lower	Upper	Unit Fit per	Description
	ol		Limit	Limit	year /
					global
	t0	arrival	50	150	days year	Peak of animal arrivals at the

Table 1 .

 1 Phenological Model Parameters. Parameter names, descriptions, and numerical ranges.

	Symb	Name	Lower	Upper	Unit Fit per	Description
	ol		Limit	Limit		year /
						global
		date	(Feb 1) (May 30)			colony
	Δt0	arrival	4	14	days global	Standard deviation of arrival dates
		date
		standard
		deviatio
		n
	m	courtshi	28	42	days global	Average duration between arrival
		p					and first departure of females
		duration
	b	female	50	100	days global	Average duration between first
		absence					departure and return of females
		duration
	Δb	female	0	14	days global	Standard deviation of duration of
		absence					female absence
		duration
		standard
		deviatio
		n
	BP	number	2.000	15.000	pairs year	Number of pairs that mate and
		of pairs					produce an egg.
	NB	non	0	1	1	year	Ratio of birds that do not find a
		breeder					breeding partner and leave the
		to					colony before the breeding period,
		breeding					to the total number of breeding
		pair ratio					pairs
	H	hatching	0	1	1	year	Ratio of pairs that repeatedly return
		success					after the breeding period, relative
		ratio					to the number of breeding pairs at
							the beginning of the breeding
							season
	F	fledging	0	1	1	year	Ratio of pairs that repeatedly return
		success					until the end of the chick rearing
		ratio					phase, relative to the number of
							breeding pairs at the beginning of
							the breeding season

Table 1 .

 1 Phenological Model Parameters. Parameter names, descriptions, and numerical ranges.

	Symb	Name	Lower	Upper	Unit Fit per	Description
	ol			Limit	Limit	year /
						global
	cmax	maximu	1d	5d	days year	Time one parent spends at the
		m time				colony per foraging trip during the
		at					guarding phase and at their first
		colony				foraging trip.
		per trip
	cmin	minimu	1d	cmax	days year	Time one parent spends at the
		m time				colony on their last foraging trip.
		at
		colony
		per trip
	smax	maximu	cmax	21d	days year	Time one parent spends away from
		m time				the colony per foraging trip during
		at	sea				guarding phase and at their first trip
		per trip				after chick emancipation. This
							includes foraging and commuting
							time
	smin	minimu	cmax	smax	days year	Time one parent spends away from
		m time				the colony per trip on their last
		at	sea				foraging trip. This includes
		per trip				foraging and commuting time
	sfem	time at	cmax	smax	days year	Duration of the first foraging trip of
		sea	at				females after hatching. It is shorter
		first trip				or equal than the following trips.
		(only
		females)

Acknowledgements

This work was supported by the Deutsche Forschungsgemeinschaft (DFG) in the framework of the priority programme 1158 "Antarctic Research with comparative investigations in Arctic ice areas" by grants (FA336/5-1, ZI1525/3-1, ZI1527/7-1), by the Institut Polaire Français Paul-Emile Victor (IPEV) within the framework of the Project 137-ANTAVIA, by the Alfred-Wegener-Institut Helmholtz-Zentrum für Polar-und Meeresforschung (AWI) within the framework of the Projects SPOT and MARE, by the Centre Scientifique de Monaco with additional support from the LIA-647 and RTPI-NUTRESS (CSM/CNRS-UNISTRA), by the Centre National de la Recherche Scientifique (CNRS) through the Programme Zone Atelier Antarctique et Terres Australes (ZATA). We thank Météo France for the meteorological data

Funding

German Research Foundation grant FA336/5-1 (BF) German Research Foundation grant ZI1525/3-1 (DZ) German Research Foundation grant ZI1527/7-1 (DZ) National Science Foundation grant 2046437 (DZ) Woods Hole Oceanographic Institution (DZ) Institut Polaire Français Paul-Emile Victor (IPEV) 137-ANTAVIA (CLB) Centre Scientifique de Monaco (CLB) Centre National de la Recherche Scientifique (CLB) CSM/CNRS-UNISTRA (LIA-647 and RTPI-NUTRESS) (CLB)

Data and materials availability:

All data are available in the main text or the supplementary materials.

Author Contributions

Conceptualization: AW, SR, BF, CLB, DZ Methodology: AW, SR, BF, AM, CM, CLB, DZ Data analysis: AW Data collection: AW, SR, AH, TB, MB, CC, DC, RC, CE, BF, AK, AM, DM, JM, SP, ES, CLB, DZ Supervision: BF, CLB, DZ Writing-original draft: AW, SR, AH, BF, CLB, DZ All authors contributed critically to the drafts and gave final approval for publication.

Competing interests:

Authors declare that they have no competing interests