DATA COLLECTION AND PROCESSING

Sea hare observations were extracted from iNaturalist using the rinat package (Barve & Hart 2014) in R version 4.3.0 (R Core Team 2023). All research-grade observations of the order Aplysiidae were queried from Australia (iNaturalist place_id = 6744) with a maximum return of 10,000 records. A data quality metric set by iNaturalist is ‘Research-grade’, which requires community agreement on species identification and the presence of date and location metadata was the primary data filter. Observations were further filtered to include only records with complete geographic coordinates (latitude and longitude) and observation dates were set for between 2010 and 2024 (which focuses on a period of substantial iNaturalist activity in Australia whilst maintaining a 15-year time series for trend analysis). Seasonal-trend decomposition using LOESS (STL) with robust fitting was applied to monthly observation counts to eliminate sampling bias associated with large growth in iNaturalist participation. Decomposition parameters were set at: time series span (Jan 2015 – Dec 2024); seasonal window – periodic; robust fitting method using iteratively weighted least squares; ~12 fold increase over 10 years; seasonal component captures recurring within-year patterns trend range ~10 (2015) to ~130 (2024); seasonal amplitude – ±10–15 observations/month; seasonal peak – Dec to Jan (consistent across years); and seasonal trough – May to Jun (consistent across years).

To evaluate whether the seasonal patterns reported here reflect biology or platform-usage artefacts, we extracted monthly aggregate counts of all Australian iNaturalist Research-Grade observations (2010–2024) from GBIF (2026), accessed via the rgbif package (Chamberlain & Boettiger 2017). This baseline (n = 5,389,576) provides a generalised measure of within-year variation in observer activity against which the seasonality of sea hare observations can be compared. If observer effort drove seasonal patterns, the monthly distributions of all-Australia observations and sea hare observations would coincide.

GEOGRAPHIC PARTITIONING

For analyses requiring spatial structure, observations were assigned to two complementary classifications. Coast: East (>140°E) and West (<125°E); a small number of observations from the Gulf of Carpentaria and central northern coast (125–140°E, n = 88, 1.3%) were excluded from coastal analyses due to insufficient sample size. Latitudinal zone: Tropical (<23°S), Subtropical (23° to 29°S), and Temperate (>29°S). Species-level analyses were restricted to taxa with ≥50 observations to ensure adequate sample sizes for statistical inference (Appendix: Table A1).

Objective 1 – Quantify seasonality on the Australian east and west coasts

