The authors did an excellent job in addressing my comments. I have nothing else to add.

The revised version of this manuscript addresses all points I made on the previous draft.

Again, I think this is a very good paper and would happily accept it for publication if I were in a position to do so.

Well done!

This is a brilliant piece of science: well written, carefully conducted, based on an outstanding dataset, using a thorough and sophisticated methodology, and presenting timely and very exciting results. I really enjoyed reading it! I have only few very minor questions and suggestions.

I could only fully understand the abstract after reading the completed manuscript. I would suggest rewording small details so that the abstract is clearer. Some points where I found difficulties:

“Community structure was mainly determined by…” Here, I found “community structure” to be somehow vague. Moreover, in a previous sentence you mentioned that network inference was used. After reading this, I was expecting you characterized the community structure via network properties, which is not the case. This may nevertheless be a matter of personal bias, but I guess others may have the same problem. I would suggest to change “community structure” by “species abundances”.

“The relative importance of these factors was mildly modulated by host plants.” This sentence was also difficult to follow to me without reading the full manuscript. I would say something like: “The relative importance of these factors mildly varied when we used particular host plant groups” o something alike. This may be again a matter of taste.

“… specialists and generalists flies almost behaved as separate communities…” I found “behaved as separate communities” difficult to understand here and when mentioned throughout the text. I would try to use a term more clearly connected with the results. Again a matter of taste and totally up to the authors to follow this suggestion.

In the second paragraph of the introduction, it may be worth mentioning that facilitation can also occur between phytophagous arthropods (e.g. Godinho et al. 2016. Oecologia 180: 161-167.)

“Since species interactions mostly occur in/on plant organs, they may be modulated by plant species identity…” Here, it is not totally clear whether “species interactions” refers to intraguild interactions or to fly-plant interactions.

Are the 8 species used all the species of Tephritidae present in the island? If so, I would explicitly state it. If not, I think it would be worth mentioning in the discussion the potential influence of other unevaluated species in the abundance of the used species and model outputs.

“Of the 12872 initial samples, only those with GPS coordinates, with at least one individual fly and belonging to one of the 21 host…” As a layman in the modelling used, it is not clear to me why you didn’t use the 0s (i.e. the samples without individuals). Moreover, how many host plants are there in the island? If there are much more than the 21 used, how do you think this could affect subsequent interpretations of assembly mechanisms? This might deserve a line in the discussion.

“… all previous models were reevaluated on the datasets excluding D. ciliatus (Models 2-0 to 2-6)” To facilitate the reading I would say: “… all previous models were reevaluated on the datasets excluding the species lacking fundamental host use estimates (D. ciliatus ; Models 2-0 to 2-6)” or something similar.

“Among plants inferred as possible hosts from species abundance patterns (i.e., those with high coefficient values), coefficients correlated positively with fly laboratory-measured fitness for specialists but not for generalists (Figure 3B)” I found this result super interesting! Still, I would better justify why you only used the plants with high coefficients. Moreover, while the modelling approach used seems evident, I would explicitly specify it (perhaps in the figure caption), explaining also the meaning of the shadowed areas in Fig. 3B. Finally, could this result be mostly driven by only one generalist and/or specialist species?

“Accounting for environmental covariates strongly improved model fit and made all residual covariances almost completely vanish, particularly among groups, suggesting that no important environmental factor structuring the community has been missed”. Completely agree. Yet, your data is temporally and spatially structured, which might contain interesting information on assembly mechanisms (e.g. the effects of dispersal processes and interannual and/or seasonal dynamics not linked to climate). While testing this may be out of the scope of your work, I also think that it would deserve a line in the discussion. At least, to avoid criticisms from autocorrelation purists, I would explicitly mention that your results suggest that the influence of temporal and spatial autocorrelation is negligible (besides the influence of other important environmental factors).

While I really enjoyed “The ghost of competition past” section in the discussion, I’m not totally sure the term perfectly fits here, at least as described by Connell (1980). To my knowledge Connell was referring to the coevolution of competitors and thus to evolutionary changes in the fundamental niches. That is, I agree that by nicely comparing realized and fundamental niches you detected a “ghost of competition”, but I would rather say that it is a current (ecological) ghost rather than a past (evolutionary) one. Perhaps it would be worth to add a few lines explaining this, mentioning also that eco-evolutionary approaches are required to truly address changes in fundamental niches due to competition and their consequences to the assembly of species.

Finally, regarding also the fundamental niches, by looking at figures 3A and 1B, it seems that in some situations flies are able to colonize plant species in which they show a very reduced fitness (close to 0, if not 0). I kept wondering how it is this possible. Perhaps, there are differences in host use among (sub)populations and fundamental niche estimates are based on individuals of a reduced number of (sub)populations. I’m just not sure, but this would perhaps deserve a brief mention in the discussion.

