knitr::opts_chunk$set(echo = TRUE, message = FALSE, 
                      warning = FALSE)

Проблема динамики численности мелких млекопитающих занимает центральное место в популяционных исследованиях. Однако закономерности и причины изменения обилия насекомоядных изучены в меньшей степени, чем для грызунов. Территория Восточно-Уральского радиоактивного следа (ВУРСа) является природным полигоном с техногенным загрязнением, на которой актуально изучение состояния и путей адаптации популяций животных.

Отловы животных проводились сотрудниками лаборатории эволюционной экологии ИЭРиЖ на импактном (район старой лежневой дороги в головной части следа) и условно контрольном (окр. пос. Метлино, юго-восточный берег оз. Кожакуль) участках по стандартной методике с использованием ловушек-плашек.

Иллюстрации

Район исследования

Объект исследования

Цель:

Описать видовой состав и динамику численности насекомоядных (бурозубок) за 20-летний период наблюдений с 2003 по 2023 гг. на территории Восточно-Уральского радиоактивного следа.

Задачи:

  1. Описать динамику обилия и видовой состав насекомоядных на контрольной и импактной территориях ВУРСа.
  2. Рассмотреть половую структуру сообществ.
  3. Сравнить изменения обилия доминирующих видов Sorex araneus Linnaeus, 1758 и S. caecutiens Laxmann, 1788 на контрольной и импактной территориях.

Таблица с исходными данными

Таблица данных 
library(readxl)
S <- read_excel("Sorex.xlsx")
formattable::formattable(S)
year Control Impact C_lovushka I_lovushka C_male C_female I_male I_female C_S.araneus I_S.araneus C_S.caec. I_S.caec.
2003 32 38 970 1500 11 21 19 19 31 32 0 6
2004 6 1 900 1456 4 2 1 0 6 1 0 0
2005 55 34 1250 1524 18 37 18 16 53 33 2 0
2006 107 144 776 1970 48 59 83 61 91 134 14 5
2007 12 1 1440 2120 9 3 1 0 12 1 0 0
2008 15 14 400 788 4 11 6 8 15 12 0 1
2009 3 24 717 3547 1 2 14 10 3 24 0 0
2010 45 17 665 882 27 18 8 9 27 7 18 7
2011 7 0 500 1300 4 3 0 0 5 0 2 0
2012 0 0 500 490 0 0 0 0 0 0 0 0
2013 2 1 400 1312 1 1 1 0 2 1 0 0
2014 21 3 400 1200 13 8 3 0 15 2 5 1
2015 5 0 400 400 4 1 0 0 4 0 1 0
2016 36 13 300 300 27 9 5 8 23 4 9 9
2017 5 0 298 300 2 3 0 0 5 0 0 0
2018 7 0 400 400 4 3 0 0 4 0 3 0
2019 28 0 300 500 16 12 0 0 21 0 6 0
2020 29 0 488 500 16 13 0 0 12 0 16 0
2021 0 0 400 400 0 0 0 0 0 0 0 0
2022 15 0 400 400 6 9 0 0 7 0 7 0
2023 51 0 400 400 16 35 0 0 20 0 28 0

Анализ данных:

1. Динамика обилия

Грызуны и бурозубки относятся к цикломорфным животным, период колебания их численности составляет обычно около трех лет. За 20 лет наблюдений отмечено 6–7 циклов.

library(readxl)
library(tidyverse)
S <- read_excel("Sorex.xlsx")

S1 <- S %>% 
    select(year:I_lovushka) %>% 
    transmute(year,
              Control = Control/C_lovushka*100, 
              Impact = Impact/I_lovushka*100
    ) %>% 
    pivot_longer(names_to = "zone", values_to = "abu", -year)

    G <- ggplot(S1,aes(x = year, y = abu, color = zone))+
    geom_point() +
    geom_line()+ 
    scale_color_manual(values = c("green", "red")) +
    scale_x_continuous(breaks = 2003:2023) +
    labs(x = NULL, y = "Обилие, экз./100 л.-с",
         title = "Обилие насекомоядных по годам, экз./100 л.-с") + 
    theme_classic() +
    theme(
        legend.position = "bottom",
        axis.text.x = element_text(angle = 45))
