Estimate Difference in percent data and numeric explanatory variable

I am trying to estimate the impact of temperature on the change in coral cover at different sites. So far, I estimated simultaneously the cover before and after a coral bleaching event for different species:

brm(bf(mvbind(before, after) ~ species,
    phi ~ (1|species)),
    family = Beta(),
    data = data)

Then I use predicted_draws() to calculate the difference of the posterior distributions for before and after.

Data exists for sites with different temperatures, I’d like to do a regression (with maximum temperature as a dependent variable) to directly estimate the change in cover. How could this be done?
I chose the above model because at each site, replicate cover measurements exist but the exact location of the replicate differed between the years. I want to include the variability between the replicates and not use the difference of the means.

Thanks a lot for your help!

g’day @aeich ,
I suppose one place to start would be to condition cover on temperature like this:

brm(bf(mvbind(before, after) ~ species + temperature,
    phi ~ (1|species)),
    family = Beta(),
    data = data)

But if you have the data to support it, maybe you could include a (species|temperature) varying effect as well.
I notice you include (1|species) for phi but not for the mean response. Have you also considered varying effects by species for the conditional mean?