random slopes by two groups (species and site)?

I have a dataset of 20 bird species counted at 8 sites, and I'd like to estimate the species-specific change in migration timing at each site. Roughly half of those species are recorded at >= 5 sites; 7-15 species are recorded depending on the site (gray tiles indicate sites where a species is not counted or there is insufficient data):

If I was just focused on a single site, I think it would be sensible to fit a random slopes model (partially pooling by species) to obtain the species-specific change in migration timing (days per year) at that particular site. In this model, 'doy' is the day-of-year, and year is 1985, 1986, etc.:

m1 <- lmer(doy ~ year + (1 + year|sppID), data = df) 

Now I'm wondering whether I could 'scale-up' this approach to also incorporate variation in siteID? If I run a linear model for each species x site combination, I not only see that different species are changing their migration timing in different directions, but that a single species can exhibit advancements/delays/no change in timing depending on the site.

My question is whether there's a coherent way to fit all the data with a single hierarchical model to obtain species-specific changes in timing that can also vary by site, and what the correct model specification would look like? I'm assuming too that some species would not be counted across enough sites to include in a model like this? It would be ok with me to drop those rare species in the dataset.

Thanks,

Jay