plotly::ggplotly(G)

Динамика обилия на контрольном и импактом участках в целом совпадают, фазы высокой и низкой численноси одинаковы.

library(readxl)
library(tidyverse)
S <- read_excel("Sorex.xlsx")
tmp <- S %>% 
    select(year:I_lovushka) %>% 
    transmute(year,
              Control = Control/C_lovushka*100, 
              Impact = Impact/I_lovushka*100
    )
wilcox.test(x = tmp$Control, y = tmp$Impact,
            alternative = "two.sided")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  tmp$Control and tmp$Impact
## W = 348, p-value = 0.001227
## alternative hypothesis: true location shift is not equal to 0

Для подтверждения различий на участках посчитали тест Вилкоксона. Значение теста составляет 348 , а соответствующее значение p равно 0.0012.Поскольку это p-значение меньше 0,05, мы можем отвергнуть нулевую гипотезу, следовательно местообитания статистически значимо различаются между собой.

library(readxl)
library(tidyverse)
S <- read_excel("Sorex.xlsx")
C <- S %>% 
  select(year:I_lovushka) %>% 
  transmute(year,
            Control = Control/C_lovushka*100, 
            Impact = Impact/I_lovushka*100
  ) %>% 
  pivot_longer(names_to = "zone", values_to = "abu", -year)

A1 <- aov(abu~zone, data=C)
summary(A1) 
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## zone         1  111.1  111.06   9.879 0.00315 **
## Residuals   40  449.7   11.24                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(A1)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = abu ~ zone, data = C)
## 
## $zone
##                     diff       lwr       upr     p adj
## Impact-Control -3.252224 -5.343488 -1.160959 0.0031455
residuals(A1) %>% shapiro.test()
## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.85264, p-value = 7.275e-05
M1 <- lm(abu~zone, data=C)
residuals(M1) %>% shapiro.test()
## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.85264, p-value = 7.275e-05
summary(M1)
## 
## Call:
## lm(formula = abu ~ zone, data = C)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.2632 -1.0110 -0.9640  0.9692  9.5255 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.2632     0.7317   5.827 8.28e-07 ***
## zoneImpact   -3.2522     1.0347  -3.143  0.00315 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.353 on 40 degrees of freedom
## Multiple R-squared:  0.1981, Adjusted R-squared:  0.178 
## F-statistic: 9.879 on 1 and 40 DF,  p-value: 0.003145
t <-data.frame(x = M1$residuals) %>% 
   ggplot(aes(x)) + 
  geom_vline(xintercept = 0, 
             color = "blue", linetype = "dashed" ) +
  geom_density() + 
  theme_bw() + 
  labs(x = "Residuals", y = "Density")
plotly::ggplotly(t)

Далее сделали дисперсионный однофакторный анализ ANOVA, который так же подтвердил, значимую разницу двух участков. Линейная модель показыввает нам что R-квадрат равен 0.1981, что указывает на то, что 20% дисперсии можно объяснить местообитанием.

library(readxl)
library(tidyverse)
S <- read_excel("Sorex.xlsx")
C <- S %>% 
  select(year:I_lovushka) %>% 
  transmute(year,
            Control = Control/C_lovushka*100, 
            Impact = Impact/I_lovushka*100
  ) %>% 
  pivot_longer(names_to = "zone", values_to = "abu", -year)
r <- cbind(C,fitted = M1$fitted.values) %>%
   ggplot()+
  geom_abline(aes(intercept = M1$coefficients[[1]], slope = M1$coefficients[[2]]))+
  geom_segment(mapping = aes(x = zone, xend = zone, y = abu, yend =fitted, ), color = "blue")+            
  geom_point(mapping = aes(x = zone, y = abu), shape = 21, fill = "white")+
  labs(x = "Местообитание", y = "Обилие", 
       title = "Зависимость обилия от местообитания")+
  theme_bw()
plotly::ggplotly(r)

