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 ( ) ± 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.