In reading the manuscript I was surprised to see the result that abundance had declined by 21% and species richness had remained reasonably constant. This seems at odds with the observation than in both plots we see an increase in the maximum abundance and richness of samples over time with no obvious change in the minima. This suggests that the shape of the distribution is changing through time.

It would seem that the method used to analyse these data may not be appropriate; the linear model assumes that the model errors are independent and identically distributed mean zero, Gaussian random variables with the same variance \hat{\sigma}^2. This implies that conditional upon the predictor variables, Y, the response, also has those distribution properties. In both cases, abundance and richness would not normally meet these assumptions; abundance and richness can't be negative and in general we see an increase in variance with increasing mean where the linear model assumes the mean and the variance are unrelated.

That there are obvious changes occurring in parts of the distribution of the of the abundance and to a lesser degree the richness data, it is important to unpack the analysis you performed and ask whether this is a real decline in the mean but a change in the variance, or whether the method you used, in potentially failing to account for the distributional properties of the data, correctly estimated the mean response. Are the results robust if you fit a suitable generalized linear model with a discrete count distribution as the error distribution? A starting point would be the Poisson GLM, with an offset term for area in the case of Abundance (where you model the total count, not the count per unit area, and then include the log(area) as an offset term), but I suspect overdispersion may be an issue and you may need to look beyond the Poisson to say the Negative Binomial.

It is not clear from the description of the data if the same sampling location appears on more than one occasion in these time series or are all points independent? If the same sites were sampled on more than one occasion, this should also be accounted for in the model. You do not state that any controlling for autocorrelation was performed. Both of these issues complicate the analysis; if sufficient sites were sampled on more than one occasion a random effect for site and time (slope) would allow you to account for the clustering of observations at the site level and for potentially different site slopes/trends. This would need a GLMM or other hierarchical model. Adding temporal autocorrelation complicates matters further, but is something to consider.

In contrast, if the mean is shown to be constant over time, then has there been a trend in variance? If so, what would account for this? Is this a sign that perhaps in the most abundant sites those sites have seen an increase in abundance but at lower abundance sites there has been a constant or even decline in total abundance?

In Figs 2c and d you show plots of the change in % abundance (composition) of 4 groups of taxa, again with trend lines derived from a linear model. Proportional data are constrained to lie on the interval 0-1 and the variance is not constant over this interval. There also seems to be some strong autocorrelation in some of the trends, which it is not clear was accounted for. If you wish to model the % of the the individual groups separately then a beta regression would be a starting model. If you want to model the changing proportions of all groups at the same time, then a dirichlet regression would be an appropriate model choice. These would account for the assumed distributions of the data.

Also, related to this, you state (bottom p 11, top p 12)

This loss was mainly driven by a decline in abundance of low-temperature dwelling species (slope ± SE = -0.60 ± 0.23, F1,22 = 6.60, P = 0.018) that is only partly mitigated by the concurrent increase in high-temperature dwelling species (slope ± SE = 0.05 ± 0.02, F1,22 = 10.46, P = 0.004), which contribute less than 5% to the total benthic invertebrate abundance (Fig. 2c)

This relates to the loss of abundance in 2a, which you attribute to the loss of low-temperature species and the lack of concomitant response in the warm-loving species, which you refer the reader to 2c to observed. Are the statistics quoted for the regression lines in 2c? If not what are they from? If they are, these do not provide evidence in support of the explanation for the loss of total abundance. If the total abundance had remained constant throughout you could still see a shift in relative composition. From the plot is looks like, as you'd expect given that the data are closed compositional values, the decline is cold-dwellers is offset by smaller increase but in all three other groups. These results speak to a change in composition over time, with a replacement of cold-dwellers by somewhat warmer-dwelling taxa, but on their own they do not account for the suggested decline in mean total abundance over time.