Распределение остатков ненормальное согласно тесту Шапиро-Уилка.

2. Половая структура сообществ

knitr::opts_chunk$set(echo = FALSE, message = FALSE, 
                      warning = FALSE)
S2 <- S %>% 
  select(-Control,- Impact, -C_lovushka, -I_lovushka,
         -C_S.araneus, -I_S.araneus, -C_S.caec. ,-I_S.caec.
  ) %>% 
transmute(year,
          male=(S$C_male/(S$C_male+S$C_female))*100 , 
          female=(S$C_female/(S$C_male+S$C_female))*100
) %>% 
  pivot_longer(names_to = "Sex", values_to = "abu", -year) 
F <- ggplot(S2,aes(x = year, y = abu, fill = Sex))+
  geom_col()+
  scale_x_continuous(breaks = 2003:2023) +
  labs(x = NULL, y = "Соотношение полов, %",
       title = "Половая структура насекомоядных на контроле, %") + 
  theme_classic() +
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(angle = 45))
  plotly::ggplotly(F)

Обобщенный график.

O <- S %>%
  select(-C_lovushka, -I_lovushka, -C_S.araneus, -I_S.araneus, -C_S.caec., -I_S.caec.)%>%
  mutate(C_male=(S$C_male/(S$C_male+S$C_female))*100 , 
         C_female=(S$C_female/(S$C_male+S$C_female))*100,
         I_male=(S$I_male/(S$I_male+S$I_female))*100 ,
         I_female=(S$I_female/(S$I_male+S$I_female))*100
  ) %>%
  pivot_longer(
    cols = c(C_male, C_female, I_male, I_female),
    names_to = 'Sex',
    values_to = 'abu'
  )
  E <- ggplot(O,aes(x = year, y = abu, fill = Sex)) +
  geom_col(position = "stack") +
  scale_x_discrete(breaks = 2003:2023) +
  labs(x = NULL, y = "Соотношение полов, %",
       title = "Половая структура насекомоядных на контроле и импакте, %") + 
  theme_classic() +
  theme(
    legend.position = "bottom")
  plotly::ggplotly(E)

Соотношение 1 : 1 почти не характерно для отловленных особей.

W1 <-sum(S$C_male)
W2 <-sum(S$C_female)
Y=W1+W2
R1 =W1/Y*100
R2=W2/Y*100
W3 <-sum(S$I_male)
W4 <-sum(S$I_female)
Y1=W3+W4
R3 =W3/Y1*100
R4=W4/Y1*100
u<- c(R1, R2, R3, R4)
print(u)
## [1] 48.02495 51.97505 54.82759 45.17241

По общим подсчетам заметна такая тенденция: на контрольном участке доля самок немного выше (52 %), чем самцов (48 %) а на импакте наоборот доля самцов выше (55 %), чем самок (45%).

3. Динамика обилия доминирующих видов Sorex araneus и S. caecutiens

Видовое разнообразие

Видовой состав насекомоядных на исследуемой территории включает 6 видов рода Sorex и 1 вид рода Neomys. К видам-доминантам относятся обыкновенная бурозубка Sorex araneus и средняя бурозубка S. caecutiens. Также отмечены малая бурозубка S. minutus, тундряная бурозубка S. tundrensis, равнозубая бурозубка S. isodon, крупнозубая бурозубка S. daphaenodon и водяная кутора Neomys fodiens. Видовой состав контрольной территории представлен 6 видами, а импактой территории 4 видами.

