See my answer on compositional response to https://stats.stackexchange.com/a/638757/284766 that includes a discussion on the book "Analyzing Compositional Data with R." Section 5.3 "Compositions as Dependent Variables" shows how to model the relative abundance of each individual species in relation to its traits. The method transforms the response variable of several columns of proportions using an isometric log-ratio transformation (ilr) into orthonormal coordinates before employing OLS linear regression. Because not all species are present in all locations, there will be some zero proportions. This will challenge the model fitting, as the transformation is based on log ratios and neither log(.../0) nor log(0/...) is defined. Chapter 7 discusses how to deal with zeroes. Also because of zero proportions, Dirichlet regression, not discussed in the book, is not directly usable. One method to circumvent zero proportions is to add a tiny amount to this proportion. See https://maartenbuis.nl/presentations/berlin10.pdf.

Because of zero proportions, the best choice might be fractional multinomial logit models. I do not think that you should use models of random or mixed effects, unless you feel that proportions measured at different locations are correlated (spatial correlation) or multiple measurements of all species are taken at each location (longitudinal data). Although each species are measured more than 250 times, but these observations are independent unless there are spatial correlation among locations or temporal correlation among measurements. Species-specific intercepts in a multinomial logit model will assess the mean unobserved deviation by species.

See my answer on compositional response to https://stats.stackexchange.com/a/638757/284766 that includes a discussion on the book "Analyzing Compositional Data with R." Section 5.3 "Compositions as Dependent Variables" shows how to model the relative abundance of each individual species in relation to its traits. The method transforms the response variable of several columns of proportions using an isometric log-ratio transformation (ilr) into orthonormal coordinates before employing OLS linear regression. Because not all species are present in all locations, there will be some zero proportions. This will challenge the model fitting, as the transformation is based on log ratios and neither log(.../0) nor log(0/...) is defined. Chapter 7 discusses how to deal with zeroes. Also because of zero proportions, Dirichlet regression, not discussed in the book, is not directly usable. One method to circumvent zero proportions is to add a tiny amount to this proportion. See https://maartenbuis.nl/presentations/berlin10.pdf.

Because of zero proportions, the best choice might be fractional multinomial logit models using fmlogit as mentioned in my answer. I do not think that you should use models of random or mixed effects, unless you feel that proportions measured at different locations are correlated (spatial correlation) or multiple measurements of all species are taken at each location (longitudinal data). Although each species are measured more than 250 times, but these observations are independent unless there are spatial correlation among locations or temporal correlation among measurements. Species-specific intercepts in a multinomial logit model will assess the mean unobserved deviation by species. If random intercepts and slopes are somehow needed, you can check out the package mlogit, but I do not think this package accepts fractional response, so you may have to crack the original script.

If you concern that some species are similar and share unobserved traits, you could consider nested logit models. See https://cran.r-project.org/web/packages/mlogit/vignettes/e2nlogit.html.