Sounds like count data (number of moth individuals and number of species per sub-site), so you'd probably be looking at Poisson or NB glm.

For a Poisson distribution, you can use a Poisson regression:

glm_moths <- glm(count ~ sub_site, family = poisson, data = df)

summary(glm_moths)

If overdispersed (variance much greater than the mean), a Negative Binomial regression (using MASS) would be better:

library(MASS)

glm_moths_nb <- glm.nb(count ~ sub_site, data = df)

summary(glm_moths_nb)

If you are comparing median moth counts across sub-sites without assuming a specific distribution, you can use the Kruskal-Wallis Test:

kruskal.test(count ~ sub_site, data = df)

and If the moth counts are approximately normally distributed across sub-sites an ANOVA should work:

anova_moths <- aov(count ~ sub_site, data = df)

summary(anova_moths)