library(tidyverse)
# calculate the standard deviation for the two fields
sd_N <- sd(data$N)
sd_C <- sd(data$C)

# This will get us the plot for field `N`
data %>%
  group_by(Year) %>%
  summarise(mean_N = mean(N)) %>%
  ggplot(aes(x = Year, y = mean_N)) +
  geom_point() +
  geom_line() +
  geom_errorbar(aes(ymin = mean_N - sd_N, ymax = mean_N + sd_N))