Statistical analysis
Each sound produced by a different O. croaticus male was assumed
to be an independent acoustic unit, and the statistical analyses were
performed by combining the sounds from multiple individuals into a
single dataset. In the intraspecific analyses, the O. croaticusindividuals were utilised as a grouping variable, to explore acoustic
variation between males. For each spectral and temporal variable of the
sound, the descriptive statistics are presented. Outliers and extremes
were detected visually from the boxplot and were eliminated from the
dataset if necessary. In order to test the assumption of normality, we
initially performed Shapiro-Wilk normality test on a raw intraspecific
dataset. Since the assumption of normality was not met for some
variables, the overall acoustic dataset was then transformed (either
using log or square root functions) followed by a Box-cox function to
estimate the transformation parameter by maximum likelihood estimation.
The acoustic dataset was re-examined for distribution fitting, and since
assumptions of normality (Shapiro - Wilk test, P <
0.05) and homoscedasticity (Bartlett test, P < 0.05) of
the variances were not achieved, we continued with non-parametric tests.
For pairwise comparisons between soniferous O. croaticus males,
we employed the non-parametric Kruskal-Wallis rank sum test Hfollowed by pairwise Dunn’s multiple comparison test with Bonferroni
correction for the P values. To investigate the mutual
relationship between mean individual acoustic variables, we utilised
non-parametric Spearman correlation tests. The association
between acoustic variables with body characteristics
(LS, W and Fulton’s K ) and water
temperature (T, in °C) was performed with Spearman correlations.
Additionally, the Chi -square (χ2) was used to
test for independence of behaviour (expressed as behavioural categories)
from sound production. In this test, the residuals from the
χ2 were used to determine which behaviours were
positively related to sound production. Kruskal-Wallis H test was
used to compare the mean behavioural variables (calling rate, behaviour
rate, n. of female nest entrances) between soniferous males. Wilcoxon
signed-rank test was performed to compare the two dependent samples,i.e., mean behavioural variables (behaviour rate and female nest
entrance) of males when they produced sounds and when they were mute.
Additionally, Wilcoxon test was used to compare the frequency and
duration of courtship and pre-spawning phases between males.
For the interspecific comparisons, the means of individual acoustic
properties of soniferous sand gobies were compared with the
non-parametric Kruskal-Wallis H test. To quantify interspecific
acoustic variability among the soniferous sand gobies (generaKnipowitschia , Orsinigobius , Pomatoschistus andNinnigobius ), we used a multivariate approach. Principal
component analysis (PCA) is a commonly used method in the bioacoustics
literature for detecting the variables that explain the most variance
among soniferous fish species or populations. PCA, based on the
correlation matrix, was performed on transformed and standardised
individual means of five sound variables (temporal: DUR, NP, PRR;
spectral: PF and FM) to assess overall acoustic variability between sand
gobies, and additionally to recognise acoustic variables explaining the
observed variance. To assess the percentage of successful classification
of the sounds assigned to the correct goby species, and to maximise the
separability among taxa, we used linear discriminant analysis (LDA). Two
different LDAs were performed, first with the complete dataset (five
acoustic variables for each species) and then removing the
temperature-dependent features (DUR and PRR). Our results were presented
as means (x̄ ) ± standard deviation (s.d. ), while the level
of significance for inter- and intraspecific comparisons was 5% (α =
0.05). Statistical analyses were performed in
STATISTICA® (v. 13.6.0., TIBCO, USA), Past (v. 4.11)
and R Studio (2022.07.0) software.