Таблица видов 
library(readxl)
T <- read_excel("table.xlsx")
formattable::formattable(T)
Контроль …2 …3 …4 …5 …6 …7 …8 …9 …10 …11 …12 …13 …14 …15 …16 …17 …18 …19 …20 …21 …22
Таксон 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023
Sorex araneus 31 6 53 91 12 15 3 27 5 0 2 15 4 23 5 4 21 12 0 7 20
Sorex caecutiens 0 0 2 14 0 0 0 18 2 0 0 5 1 9 0 3 6 16 0 7 28
Sorex isodon 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Sorex tundrensis 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
Sorex minutus 0 0 0 2 0 0 0 0 0 0 0 1 0 3 0 0 0 0 0 1 3
Sorex daphaenodon 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
Neomys fodiens 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
Итого 32 6 55 108 12 15 3 45 7 0 2 21 5 36 5 7 28 29 0 15 51
Импакт NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Таксон 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 2022
Sorex araneus 32 1 33 134 1 12 24 7 0 0 1 2 0 4 0 0 0 0 0 0 0
Sorex caecutiens 6 0 0 5 0 1 0 7 0 0 0 1 0 9 0 0 0 0 0 0 0
Sorex isodon 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Sorex tundrensis 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Sorex minutus 0 0 0 5 0 1 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0
Sorex daphaenodon 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Neomys fodiens 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Итого 38 1 34 144 1 14 24 17 0 0 1 3 0 13 0 0 0 0 0 0 0 

H <- S %>% 
  transmute(year,
            S.ar_C=(S$C_S.araneus/S$C_lovushka)*100 , 
            S.ar_I=(S$I_S.araneus/S$I_lovushka)*100,
            S.caec_C=(S$C_S.caec./S$C_lovushka)*100 , 
            S.caec_I=(S$I_S.caec./S$I_lovushka)*100
   )%>% 
  pivot_longer(names_to = "species", values_to = "abu", -year)
  l <- ggplot(H,aes(x = year, y = abu, color = species))+
  geom_point() +
  geom_line()+ 
  scale_color_manual(values = c("green", "red", "blue","orange")) +
  scale_x_continuous(breaks = 2003:2023) +
  labs(x = NULL, y = "Обилие, экз./100 л.-с",
       title = "Динамика обилия доминирующих видов Sorex araneus и S. caecutiens, экз./100 л.-с") + 
  theme_classic() +
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(angle = 45))
 plotly::ggplotly(l)

Фазы обилия животных двух видов бурозубок в целом довольно синхронны на участках, что может косвенно свидетельствовать об отсутствии конкуренции между ними, в том числе пищевой, так как обыкновенная бурозубка значительно крупнее средней, следовательно, способна охотиться на более крупных насекомых, и пищевые ресурсы этих видов не пресекаются. Более низкую встречаемость средней бурозубки в отловах можно объяснить опять же ее меленькими размерами, и тем, что она реже попадается в ловушки-давилки. Однако, нельзя отрицать и в целом более низкую численность средней бурозубки по сравнению с обыкновенной.

J <- S %>% 
  transmute(year,
            S.ar_C=(S$C_S.araneus/S$C_lovushka)*100 , 
            S.ar_I=(S$I_S.araneus/S$I_lovushka)*100,
            S.caec_C=(S$C_S.caec./S$C_lovushka)*100 , 
            S.caec_I=(S$I_S.caec./S$I_lovushka)*100
  )  
cor(J[, c('S.ar_C', 'S.ar_I', 'S.caec_C', 'S.caec_I')])
##             S.ar_C     S.ar_I   S.caec_C  S.caec_I
## S.ar_C   1.0000000 0.74529900 0.48606174 0.4427044
## S.ar_I   0.7452990 1.00000000 0.01233548 0.1883452
## S.caec_C 0.4860617 0.01233548 1.00000000 0.2726956
## S.caec_I 0.4427044 0.18834518 0.27269555 1.0000000

Рассчитали коэфициент Пирсона, по результатам можно сделать вывод, что динамики обилия Sorex araneus и S. caecutiens на контроле имеют среднюю связь (r=0.49), а на импакте связь отсутствуют (r=0.18). Это дает основание полагать, что между этими видами нет явной конкуренции.

tmp1 <- S %>% 
  transmute(year,
            Aran=(S$C_S.araneus/S$C_lovushka)*100 , 
            Caec=(S$C_S.caec./S$C_lovushka)*100) 
wilcox.test(x = tmp1$Aran, y = tmp1$Caec,
            alternative = "two.sided")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  tmp1$Aran and tmp1$Caec