Hope this is of any help and congratulations for this excellent study!

Joaquín Calatayud

Facon et al.: Joint species distributions reveal the combined effects of host plants, …

This study investigates fruit fly communities on different host fruits on la Réunion over an 18 year period. Abundances of fruit flies were predicted using both climatic variables and host plant species in a joint species distribution model. This approach yields correlations among fruit fly species as a side-effect (the variance-covariance matrix of the residuals of the k x k species, also sometimes called the association matrix), which may represent species interactions.

Accounting for environment, the association matrix became almost diagonal, indicating no apparent associations among fruit fly species. This is in contrast to the raw-data association matrix, which showed a strong difference in host use between generalists and specialists. The distinction between generalists and specialists was based on a fly x host laboratory choice experiment (published elsewhere by members of the same group). In the discussion, the authors carefully interpret the few remaining (all negative) correlations in the residuals as indication of competition among some generalists and some specialists. In particular, they expand on the problem of not seeing many competitive interactions due to evolution having led to niche displacement among similar species (aka the ghost of competition past).

Overall, I find this one of the best jSDM studies I have seen in the literature so far. The relatively small number of species (only fruit fly species) and the huge number of observations (5000 samples with nearly 100,000 individuals) is hard to improve upon. Of course, there are always a few things that remain unclear (at least to me). For example, competition is likely to leave an effect only when resources are scarce, which requires either high population densities of the flies, and/or low availability of fruits. The authors do not provide data on either. Thus, for any given sample, it is unclear whether we would actually expect any trace of competition.

Also I am not very impressed by the model selection approach used, although admittedly this is almost a philosophical issue and their practice is in line with common analytical strategies (more on this below).

I am particularly happy to see the substantial effect of relatively coarse environmental predictors on the association matrix. It could be argued that the small remaining covariance could actually be explained by other predictors, such as acidity of the fruit or something like that. If so, no covariance would indicate no competition, expounding the problem of witnessing the outcome of hundreds of generations undergoing niche separation without being able to see the selection in action. Only extensive laboratory studies with monospecific and paired flies over generations could test for a widening and hence overlapping host niches when released from the invisible competition.

As far as observation studies go, I don’t think there is anything more we could ask for.

Model selection bias

I wouldn’t quite call it a secret, but it seems that the well-known effect of model selection causing a bias in model parameter estimates is unknown to ecologists. If, however, we are interested in model parameters, not only their prediction, we need to be aware of this. In a nutshell, the selection of models leads to a final (few) “best models”, which are then interpreted statistically (“significant effects”) and their parameters are estimated alongside their standard error. Now, since the computer does not “know” that this best model is the result of investigating dozens to hundreds of models on the same data before, it naively assumes this is the only model fitted to the data. It “ignores” the model selection uncertainty that results from variability of the data. As a result, the standard errors are too small, as they do not take into account the model selection procedure’s uncertainty. Also, the estimates are biased, as all models that show non-significant effects of this predictor are removed through model selection, making the remaining models more likely to have large (absolute) estimates than small.

It is relatively simple to show this through simulation ([http://joshualoftus.com/post/model-selection-bias-invalidates-significance-tests]), but Harrell (2001, Springer: Regression Modeling Strategies) writes about it, and it is found (as obvious introductory statement) in [https://pubmed.ncbi.nlm.nih.gov/18189162] or “well known” in [https://academic.oup.com/aje/article/167/5/523/212182]. Parameter shrinkage has been advocated as strategy to counteract the selection bias (see last link for review).

Now, in the present study I do not actually think that model selection bias is a big issue: the data set is large (reducing the problem of model selection leading to different models during bootstrapping); also, the BIC difference between models is very large, indicating little ambiguity in the model ranking. Still, I would have preferred a presentation of few fuller models than the model selection outcome (as indeed Burnham & Anderson themselves argue for).

Model diagnostics

I did not find any statements on how well the model structure meets distributional assumptions and independence of residuals. While I guess that the PLN approach is relatively robust to overdispersion, I think it can be expected of the authors to provide a statement on whether the data were actually well modelled assuming a Poisson distribution. My experience is that fast reproducing and flying beasts, such as fruit flies, tend to clump, requiring a negative binomial to represent the variance in the data. (The DHARMa package may not readily work on this model type, but the principle applies and the authors should be able to simulate the residuals themselves, as demonstrated in that package’s vignette.)

Similarly, the spatial nature of the sampling warrants an assessment of spatial autocorrelation in the model residuals. In fact, it could be that the remaining association signal can be partly explained by spatial effects.