Analyses used all observations partitioned by coast as defined above (n = 5,508 East, n = 1,033 West; 88 Gulf of Carpentaria/central northern observations excluded). Observations for a total of 15 taxa were used from five genera (Table 1): eight species of Aplysia Linnaeus, 1767 (A. argus Rüppell et Leuckart, 1830, A. concava G. B. Sowerby I, 1833, A. extraordinaria (J. K. Allen, 1932), A. gigantea G. B. Sowerby II, 1869, A. juliana Quoy et Gaimard, 1832, A. oculifera A. Adams et Reeve, 1850, A. reticulata Eales, 1960 and A. sydneyensis G. B. Sowerby II, 1869), two species of Bursatella Blainville, 1817 (B. hirsuta Nimbs et N. G. Wilson, 2020 and B. leachii Blainville, 1817), two species of Dolabrifera J. E. Gray, 1847 (D. brazieri G. B. Sowerby II, 1870 and D. dolabrifera (Rang, 1828)), one species of Stylocheilus A. A. Gould, 1852 (S. striatus (Quoy et Gaimard, 1832)) and the monospecific Dolabella Lamarck, 1801 (D. auricularia ([Lightfoot], 1786) and Syphonota H. Adams et A. Adams, 1854 (S. geographica (A. Adams et Reeve, 1850)). Nomenclature used in this study was validated using the World Register of Marine Species (WoRMS 2026).

Table 1

Dataset summary by species. N = number of research-grade iNaturalist observations per species by latitudinal zone. Tropical: <23°S; Subtropical: 23° to 29°S; Temperate: >29°S

Objective 2 – Test for species-specific differences in phenology (temporal niche partitioning)

Analyses used the 15 species with ≥50 observations (Figs 1–15, Table 1). Circular statistics were computed both for pooled species data and within latitudinal zones (as defined above) for each species (Agostinelli & Lund 2017). Watson’s U² compared seasonal distributions among zones within species. Spectral analysis of detrended monthly counts corroborated dominant 12‑month periodicity where it was informative.

Figs 1–15

Image of Australian sea hare taxa analysed in this study: 1 – Aplysia argus, Caloundra, QLD; 2 – Aplysia concava, North Solitary Island, NSW; 3 – Aplysia extraordinaria, Sandy Beach, NSW; 4 – Aplysia gigantea North Mole, WA; 5 – Aplysia juliana, Bare Bluff, NSW; 6 – Aplysia sydneyensis, Clifton Gardens, NSW; 7 – Aplysia reticulata, Port Headland, WA; 8 – Aplysia oculifera, Garden Island, WA; 9 – Bursatella hirsuta, Woodman Point, WA; 10 – Bursatella leachii, Sawtell, NSW; 11 – Dolabella auricularia, Diggers Headland, NSW; 12 – Dolabrifera brazieri, Woolgoolga, NSW; 13 – Dolabrifera dolabrifera, Woolgoolga, NSW; 14 – Stylocheilus striatus, Sandy Beach, NSW; and 15 – Syphonota geographica, Nelson Bay, NSW. Photographs: 1–3, 5–6, 9–15 by Matt John Nimbs; 4 by Tim Karnasuta* (iNaturalist), 7 by Ilze Keevey* (iNaturalist) and 8 by Tobias Westmeier [twnature]* (iNaturalist). *Reproduced under creative commons license CC-BY-NC

Month of observation was converted to radians (θ = month × 2π / 12) to treat time as angular data. For each species and region, the following metrics were calculated: (1) mean direction (μ), representing the average peak month; (2) concentration parameter (r), measuring the strength of seasonality (0 = uniform distribution [no seasonality], 1 = all observations in one month [very high seasonality]); and (3) Rayleigh’s test for non-uniformity, which tests the null hypothesis of a uniform circular distribution (Rayleigh 1919). Watson’s two-sample U² test was used to compare seasonal distributions between regions and time periods. This non-parametric test evaluates whether two circular samples are drawn from the same underlying distribution and is robust to deviations from the von Mises distribution (Jammalamadaka & SenGupta 2001).

DATASET CHARACTERISTICS

The retrieved dataset comprised 6,629 research-grade observations spanning 15 years [2010–2024], from 15 species. Observations were strongly biased toward eastern Australia (n = 5,508, 83%) compared to the west coast (n = 1,033, 16%). The most frequently observed taxa were A. argus (n = 1,829), D. auricularia (n = 886), and A. juliana (n = 710). Observations increased strongly from 2010 (n = 13) to 2024 (n = 1,456), reflecting rapid growth in iNaturalist participation in Australia.

Objective 1 – Seasonality across the eastern and western Australian coasts

When observations were aggregated by coast, east coast data showed highly significant seasonality (Rayleigh’s Z = 529.2, N = 5,508, p < 0.0001) (Appendix: Table A1) with a clear summer peak in December (μ = 0.01 radians, concentration r = 0.31). In contrast, west coast observations exhibited weak seasonality (Rayleigh’s Z = 1.7, N = 1,033, p < 0.001; r = 0.04) with observations distributed almost uniformly throughout the year, peaking nominally in September but with maxima also in April and January. Watson’s two-sample test confirmed that east and west coast exhibited significant seasonal distributions (N¹ = 5,508; N² = 1,033; U² = 0.89, p < 0.001) (Appendix: Table A2a). However, analysis by latitudinal zone (Fig. 16) revealed that this apparent coast-level difference largely reflects latitudinal composition of observations on each coast, with tropical populations peaking in August–September and temperate populations in December–January regardless of coast longitude (Fig. 16; Appendix: Table A3).

Fig. 16

Heatmap of observation densities (month × latitude) for data retrieved from iNaturalist

Objective 2 – Species-specific differences in phenology (temporal niche partitioning)

Species-level analyses initially revealed substantial variation in seasonal concentration when data were pooled across all locations (Table 2). However, examination of geographic structure within species revealed that apparent weak seasonality patterns in several species was an artifact of aggregating observations across tropical, subtropical, and temperate regions.

Table 2

Species-specific differences in phenology (pooled across latitudinal zones) of Australian sea hares. N = number of observations, μ = mean direction (peak month in radians), r = concentration parameter [closeness of fit to circular model] (thus 1 = all obs in 1 month, 0 = no seasonality), Rayleigh p = significance statistic from test of uniformity (evaluates whether observations are distributed uniformly around the year (null hypothesis) or show significant clustering at particular times). Test statistic Z = n × r². * Species marked with an asterisk show significant geographic structure (see Table 3)

When analysed by latitudinal zone (Tropical <23°S, Subtropical 23° to 29°S, Temperate >29°S), most species exhibited stronger patterns of regional seasonality than when pooled (Table 3). For example, A. oculifera, which exhibited a non-significant year-round pattern when data were pooled (r = 0.14, Z = 1.2, p = 0.45, N = 62), showed some significant seasonal patterns when analysed by region: temperate observations exhibited a significant January concentration (r = 0.65, Z = 8.0, p < 0.001, N = 19 interpreted with caution given the small sample size) while tropical observations showed weaker non-significant peaks in August (r = 0.20, Z = 1.7, p = 0.19, N = 42). Similarly, A. argus, exhibited highly significant but weak overall seasonality (r = 0.18, Z = 59.3, p < 0.001, N = 1,829) but when analysed per region, patterns showed highly significant tropical observations peaking in October (r = 0.41, Z = 39.9, p < 0.001, N = 237) and again for temperate observations in January (r = 0.31, Z = 102.0, p < 0.001, N = 1,062).

Table 3

Species-specific differences by latitudinal zone (revealing geographic structure) for taxa that exhibit mixed latitudinal phenology. Division into latitudinal zones (Tropical: <23°S; Subtropical: 23° to 29°S; Temperate: >29°S) shows patterns that were obscured in pooled data. Only zones with ≥10 observations shown. N = number of observations, r = concentration parameter [closeness of fit to circular model] (thus 1 = all obs in 1 month, 0 = no seasonality). Test statistic Z = n × r²

Species showing consistent phenology across regions included B. hirsuta (January peak, r = 0.90, all observations temperate), A. extraordinaria (January peak, r = 0.64, predominantly temperate), and S. striatus (January–February peaks across zones, though with some regional variation in timing). In contrast, several widespread species showed markedly different peak timing across latitudes. A. concava peaked in July–August in tropical regions (r = 0.56) but in December in temperate zones (r = 0.50), a difference of approximately 5 months. D. auricularia showed June peaks in the tropics (r = 0.44) compared to January in temperate areas (r = 0.37).

The apparent “winter peak” of A. reticulata (August, pooled data) was driven primarily by tropical observations (n = 228, August peak, r = 0.57), while temperate populations showed autumn peaks (March, n = 46, r = 0.53). In the tropics, August represents the late dry season rather than winter, highlighting the importance of interpreting phenological patterns in their appropriate regional context.

Within-year periodicity confirming seasonality (Objectives 1–2)

Spectral analysis of detrended data confirmed dominant 12-month periodicity for temperate latitudinal bands, consistent with the unimodal seasonal pattern evident in Fig. 17 (Appendix: Table A4). When data were pooled across latitudes for the west coast, a six-month (semi-annual) periodicity appeared, consistent with the multimodal pattern that showed local peaks in September, April, and January. However, this signal represents a composite of offset annual peaks from tropical and temperate observations rather than a genuine six-month biological cycle, consistent with the latitudinal structure evident in both the species-level (Table 3) and assemblage-level (Fig. 17) analyses. Species-specific spectral analyses supported this interpretation: highly seasonal species (e.g., A. juliana, and D. auricularia) showed strong 12-month peaks regardless of whether observations were made on the east or west coast (Appendix: Table A4).

Fig. 17

Species specific seasonal patterns for selected Australian sea hares

LATITUDINAL DRIVERS OF PHENOLOGY

The most prominent pattern in sea hare phenology was the presence of latitudinal gradients in peak timing, with temperate populations (south of 29°S) peaking in December–January and tropical populations (north of 23°S) peaking in August–September (Fig. 18). This gradient is likely to reflect fundamental differences in environmental seasonality across latitude. In temperate Australian waters, sea surface temperatures can fluctuate by as much as 4–6 °C between summer and winter (Ridgway & Dunn 2003, Schaeffer & Roughan 2017), creating a strong environmental cue for synchronised reproduction in ectothermic marine invertebrates. Sea surface temperatures above 18–20 °C are generally required for successful spawning and larval development in temperate Aplysia species (Kandel 1979, Susswein et al. 1983), and the timing of thermal thresholds is likely to constrain population-level phenology to the summer months. However, links between temperature and phenology are inferred from the coincidence between observed peaks and known thermal optima rather than directly tested against sea surface temperature data, thus the mechanistic interpretation remains speculative.

Fig. 18

Monthly proportion of total aplysiid observations by latitudinal bands. Colours represent four degree latitudinal bands, line style indicates zone type (tropical, subtropical & temperate). To better visualise peaks the x-axis is summer centred

In tropical regions, there is relatively little annual sea surface temperature variation, but macroalgal communities nonetheless show pronounced seasonality. On Western Australian tropical reefs, where most Australian A. reticulata records originate, canopy-forming Sargassum spp. reach peak biomass in summer and senesce through winter, while understory brown algae such as Dictyopteris and Lobophora reach peak biomass in winter when the Sargassum canopy is at its minimum (Fulton et al. 2014). The August–September peaks in tropical sea hare observations coincide with this winter understory-algae peak, providing a plausible food-availability mechanism that is consistent with the algivory typical of Aplysiidae (Nimbs et al. 2017a). These contrasting drivers – temperature-driven seasonality in temperate waters versus macroalgal community turnover in the tropics are likely to produce the latitudinal gradient evident in Fig. 17.

The apparent east–west coast differences in phenology (Appendix: Table A1) are a product of the latitudinal gradient interacting with geographic distribution of observations. East coast data are dominated by observations from temperate New South Wales and Victoria, producing a strong December peak whereas data for the west coast spans a narrower latitudinal range, and mixed latitudinal observations with offset peak timing has acted to produce a weak seasonality signal and multimodal pattern. Oceanographic differences between coasts – the strongly seasonal East Australian Current versus the thermally stable Leeuwin Current (Feng et al. 2003, Wernberg et al. 2012) – are likely to contribute to the patterns, but latitudinal sample composition is sufficient to explain most of the observed contrast.

SPECIES-SPECIFIC PATTERNS AND ECOLOGICAL IMPLICATIONS

The gradient from highly seasonal (r > 0.8) to year-round (r < 0.2) species suggests temporal niche partitioning within Australian sea hare assemblages. Tightly seasonal species like B. hirsuta (Nimbs & Wilson 2020, Wells et al. 2021) and A. gigantea may be specialists with narrow thermal or dietary requirements, restricting their activity to optimal summer conditions. In contrast, taxa with weak seasonality (A. oculifera and B. leachii) may be thermal generalists capable of exploiting resources across a broader temperature range or may comprise populations with staggered cohorts producing year-round reproductive output.

The August peak in observations of A. reticulata is notable, as it suggests this species may occupy a distinct temporal or latitudinal niche. Winter-peaking marine invertebrates are relatively uncommon in Australian waters, and this pattern may reflect specialised food requirements coinciding with particular macroalgal communities (Pennings 1990). Alternatively, the apparent winter peak in this taxon may have resulted from higher detectability during cooler months if this species is cryptic or less active during summer.

The clustering of peak months among common species (November–December) raises questions about potential competitive interactions. If multiple sea hare species compete for similar macroalgal resources, such temporal overlap could intensify resource competition during population blooms. However, sea hares exhibit diverse dietary preferences (Pennings 1990, Rogers et al. 1995, Nimbs et al. 2017a), and spatial segregation by habitat or depth would reduce direct competition even when phenological patterns overlap.

GEOGRAPHIC STRUCTURE IN SPECIES-LEVEL PATTERNS

Pooled species-level analyses can underestimate seasonal concentration when populations exhibit different phenologies across latitudinal gradients. Several species classified as weakly seasonal based on pooled data (A. oculifera, A. argus, D. auricularia) showed moderate to strong seasonality within geographic regions (Table 3). As with the east and west coast comparative analysis, this pattern arose because tropical and temperate observations often peak at different times of year – sometimes differing by 5–7 months – causing their signals to partially cancel when combined.

For example, A. oculifera appeared essentially year-round in pooled analyses (r = 0.14, non-significant) but exhibited a January concentration in temperate regions (r = 0.65). Similarly, A. concava showed opposite seasonal peaks in tropical (July–August, late dry season/early wet season) versus temperate (December, summer) regions, likely reflecting responses to different temperature and productivity regimes.

These patterns highlight the need to consider latitudinal range when characterising phenology in widespread species, as pooling observations across large geographic areas may generate assemblage-level patterns that do not represent local population dynamics. The regional structure also suggests that sea hare phenology is coupled to local environmental conditions rather than being determined solely by phylogenetic constraints or life history traits.

METHODOLOGICAL CONSIDERATIONS AND THE IMPORTANCE OF DETRENDING

Exponential growth in iNaturalist participation has been documented globally (Callaghan et al. 2019), and failure to account for this trend can produce spurious periodicities through the interaction of trend and seasonal components (Grolemund & Wickham 2011). The pronounced increase in observation numbers from 2010 (n = 13) to 2024 (n = 1,456) introduces challenges for detecting phenological shifts. Additionally, the youth of iNaturalist means that the 15-year reliable data period is relatively short for climate-driven phenological studies, which typically require 30+ years of data to distinguish directional trends from natural variability (Edwards & Richardson 2004, Poloczanska et al. 2013). Strong increases in sample sizes over recent years may introduce bias when using trend analyses if observer behaviour or coverage changes systematically over time (Courter et al. 2013).

An additional concern is whether seasonal patterns might reflect within-year variation in observer activity rather than actual sea hare biology. Observer effort and sea hare seasonality were decoupled: all-Australia iNaturalist uploads peaked in September–October with a June minimum (Appendix: Table A7), whereas temperate Aplysiidae observations peaked in January (18.0% of records vs 8.2% baseline) and tropical sea hares in August–September (up to 22.3% vs 12.0% baseline). The mismatch indicates that the seasonal patterns reported here reflect biology rather than observer effort bias. Depth of observation is not recorded by iNaturalist and cannot be reliably inferred from coordinates given the variable bathymetry of nearshore Australian waters; the dataset therefore aggregates intertidal and shallow subtidal records without a recoverable per-observation depth, which is treated here as a component of broader observer-effort variation.

Despite these limitations, citizen science data offer unique advantages for phenological research. The geographic and temporal coverage achievable through platforms like iNaturalist far exceeds what is feasible through traditional field surveys (Dickinson et al. 2012, Pocock et al. 2015). For conspicuous, taxonomically-stable species like sea hares, research-grade identifications are generally reliable, particularly at the genus or family level. Future studies could strengthen phenological inference by incorporating environmental covariates (e.g., satellite-derived sea surface temperature) to normalise observations and account for spatiotemporal variation in observer effort (Szabo et al. 2010, Callaghan et al. 2019).

EVIDENCE FOR DIRECTIONAL SHIFTS

Australian coastal waters have warmed significantly over recent decades, with the eastern region showing particularly rapid temperature increases (Ridgway 2007, Wu et al. 2012). Circular regression of observation timing against year provided mixed evidence for directional shifts: two of eight tested species exhibited statistically significant directional shifts over 2010–2024. A. sydneyensis, a well-documented temperate NSW endemic, shifted earlier by approximately one month per decade – a direction consistent with regional warming and broadly comparable in magnitude to shifts reported for other marine taxa (Poloczanska et al. 2013). D. auricularia, with its largest populations in the tropics, shifted later by approximately 1.5 months per decade, potentially reflecting altered monsoon dynamics rather than a simple temperature response. The remaining six species showed no detectable directional shifts. Both significant effects attenuated to non-significance under Bonferroni correction in a 2015–2024 sensitivity window, indicating that the full-period results are borderline for statistical power. This mixed result suggests three non-exclusive interpretations: (1) phenological responses to warming are species-specific and depend on local environmental drivers (temperature for temperate species, monsoon dynamics for tropical species); (2) the observation period (15 years) is at the lower end of what is typically required to detect climate-driven phenological signals (Edwards & Richardson 2004, Poloczanska et al. 2013); and (3) other environmental factors (photoperiod, food availability) may constrain phenological responses regardless of temperature. No quantitative historical phenological data are available for Australian sea hares against which to benchmark the present results; previous Australian work has described boom–bust population dynamics linked to food availability (Nimbs et al. 2017a) without resolving the seasonal structure documented here.

The baseline phenological data established here provide a foundation for detecting future shifts as warming continues. Continued monitoring through citizen science platforms, ideally integrated with oceanographic time series data, would enable more robust detection of climate-driven changes. Particular attention should be paid to species at the margins of their thermal ranges, such as tropical species reaching their southern limits along the east coast, i.e. the rarely-observed A. reticulata, which may show earlier phenological responses to warming.

The latitudinal gradient in phenology raises questions about how populations at different latitudes may respond to future environmental change. Temperate populations, with tightly synchronised seasonal reproduction, may be more vulnerable to phenological mismatches if warming decouples temperature cues from food availability (Edwards & Richardson 2004, Winder & Schindler 2004). Tropical populations, with phenology linked to wet/dry season dynamics rather than temperature alone, may respond differently to warming. Subtropical populations at the interface between these regimes warrant particular attention as potential early indicators of climate-driven phenological shifts.