It is very rare to find microbiome studies where variance inflation and random forest algorithms are implemented. This is a good approach!
However, the community composition of the control treatment EFB-0 is quite different from what has been reported in normal honeybee workers: there is an overabundance of Gilliamella apicola and Bartonella bacilliformis (which I think is actually Bartonella apicola ), and the abundance of the Lactobacilli is remarkably low (~30%). This clade is well known to dominate the hindgut. I could not find this discussed in the manuscript.
This leads me to the following question: are the relative abundances normalized? this should be done to compare between samples, since a difference in the amount of sample from where DNA was extracted can have an impact in the sequencing results. This should also be normalized for 16SrRNA copy number.
The normalization step could be the cause of the variation attributable to Julian days (each sample would be different on its own despite the treatment). And more importantly, since OTU relative abundance is Compositional Data, it would also affect the results on the effect on core-gut abundances: an increase in Melissococcus is necessarily followed by a decrease in everything else.
This should be addressed in a final version, since the relative abundance would affect all other downstream statistical analyses.