I think the below would work? If you need to join on multiple character columns that exist in both, can do by = c("var1", "var2", ...)

  result_df <-  df2 %>%  group_by(species) %>% 
  summarize(mean_variable1 = mean(variable1)
            ,mean_variable2 = mean(variable2)) %>% 
  inner_join(df1, by = c("species"))