## W = 324, p-value = 0.008917
## alternative hypothesis: true location shift is not equal to 0
P <- S %>% 
  transmute(year,
            Aran=(S$C_S.araneus/S$C_lovushka)*100 , 
            Caec=(S$C_S.caec./S$C_lovushka)*100) %>%
  pivot_longer(names_to = "species", values_to = "abu", -year)
A2 <- aov(abu~species, data=P)
summary(A2) 
##             Df Sum Sq Mean Sq F value Pr(>F)  
## species      1  33.21   33.21   5.561 0.0233 *
## Residuals   40 238.87    5.97                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(A2)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = abu ~ species, data = P)
## 
## $species
##                diff       lwr        upr     p adj
## Caec-Aran -1.778343 -3.302527 -0.2541581 0.0233459
residuals(A2) %>% shapiro.test()
## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.84608, p-value = 4.999e-05
M2 <- lm(abu~species, data=P)
residuals(M2) %>% shapiro.test()
## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.84608, p-value = 4.999e-05
summary(M2)
## 
## Call:
## lm(formula = abu ~ species, data = P)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9378 -1.1808 -0.9545  0.8334  8.7890 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.9378     0.5333   5.509 2.31e-06 ***
## speciesCaec  -1.7783     0.7541  -2.358   0.0233 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.444 on 40 degrees of freedom
## Multiple R-squared:  0.122,  Adjusted R-squared:  0.1001 
## F-statistic: 5.561 on 1 and 40 DF,  p-value: 0.02335
n <- data.frame(x = M2$residuals) %>% 
  ggplot(aes(x)) + 
  geom_vline(xintercept = 0, 
             color = "blue", linetype = "dashed" ) +
  geom_density() + 
  theme_bw() + 
  labs(x = "Residuals", y = "Density")
 plotly::ggplotly(n)
P <- S %>% 
  transmute(year,
            Aran=(S$C_S.araneus/S$C_lovushka)*100 , 
            Caec=(S$C_S.caec./S$C_lovushka)*100) %>%
  pivot_longer(names_to = "species", values_to = "abu", -year)
f <- cbind(P,fitted = M2$fitted.values) %>%
  ggplot()+
  geom_abline(aes(intercept = M2$coefficients[[1]], slope = M2$coefficients[[2]]))+
  geom_segment(mapping = aes(x = species, xend =species , y = abu, yend =fitted, ), color = "blue")+            
  geom_point(mapping = aes(x = species, y = abu), shape = 21, fill = "white")+
  labs(x = "Вид", y = "Обилие", 
       title = "Зависимость обилия от вида")+
  theme_bw()
 plotly::ggplotly(f)

Проверили с помощью критерия Вилкоксона достоверно ли различаются виды на контроле и импакте. На контроле виды достоверно различаются, далее сделали дисперсионный анализ ANOVA, проверили на нормальность остатки, линейная модель показыввает нам что R-квадрат равен 0.122, что указывает на то, что только 12% дисперсии можно объяснить видовой принадлеженостью. Однако распределение остатков ненормальное согласно тесту Шапиро-Уилка.

На импактоном участке, значимых различий между видами не выявлено, согласно критерию Вилкоксона.

tmp2 <- S %>% 
  transmute(year,
            Aran=(S$I_S.araneus/S$I_lovushka)*100 , 
            Caec=(S$I_S.caec./S$I_lovushka)*100) 

wilcox.test(x = tmp2$Aran, y = tmp2$Caec,
            alternative = "two.sided")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  tmp2$Aran and tmp2$Caec
## W = 276.5, p-value = 0.1161
## alternative hypothesis: true location shift is not equal to 0

Выводы

  1. Динамика обилия насекомоядных сходна на контрольной и импактной территориях ВУРСа за изученный 20-летний период , фазы высокой и низкой численности совпадают.
  2. Половая структура неоднородна от года к году, но было выявлено, что на контрольном участке преобладают в основном самки, а на импакте – самцы.
  3. У доминирующих видов обыкновенной и средней бурозубок - динамика обилия имеет синхронный характер, что косвенно может свидетельствовать об отсутствии конкуренции между ними.