Fish growth and survival

Fish Growth:

Average fish growth (mass, g) over time by bird treatment and study year:

Average specific growth rate of PIT tagged fish by treatment and study year, calculated based on Crane et al 2019, Reviews in Aquaculture:

Standard mass and length specific growth rates
BirdTreatment StudyYear AvgMassSpecificGrowthRate AvgLengthSpecificGrowthRate n
Birds Year1 0.6126579 0.1518465 809
Birds Year2 0.2683540 -0.0160723 760
NoBirds Year1 0.5506080 0.0825875 29
NoBirds Year2 0.2684678 0.0028404 22

Fish Survival:

Average survival by treatment and study year:

Study Year Bird Treatment Mean survival (recaptures)
1 Birds 24%
1 No Birds 69%
2 Birds 33%
2 No Birds 69%

Next steps:

Mark-recapture models to better estimate survival, and use these estimates as predictors in daphnia & GHG models.

Progress to-date:

  • Created full datasets for year 1 and year 2 in capture-history format

  • Read up on mark-recapture models

  • Familiarized/practiced with running CJS models in “marked” R Package

Challenges:

  • Low detection may cause biased (low) survival estimates, inflate uncertainty and create identifiability issues

  • Ways to statsitically work through this: let detection probability vary by sampling occasion in the model, include covariates in the model that may affect detection (effort, temp, depth, etc.), use Bayesian or robust models

Daphnia spp. abundance

Visualizations, means & year comparisons

Daphnia

Grey dashed lines represent dates of fish introduction in year 1 and year 2 of study.

Year 1: average total daphnia abundances

Mean estimated daphnia abundance (ind/m^3) in Year 1
treatment Pre_Post mean_abundance sd_abundance
Birds Fish Post 786.4473 1045.2053
Birds Fish Pre 1095.5727 713.8342
Birds NoFish Post 1257.6464 1923.3930
Birds NoFish Pre 464.1978 376.9683
NoBirds Fish Post 383.3047 421.6185
NoBirds Fish Pre 669.0987 687.4315
NoBirds NoFish Post 825.3686 953.7184
NoBirds NoFish Pre 760.3012 488.8243

Year 2: average total daphnia abundances

Mean estimated daphnia abundance (ind/m^3) in Year 2
treatment Pre_Post mean_abundance sd_abundance
Birds Fish Post 1821.623 1882.926
Birds Fish Pre 3017.717 2509.496
Birds NoFish Post 2257.243 2324.635
Birds NoFish Pre 4088.981 5140.673
NoBirds Fish Post 1805.652 2028.517
NoBirds Fish Pre 5919.426 5709.769
NoBirds NoFish Post 2263.500 2187.454
NoBirds NoFish Pre 2821.677 3182.385

Year Comparison: total daphnia abundances

Based on looking at the raw data and averages, we hypothesize that there were significantly higher daphnia abundances in year 2 compared to year 1. Below is a generalized linear mixed effects model using a negative binomial distribution investigating daphnia abundance by study year, with plot ID as a random effect to account for repeated sampling. Here, we find evidence for our hypothesis, as daphnia abundance in year 2 is significantly higher than in year 1 (p >> 0.001). Model diagnostics and residuals were checked using the DHARMa package.

##  Family: nbinom2  ( log )
## Formula:          total ~ StudyYear + (1 | PlotID)
## Data: daphnia
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    3644.7    3658.2   -1818.4    3636.7       213 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  PlotID (Intercept) 0.1182   0.3439  
## Number of obs: 217, groups:  PlotID, 16
## 
## Dispersion parameter for nbinom2 family (): 1.02 
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      6.6188     0.1369   48.36   <2e-16 ***
## StudyYearYear2   1.2472     0.1395    8.94   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Zooplankton

Grey dashed lines represent dates of fish introduction in year 1 and year 2 of study.

Year 1: average zooplankton abundance

Mean estimated zooplankton abundance (ind/m^3) in Year 1
treatment Pre_Post mean_abundance sd_abundance
Birds Fish Post 5263.1297 4294.01404
Birds Fish Pre 1071.4005 420.16931
Birds NoFish Post 4166.6716 4686.61146
Birds NoFish Pre 529.1473 286.46029
NoBirds Fish Post 2331.9024 1565.22807
NoBirds Fish Pre 724.0964 113.92453
NoBirds NoFish Post 6154.0283 14968.74086
NoBirds NoFish Pre 899.3866 97.30572

Year 2: average zooplankton abundance

Mean estimated zooplankton abundance (ind/m^3) in Year 2
treatment Pre_Post mean_abundance sd_abundance
Birds Fish Post 12725.954 8947.446
Birds Fish Pre 5805.877 4135.730
Birds NoFish Post 12381.130 7807.918
Birds NoFish Pre 6397.340 5833.891
NoBirds Fish Post 14516.809 10333.944
NoBirds Fish Pre 10384.926 8830.063
NoBirds NoFish Post 17499.054 11396.432
NoBirds NoFish Pre 4853.957 3573.115

Year comparison: zooplankton abundance

As with daphnia, we hypothesized total zooplankton abundances would be significantly higher in year 2 compared to year 1. Below is a generalized linear mixed effects model using a negative binomial distribution investigating zooplankton abundance by study year, with plot ID as a random effect to account for repeated sampling. Here, we similarly find evidence for our hypothesis, as zooplankton abundance in year 2 is significantly higher than in year 1 (p >> 0.001). Model diagnostics and residuals were checked using the DHARMa package.

##  Family: nbinom2  ( log )
## Formula:          DensityRounded ~ StudyYear + (1 | PlotID)
## Data: zoop
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    4296.6    4310.1   -2144.3    4288.6       213 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  PlotID (Intercept) 0.05713  0.239   
## Number of obs: 217, groups:  PlotID, 16
## 
## Dispersion parameter for nbinom2 family (): 1.33 
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      8.2125     0.1132   72.52   <2e-16 ***
## StudyYearYear2   1.1430     0.1263    9.05   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Year 1: Models

Model comparison: daphnia ~ pre/post fish introductions:

Below is a table of model comparisons for daphnia abundance models. We investigated total daphnia abundance by fish treatment, bird treatment and pre/post fish introductions. We also included a post-hoc analysis investigating daphnia abundance after removing the final sampling date, due to patterns we saw in the raw data (see raw data/visualization tab) where daphnia abundance increased drastically in year 1. We hypothesize this may be due to fish moralities towards the end of the study, due to either Saprolegnia spp. infections or bird predation. All models are GLMMs with negative binomial distributions and log link functions, and include water temperature, dissolved oxygen and depth as covariates. Plot ID and sampling occassion are included as random effects for all models.

In the model including all sampling dates, we see a marginal effect of fish x pre/post, with a decrease in total daphnia abundance in fish plots after fish were introduced. When we exclude the final sampling date, this effect is significant (p < 0.05).

Daphnia models, year 1 (full flooding period and post-hoc analysis removing final sampling day)
Response Model Log Likelihood Comparison X2 P-value
Daphnia abundance: year 1, including final sampling date M1: fish x birds x pre/post -674.60 M1 vs M2 1.96 0.162

M2: (fish x pre/post) +

(birds x pre/post) +

(fish x birds)

 -675.68 M2 vs M3 1.17 0.279

M3: (fish x pre/post) +

(birds x pre/post)

 -676.26 M3 vs M4 1.53 0.216

M4: (fish x pre/post) +

birds

 -677.03 M4 vs M5 1.96 0.161
M5: fish x pre/post -678.01 M5 vs M6 3.48 0.062 ·
M6: fish + pre/post -679.75 M6 vs M7 0.07 0.797
M7: fish only -679.75 M7 vs M8 0.32 0.569
M8: no fish or pre/post -679.95 NA NA NA
Daphnia abundance: year 1, omitting final sampling date M1: fish x birds x pre/post -536.28 M1 vs M2 2.89 0.089 ·
M2: (fish x pre/post) + (birds x pre/post) + (fish x birds) -537.72 M2 vs M3 0.37 0.541
M3: (fish x pre/post) + (birds x pre/post) -537.91 M3 vs M4 1.48 0.223
M4: (fish x pre/post) + birds -538.65 M4 vs M5 1.16 0.282
M5: fish * pre/post -539.23 M5 vs M6 5.81 0.016 *
M6: fish + pre/post -542.12 M6 vs M7 0.44 0.507
M7: fish only -542.35 M7 vs M8 0.13 0.723
M8: no fish or pre/post -542.35 NA NA NA

Model outputs: daphnia ~ pre/post fish introductions:

Fish*pre/post, including final sampling date:
##  Family: nbinom2  ( log )
## Formula:          DensityRounded ~ FishTreatment * Pre_Post + MeanTempScaled +  
##     MeanDOScaled + MeanDepthScaled + (1 | PlotID) + (1 | SamplingOccasion)
## Data: totaldaphnia1
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##      1376      1401      -678      1356        80 
## 
## Random effects:
## 
## Conditional model:
##  Groups           Name        Variance Std.Dev.
##  PlotID           (Intercept) 0.1874   0.4329  
##  SamplingOccasion (Intercept) 0.2624   0.5123  
## Number of obs: 90, groups:  PlotID, 16; SamplingOccasion, 6
## 
## Dispersion parameter for nbinom2 family ():  1.4 
## 
## Conditional model:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                      6.27119    0.35492  17.670   <2e-16 ***
## FishTreatmentNoFish              0.41303    0.33000   1.252   0.2107    
## Pre_PostPre                      0.24510    0.58153   0.421   0.6734    
## MeanTempScaled                  -0.02537    0.24787  -0.102   0.9185    
## MeanDOScaled                    -0.43242    0.18054  -2.395   0.0166 *  
## MeanDepthScaled                  0.07517    0.14228   0.528   0.5972    
## FishTreatmentNoFish:Pre_PostPre -0.81462    0.42965  -1.896   0.0580 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## FishTreatment = Fish:
##  contrast   ratio    SE  df null z.ratio p.value
##  Post / Pre 0.783 0.455 Inf    1  -0.421  0.6734
## 
## FishTreatment = NoFish:
##  contrast   ratio    SE  df null z.ratio p.value
##  Post / Pre 1.767 1.020 Inf    1   0.987  0.3236
## 
## Tests are performed on the log scale

Fish*pre/post, omitting final sampling date:
##  Family: nbinom2  ( log )
## Formula:          DensityRounded ~ FishTreatment * Pre_Post + MeanTempScaled +  
##     MeanDOScaled + MeanDepthScaled + (1 | PlotID) + (1 | SamplingOccasion)
## Data: totaldaphnia1_filtered
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1098.5    1121.5    -539.2    1078.5        64 
## 
## Random effects:
## 
## Conditional model:
##  Groups           Name        Variance  Std.Dev. 
##  PlotID           (Intercept) 2.214e-01 4.706e-01
##  SamplingOccasion (Intercept) 1.827e-09 4.274e-05
## Number of obs: 74, groups:  PlotID, 16; SamplingOccasion, 5
## 
## Dispersion parameter for nbinom2 family (): 1.47 
## 
## Conditional model:
##                                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                      5.792047   0.261614  22.140  < 2e-16 ***
## FishTreatmentNoFish              0.590500   0.374425   1.577  0.11478    
## Pre_PostPre                      0.666370   0.336810   1.978  0.04788 *  
## MeanTempScaled                  -0.193672   0.113435  -1.707  0.08776 .  
## MeanDOScaled                    -0.388451   0.134183  -2.895  0.00379 ** 
## MeanDepthScaled                 -0.005061   0.139696  -0.036  0.97110    
## FishTreatmentNoFish:Pre_PostPre -1.090539   0.440112  -2.478  0.01322 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## FishTreatment = Fish:
##  contrast   ratio    SE  df null z.ratio p.value
##  Post / Pre 0.514 0.173 Inf    1  -1.978  0.0479
## 
## FishTreatment = NoFish:
##  contrast   ratio    SE  df null z.ratio p.value
##  Post / Pre 1.528 0.532 Inf    1   1.218  0.2234
## 
## Tests are performed on the log scale

Model comparison: zooplankton ~ pre/post fish introductions:

Below is a table of model comparisons for total zooplankton abundance models. We followed the same model selection structure as in our daphnia models. We hypothesized that we would see a stronger overall effect when we subset the data by family = daphniidae, as we hypothesize fish will preferentially consume these large bodied cladocerans, and because daphnia are known to exert strong top-down control of methane oxidizing bacteria.

Here, when we include the final sampling date, we see a marginal effect of pre/post fish introduction. When we omit the final sampling date, we see a marginal effect of fish x birds interaction, birds x pre/post, and a significant effect of pre/post fish introductions. However, we don’t see clear effects of fish suppressing overall zooplankton abundance.

Full zooplankton models, year 1 (full flooding period and post-hoc analysis removing final sampling day). There is no evidence in either model set that treatment affects total zooplankton abundance. There is evidence, however, that total zooplankton abundance grew post fish introduction - especially when the final sampling date was removed.
Response Model Log Likelihood Comparison X2 P-value
Zooplankton abundance: year 1, including final sampling date M1: fish x birds x pre/post -808.93 M1 vs M2 0.26 0.607

M2: (fish x pre/post) +

(birds x pre/post) +

(fish x birds)

 -808.83 M2 vs M3 2.10 0.147

M3: (fish x pre/post) +

(birds x pre/post)

 -810.02 M3 vs M4 0.65 0.421

M4: (fish x pre/post) +

birds

 -810.34 M4 vs M5 1.75 0.186
M5: fish x pre/post -811.21 M5 vs M6 0.55 0.459
M6: fish + pre/post -811.49 M6 vs M7 2.97 0.084 ·
M7: fish only -812.97 M7 vs M8 0.01 0.936
M8: no fish or pre/post -812.98 NA NA NA
Zooplankton abundance: year 1, omitting final sampling date M1: fish x birds x pre/post -639.53 M1 vs M2 0.49 0.482
M2: (fish x pre/post) + (birds x pre/post) + (fish x birds) -639.78 M2 vs M3 2.85 0.091 ·
M3: (fish x pre/post) + (birds x pre/post) -640.07 M3 vs M4 0.36 0.551
M4: (fish x pre/post) + birds -641.39 M4 vs M5 3.72 0.054 ·
M5: fish x pre/post -643.24 M5 vs M6 0.02 0.880
M6: fish + pre/post -643.26 M6 vs M7 12.71 <0.001***
M7: fish only -649.61 M7 vs M8 0.031 0.861
M8: no fish or pre/post -649.63 NA NA NA

Year 2: Models

Model comparison: daphnia ~ pre/post fish introductions:

Below is a model comparison for year 2 daphnia models. Unlike in year 1, we are not including a post-hoc analysis removing the final sampling date, because the raw data from year 2 doesn’t show a pattern of large daphnia increases at the end of the field season, as we saw in year 1. All models are GLMMs with negative binomial distributions and log link functions, and include water temperature, dissolved oxygen and depth as covariates. Plot ID and sampling occassion are included as random effects for all models.

Here, we see some evidence of an effect of the three-way interaction between fish, birds, and pre/post fish introduction. We also see a marginal effect of the fish & pre/post interaction. When investigating post-hoc comparisons of estimated marginal means, we see evidence of fish reducing daphnia after introduction, but only in the absence of birds.

Daphnia model selection, year 2
Response Model Log Likelihood Comparison X2 P-value
Daphnia abundance: year 2 M1: fish x birds x pre/post -1114.1 M1 vs M2 4.39 0.036 *
M2: (fish x pre/post) + (birds x pre/post) + (fish x birds) -1116.3 Ms vs M3 0.32 0.573
M3: (fish * pre/post) + (birds * pre/post) -1116.4 M3 vs M4 0.33 0.566
M4: (fish x pre/post) + birds -1116.6 M4 vs M5 0.008 0.927
M5: fish x pre/post -1116.6 M5 vs M6 3.56 0.059 ·
M6: fish + pre/post -1118.4 M6 vs M7 0.29 0.590
M7: fish only -1118.5 M7 vs M8 0.01 0.926
M8: no fish or pre/post -1118.5 NA NA NA 
##  Family: nbinom2  ( log )
## Formula:          
## DensityRounded ~ FishTreatment * BirdTreatment * Pre_Post + MeanTempScaled +  
##     MeanDOScaled + MeanDepthScaled + (1 | PlotID) + (1 | SamplingOccasion)
## Data: totaldaphnia2
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    2247.6    2287.4   -1109.8    2219.6       113 
## 
## Random effects:
## 
## Conditional model:
##  Groups           Name        Variance Std.Dev.
##  PlotID           (Intercept) 0.1022   0.3197  
##  SamplingOccasion (Intercept) 0.2575   0.5074  
## Number of obs: 127, groups:  PlotID, 16; SamplingOccasion, 8
## 
## Dispersion parameter for nbinom2 family (): 1.52 
## 
## Conditional model:
##                                                      Estimate Std. Error
## (Intercept)                                           7.48124    0.35918
## FishTreatmentNoFish                                   0.19126    0.35244
## BirdTreatmentNoBirds                                 -0.20758    0.36538
## Pre_PostPre                                           0.42916    0.52548
## MeanTempScaled                                        0.06541    0.13963
## MeanDOScaled                                         -0.05650    0.10124
## MeanDepthScaled                                      -0.11337    0.12062
## FishTreatmentNoFish:BirdTreatmentNoBirds              0.13451    0.50813
## FishTreatmentNoFish:Pre_PostPre                      -0.09691    0.44004
## BirdTreatmentNoBirds:Pre_PostPre                      0.62363    0.43799
## FishTreatmentNoFish:BirdTreatmentNoBirds:Pre_PostPre -0.94658    0.61677
##                                                      z value Pr(>|z|)    
## (Intercept)                                           20.829   <2e-16 ***
## FishTreatmentNoFish                                    0.543    0.587    
## BirdTreatmentNoBirds                                  -0.568    0.570    
## Pre_PostPre                                            0.817    0.414    
## MeanTempScaled                                         0.468    0.639    
## MeanDOScaled                                          -0.558    0.577    
## MeanDepthScaled                                       -0.940    0.347    
## FishTreatmentNoFish:BirdTreatmentNoBirds               0.265    0.791    
## FishTreatmentNoFish:Pre_PostPre                       -0.220    0.826    
## BirdTreatmentNoBirds:Pre_PostPre                       1.424    0.154    
## FishTreatmentNoFish:BirdTreatmentNoBirds:Pre_PostPre  -1.535    0.125    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## FishTreatment = Fish, BirdTreatment = Birds:
##  contrast   ratio    SE  df null z.ratio p.value
##  Post / Pre 0.651 0.342 Inf    1  -0.817  0.4141
## 
## FishTreatment = NoFish, BirdTreatment = Birds:
##  contrast   ratio    SE  df null z.ratio p.value
##  Post / Pre 0.717 0.362 Inf    1  -0.658  0.5107
## 
## FishTreatment = Fish, BirdTreatment = NoBirds:
##  contrast   ratio    SE  df null z.ratio p.value
##  Post / Pre 0.349 0.178 Inf    1  -2.060  0.0394
## 
## FishTreatment = NoFish, BirdTreatment = NoBirds:
##  contrast   ratio    SE  df null z.ratio p.value
##  Post / Pre 0.991 0.502 Inf    1  -0.018  0.9854
## 
## Tests are performed on the log scale

Model comparison: zooplankton ~ pre/post fish introductions:

Below is a model comparison for year 2 models using our full zooplankton dataset as the response variable. Again, we see a significant effect of the three-way fish x bird x pre/post interaction, and we see significant effects of the fish x pre/post interaction. Model structure was the same as models above.

When running post-hoc comparisons, we see evidence that zooplankton abundance increased significantly post fish introduction in all treatments but one - fish/no birds.

Full zooplankton model selection, year 2
Response Model Log Likelihood Comparison X2 P-value
Zooplankton abundance: year 2 M1: fish x birds x pre/post -1248.9 M1 vs M2 4.75 0.029 *
M2: (fish x pre/post) + (birds x pre/post) + (fish x birds) -1251.3 Ms vs M3 0.70 0.403
M3: (fish * pre/post) + (birds * pre/post) -1251.6 M3 vs M4 0.03 0.860
M4: (fish x pre/post) + birds -1251.7 M4 vs M5 0.41 0.519
M5: fish x pre/post -1251.9 M5 vs M6 4.23 0.039 *
M6: fish + pre/post -1254.0 M6 vs M7 4.53 0.033 *
M7: fish only -1256.2 M7 vs M8 0.19 0.661
M8: no fish or pre/post -1256.3 NA NA NA 
##  Family: nbinom2  ( log )
## Formula:          
## DensityRounded ~ FishTreatment * BirdTreatment * Pre_Post + scale(MeanTemp) +  
##     scale(MeanDO) + scale(MeanDepth) + (1 | PlotID) + (1 | SamplingOccasion)
## Data: totalzoop2
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    2525.8    2565.4   -1248.9    2497.8       111 
## 
## Random effects:
## 
## Conditional model:
##  Groups           Name        Variance Std.Dev.
##  PlotID           (Intercept) 0.01754  0.1325  
##  SamplingOccasion (Intercept) 0.15588  0.3948  
## Number of obs: 125, groups:  PlotID, 16; SamplingOccasion, 8
## 
## Dispersion parameter for nbinom2 family ():  3.1 
## 
## Conditional model:
##                                                      Estimate Std. Error
## (Intercept)                                           9.47386    0.25878
## FishTreatmentNoFish                                   0.05083    0.20747
## BirdTreatmentNoBirds                                  0.02129    0.22365
## Pre_PostPre                                          -0.97545    0.46196
## scale(MeanTemp)                                       0.17903    0.20966
## scale(MeanDO)                                        -0.07859    0.06378
## scale(MeanDepth)                                     -0.05060    0.07768
## FishTreatmentNoFish:BirdTreatmentNoBirds              0.15248    0.30646
## FishTreatmentNoFish:Pre_PostPre                       0.01778    0.29763
## BirdTreatmentNoBirds:Pre_PostPre                      0.43206    0.30873
## FishTreatmentNoFish:BirdTreatmentNoBirds:Pre_PostPre -0.93663    0.42558
##                                                      z value Pr(>|z|)    
## (Intercept)                                            36.61   <2e-16 ***
## FishTreatmentNoFish                                     0.24   0.8065    
## BirdTreatmentNoBirds                                    0.10   0.9242    
## Pre_PostPre                                            -2.11   0.0347 *  
## scale(MeanTemp)                                         0.85   0.3931    
## scale(MeanDO)                                          -1.23   0.2179    
## scale(MeanDepth)                                       -0.65   0.5148    
## FishTreatmentNoFish:BirdTreatmentNoBirds                0.50   0.6188    
## FishTreatmentNoFish:Pre_PostPre                         0.06   0.9524    
## BirdTreatmentNoBirds:Pre_PostPre                        1.40   0.1617    
## FishTreatmentNoFish:BirdTreatmentNoBirds:Pre_PostPre   -2.20   0.0278 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## FishTreatment = Fish, BirdTreatment = Birds:
##  contrast   ratio    SE  df null z.ratio p.value
##  Post / Pre  2.65 1.230 Inf    1   2.112  0.0347
## 
## FishTreatment = NoFish, BirdTreatment = Birds:
##  contrast   ratio    SE  df null z.ratio p.value
##  Post / Pre  2.61 1.190 Inf    1   2.103  0.0355
## 
## FishTreatment = Fish, BirdTreatment = NoBirds:
##  contrast   ratio    SE  df null z.ratio p.value
##  Post / Pre  1.72 0.822 Inf    1   1.138  0.2553
## 
## FishTreatment = NoFish, BirdTreatment = NoBirds:
##  contrast   ratio    SE  df null z.ratio p.value
##  Post / Pre  4.32 2.030 Inf    1   3.101  0.0019
## 
## Tests are performed on the log scale

Estimated MOB activity

Justification/hypotheses:

Here, we are analyzing three stable isotope metrics: 1) stable isotope ratios in porewater, 2) stable isotope ratios in headspace, and 3) the difference between headspace and porewater ratios. All stable isotope analyses are investigating \(\delta\) 13C ( \(\delta\) 13CVPDB), the ratio of carbon-13 to carbon-12 isotopes measured relative to the Vienna Peedee belemnite (VPDB) standard, which is defined as a ratio of exactly 0‰. All porewater and headspace samples were taken in the field and processed at the UC Davis stable isotope facility.

Porewater and headspace isotopic ratios will give us inference into methanogenesis and methane oxidation in different areas of our system. In short, methanogens strongly prefer 12C and thus newly produced methane is depleted and should exhibit a very negative \(\delta\) 13C value (typical rice field methane values are reported between -55‰ and -70‰). Methane oxidizing bacteria, however, preferentially consume 12CH4, removing the lighter isotopes and leaving methane enriched in 13C. Here, \(\delta\) 13C value should be higher (i.e., less negative) if there are more methane oxidizing bacteria in the system. Our hypotheses for our three stable isotope metrics are as follows:

  1. Porewater: samples are taken within the soil (i.e., the anaerobic production zone). We predict that fish primarily increase methane oxidation in oxic zones - because fish should release methane oxidizing bacteria from zooplankton predation pressure in the water column. This means we should see a change in methane oxidation, not production. Thus, we should predict to see little or no differences in porewater \(\delta\) 13C by treatment.

  2. Headspace: Samples are taken as gas is leaving the system during upward transport into the atmosphere. As MOBs oxidize methane, the remaining methane will diffuse into the headspace. Here, the remaining methane should be isotopically heavier because MOBs preferentially oxidize 12CH4 - leaving behind heavier 13CH4. This is evidence that oxidation is occurring between the soil and the atmosphere, and so we hypothesize there will be heavier isotopic signatures in our fish plots when compared to our plots without fish.

  3. Difference between headspace and porewater: Here, because we expect that headspace \(\delta\) 13C should be heavier (less negative) than porewater \(\delta\) 13C, we hypothesize that the difference between the two values (𝛥 \(\delta\) 13C = headspace \(\delta\) 13C - porewater \(\delta\) 13C) should be larger if fish are initiating a trophic cascade that results in shifting oxidation dynamics as methane diffusion is occuring in the oxic zone of the water column. In sum, we would expect to see larger, positive 𝛥 \(\delta\) 13C values in fish plots when compared to plots without fish.

Each tab displays plots and models investigating the three metrics:

Porewater analyses:

Plotting full year

Plotting winter flooding period only

The following graphs display \(\delta\) C13 isotopic ratios by treatment when fish were on fields, shown in boxplots because of low sample size, particularly in year 1 where only two samples per plot were taken towards the end of the winter flooded period.

Year 1 Models:

Below is a table of model comparison results. All models are linear models investigating \(\delta\) 13C values in porewater by treatment and date. For year 1, we have only two sampling dates during the winter flooded period. Due to repeated sampling, we first tried a linear mixed effects model with plot ID as a random effect. However, this produced a singular fit (plot level variance was ~0). So, we dropped the random effect. Outliers were checked and removed from the dataset. Model diagnostics were checked with the DHARMa package. Estimated marginal means comparisons were done using the emmeans package.

Response Model Comparison F statistic P-value
Porewater \(\delta\) 13C M1: fish treatment * bird treatment M1 vs M2 0.182 0.675
M2: fish treatment + bird treatment M2 vs M3 0.049 0.827
M3: fish treatment M3 vs M4 5.096 0.037 *
M4: no fish NA NA NA

We see a marginal effect of fish treatment:

## 
## Call:
## lm(formula = C13PW ~ FishTreatment + Date, data = porewater_no_out)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9679 -0.8266 -0.3605  0.5541  2.2121 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)  
## (Intercept)         898.53321  758.12854   1.185   0.2514  
## FishTreatmentNoFish  -1.29940    0.57559  -2.258   0.0366 *
## Date                 -0.04826    0.03837  -1.258   0.2246  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.303 on 18 degrees of freedom
## Multiple R-squared:  0.2987, Adjusted R-squared:  0.2208 
## F-statistic: 3.833 on 2 and 18 DF,  p-value: 0.04104

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(M1)
## W = 0.94176, p-value = 0.2361

Plotted model results:

Jittered points display raw data, overlayed with estimated marginal means from model 1, with error bars representing upper and lower 95% confidence limits.

##  FishTreatment emmean    SE df lower.CL upper.CL
##  Fish           -54.9 0.416 18    -55.8    -54.1
##  NoFish         -56.2 0.394 18    -57.1    -55.4
## 
## Results are averaged over the levels of: Date 
## Confidence level used: 0.95

Year 2 Models:

Below is a table of model comparison results. All models are linear mixed effects models investigating \(\delta\) 13C values in porewater by treatment. Including date caused models to have convergence issues. For year 2, we included plot ID as a random effect to account for repeated sampling. Model diagnostics were checked with the DHARMa package. Estimated marginal means comparisons were done using the emmeans package.

Response Model Log-liklihood Comparison X2 P
Porewater \(\delta\) 13C M1: fish treatment * bird treatment -82.80 M1 vs M2 0.908 0.341
M2: fish treatment + bird treatment -83.25 M2 vs M3 1.673 0.196
M3: fish treatment -84.09 M3 vs M4 0.230 0.597
M4: no fish -84.23 NA NA NA

We see no evidence of treatment effect in year 2.

Headspace analyses:

Plotting full year

Plotting winter flooding period

Year 1 Models:

Below is a table of model comparison results. All models are linear mixed effects models investigating \(\delta\) 13C values in headspace gas by treatment and date. We included plot ID as a random effect to account for repeated sampling. Model diagnostics were checked with the DHARMa package. Estimated marginal means comparisons were done using the emmeans package.

Response Model Log-liklihood Comparison X2 P
Headspace \(\delta\) 13C M1: fish treatment * bird treatment -59.42 M1 vs M2 3.464 0.063 ·
M2: fish treatment + bird treatment -61.15 M2 vs M3 0.002 0.967
M3: fish treatment -61.16 M3 vs M4 2.826 0.093 ·
M4: no fish -62.57 NA NA NA

We see a marginal effect of the model that includes bird*fish interaction:

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: C13Gas ~ (BirdTreatment * FishTreatment) + Date + (1 | PlotID)
##    Data: isotope_yr1_fish
## 
## REML criterion at convergence: 107.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -0.98047 -0.57211 -0.09588  0.53130  1.18175 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  PlotID   (Intercept) 10.68    3.268   
##  Residual             12.94    3.597   
## Number of obs: 21, groups:  PlotID, 15
## 
## Fixed effects:
##                                           Estimate Std. Error        df t value
## (Intercept)                              5124.5597  2528.4550    6.4357   2.027
## BirdTreatmentNoBirds                        4.6487     3.6170   12.2573   1.285
## FishTreatmentNoFish                         9.2978     3.9859   14.9445   2.333
## Date                                       -0.2622     0.1280    6.4357  -2.048
## BirdTreatmentNoBirds:FishTreatmentNoFish   -8.3750     4.9909   12.4243  -1.678
##                                          Pr(>|t|)  
## (Intercept)                                0.0859 .
## BirdTreatmentNoBirds                       0.2225  
## FishTreatmentNoFish                        0.0341 *
## Date                                       0.0833 .
## BirdTreatmentNoBirds:FishTreatmentNoFish   0.1183  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) BrdTNB FshTNF Date  
## BrdTrtmntNB  0.265                     
## FshTrtmntNF  0.405  0.654              
## Date        -1.000 -0.265 -0.406       
## BrdTNB:FTNF -0.280 -0.748 -0.781  0.280

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(M1)
## W = 0.91743, p-value = 0.07712
##  FishTreatment BirdTreatment emmean   SE    df lower.CL upper.CL
##  Fish          Birds          -56.4 2.98 15.14    -62.7    -50.0
##  NoFish        Birds          -47.1 2.47 12.60    -52.4    -41.7
##  Fish          NoBirds        -51.7 2.07  7.31    -56.6    -46.9
##  NoFish        NoBirds        -50.8 2.37 11.70    -56.0    -45.6
## 
## Results are averaged over the levels of: Date 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

We also see a marginal effect of the model that includes fish treatment only:

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: C13Gas ~ FishTreatment + Date + (1 | PlotID)
##    Data: isotope_yr1_fish
## 
## REML criterion at convergence: 119.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.18170 -0.65964  0.07057  0.56541  1.25093 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  PlotID   (Intercept) 10.57    3.251   
##  Residual             13.95    3.735   
## Number of obs: 21, groups:  PlotID, 15
## 
## Fixed effects:
##                      Estimate Std. Error        df t value Pr(>|t|)
## (Intercept)         3857.2732  2499.6337    7.2541   1.543    0.165
## FishTreatmentNoFish    4.0274     2.4949   12.3814   1.614    0.132
## Date                  -0.1979     0.1265    7.2551  -1.564    0.160
## 
## Correlation of Fixed Effects:
##             (Intr) FshTNF
## FshTrtmntNF  0.305       
## Date        -1.000 -0.306

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(M2)
## W = 0.95673, p-value = 0.4529
##  FishTreatment emmean   SE   df lower.CL upper.CL
##  Fish           -53.2 1.76 11.1    -57.0    -49.3
##  NoFish         -49.1 1.76 14.4    -52.9    -45.4
## 
## Results are averaged over the levels of: Date 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

Year 2 Models:

Below is a table of model comparison results. All models are linear mixed effects models investigating \(\delta\) 13C values in headspace gas by treatment and date. For year 2, we included plot ID as a random effect to account for repeated sampling. Model diagnostics were checked with the DHARMa package. Estimated marginal means comparisons were done using the emmeans package.

Response Model Log-liklihood Comparison X2 P
Headspace \(\delta\) 13C M1: fish treatment * bird treatment -144.33 M1 vs M2 1.088 0.197
M2: fish treatment + bird treatment -144.88 M2 vs M3

3.819

R

 0.148
M3: fish treatment -146.79 M3 vs M4 0.463 0.496
M4: no fish -147.79 NA NA NA

We see no evidence of treatment effect in year 2.

𝚫 \(\delta\) C13 analyses:

Plotting full year

Plotting winter flooding period

Year 1 Models:

Below is a table of model comparison results. All models are linear mixed effects models investigating 𝛥\(\delta\) 13C values in headspace gas by treatment and date. We included plot ID as a random effect to account for repeated sampling. Model diagnostics were checked with the DHARMa package. Estimated marginal means comparisons were done using the emmeans package.

Response Model Log-liklihood Comparison X2 P
𝛥\(\delta\) 13C M1: fish treatment * bird treatment -65.872 M1 vs M2 1.57 0.210
M2: fish treatment + bird treatment -66.657 M2 vs M3 0.202 0.653
M3: fish treatment -66.753 M3 vs M4 3.165 0.075 ·
M4: no fish -68.758 NA NA NA

We see a marginal effect of fish treatment:

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: delta ~ FishTreatment + scale(Date) + (1 | PlotID)
##    Data: isotope_yr1_fish
## 
## REML criterion at convergence: 124.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9577 -0.5551  0.2666  0.6970  1.5415 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  PlotID   (Intercept)  3.521   1.876   
##  Residual             36.003   6.000   
## Number of obs: 21, groups:  PlotID, 15
## 
## Fixed effects:
##                     Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)           1.3404     2.0124 11.5358   0.666    0.518
## FishTreatmentNoFish   5.1006     2.9836 13.5720   1.710    0.110
## scale(Date)          -0.4902     1.4516 12.9363  -0.338    0.741
## 
## Correlation of Fixed Effects:
##             (Intr) FshTNF
## FshTrtmntNF -0.716       
## scale(Date)  0.238 -0.334

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(M1)
## W = 0.96787, p-value = 0.6856
##  FishTreatment emmean   SE    df lower.CL upper.CL
##  Fish            1.36 2.08  9.06    -3.35     6.07
##  NoFish          6.46 2.22 13.86     1.70    11.23
## 
## Results are averaged over the levels of: Date 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

Year 2 Models:

Below is a table of model comparison results. All models are linear mixed effects models investigating 𝛥\(\delta\) 13C values in headspace gas by treatment and date. We included plot ID as a random effect to account for repeated sampling. Model diagnostics were checked with the DHARMa package. Estimated marginal means comparisons were done using the emmeans package.

Response Model Log-liklihood Comparison X2 P
𝛥\(\delta\) 13C M1: fish treatment * bird treatment -162.72 M1 vs M2 0.297 0.586
M2: fish treatment + bird treatment -162.87 M2 vs M3 0.037 0.848
M3: fish treatment -162.88 M3 vs M4 0.046 0.831
M4: no fish -162.91 NA NA NA

We see no evidence of treatment effect in year 2.

GHG flux

CH4 Timeseries

Full year

Displaying entire year of methane data as collected by LI-COR trace gas analyzers for year 1 and 2 of study.

Winter flooding period only:

Displaying only when fish are on fields.

Year 1 Models

Model comparison: methane flux ~ treatment:

Below is a table comparing methane flux models, investigating methane flux by fish/bird treatment. For year 1, we don’t have enough data pre- fish introductions to compare pre/post. We begin by establishing there were no treatment differences in our pre-fish sampling date, using a simple ANOVA. Then, we utilized a linear mixed effects model that includes only the period post- fish introductions to investigate how fish affected methane emissions. Methane emissions were log transformed to ensure normal distribution of residuals. Plot ID and sampling occasion were included as random effects in all models. Model diagnostics were run with DHARMa and performance packages (residual diagnostics, normality, collinearity, heteroscedasity). Covariates (water temperature, PO4, NH4) were included if no collinearity was detected (VIF < 3). The emmeans package was used to compare effects.

ANOVA of pre- fish introduction methane flux by treatment:

anova(log(methane flux) ~ fish treatment * bird treatment)

##                             Df Sum Sq Mean Sq F value Pr(>F)
## FishTreatment                1  0.084  0.0845   0.041  0.843
## BirdTreatment                1  1.725  1.7254   0.835  0.380
## FishTreatment:BirdTreatment  1  0.698  0.6982   0.338  0.573
## Residuals                   11 22.726  2.0660

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(anova1)
## W = 0.93518, p-value = 0.3256

Model comparisons: post- fish introduction methane flux by treatment:

Response Model Log Likelihood Comparison X2 P
Methane flux fish treatment * bird treatment -82.08 M1 vs M2 1.87 0.171
fish treatment + bird treatment -83.02 M2 vs M3 0.01 0.925
fish treatment -83.02 M3 vs M4 9.75 0.002*
no fish -87.90 NA NA NA

Year 1: methane flux post fish introduction

lmer(log(methane flux) ~ fish treatment + water temperature + PO4 + NH4 + (1 | PlotID) + (1 | Sampling Occasion))

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logFlux ~ FishTreatment + scale(MeanTemp) + scale(PO4) + scale(NH4) +  
##     (1 | PlotID) + (1 | SamplingOccasion)
##    Data: GHG_Yr1_post
## 
## REML criterion at convergence: 169.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4388 -0.6833 -0.2333  0.4707  2.0156 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  PlotID           (Intercept) 0.4214   0.6491  
##  SamplingOccasion (Intercept) 0.1584   0.3979  
##  Residual                     2.0129   1.4188  
## Number of obs: 46, groups:  PlotID, 16; SamplingOccasion, 3
## 
## Fixed effects:
##                     Estimate Std. Error      df t value Pr(>|t|)   
## (Intercept)           2.6455     0.4475  2.8760   5.912   0.0109 * 
## FishTreatmentNoFish   1.8087     0.5313 10.7152   3.404   0.0061 **
## scale(MeanTemp)      -0.3016     0.3463  1.0819  -0.871   0.5348   
## scale(PO4)            0.3779     0.2557 39.4412   1.478   0.1473   
## scale(NH4)           -0.2375     0.2212 36.3136  -1.074   0.2900   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) FshTNF sc(MT) s(PO4)
## FshTrtmntNF -0.593                     
## scal(MnTmp)  0.179 -0.004              
## scale(PO4)  -0.022 -0.043 -0.108       
## scale(NH4)   0.015  0.017 -0.115 -0.182

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(M1)
## W = 0.95629, p-value = 0.08219
##  contrast      estimate    SE   df t.ratio p.value
##  Fish - NoFish    -1.81 0.532 12.9  -3.397  0.0048
## 
## Degrees-of-freedom method: kenward-roger

Plotted model results:

Jittered points display raw data, overlayed with estimated marginal means from model 1, with error bars representing upper and lower 95% confidence limits.

##  FishTreatment emmean    SE   df lower.CL upper.CL
##  Fish            2.75 0.441 3.36     1.43     4.07
##  NoFish          4.56 0.440 3.36     3.24     5.88
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

Year 2 Models

Model comparison: methane flux ~ treatment, pre/post fish introduction:

Below is a table comparing methane flux models, investigating methane flux by fish/bird treatment. Here, we utilized a linear mixed effects model including pre and post fish introductions to investigate how fish affected methane emissions. Methane emissions were log transformed to ensure normal distribution of residuals. Plot ID and sampling occasion were included as random effects in all models. Model diagnostics were run with DHARMa and performance packages (residual diagnostics, normality, collinearity, heteroscedasity). Covariates (redox, NO3, NH4) were included if no collinearity was detected (VIF < 3).

Response Model Log Likelihood Comparison X2 P
Methane flux fish treatment * bird treatment * pre/post -117.86 M1 vs M2 0.01 0.926
(fish treatment * pre/post) + (bird treatment * pre/post) + (fish treatment * bird treatment) -117.86 M2 vs M3 0.54 0.462
(fish treatment * pre/post) + (bird treatment * pre/post) -118.14 M3 vs M4 0.14 0.706
(fish treatment * pre/post) + bird treatment -118.21 M4 vs M5 0.04 0.838
fish treatment * pre/post -118.23 M5 vs M6 1.02 0.313
fish treatment + pre/post -118.74 M6 vs M7 4.60 0.032 *
fish treatment -121.04 M7 vs M8 1.58 0.208
no fish treatment or pre/post -121.83 NA NA NA

Model 1:

lmer(log(methane flux) ~ fish treatment + pre/post + redox + NO3 + NH4 + (1|Plot ID) + (1|Sampling Occassion)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logFlux ~ FishTreatment + Pre_Post + scale(Redox) + scale(NO3) +  
##     scale(NH4) + (1 | PlotID)
##    Data: GHG_Yr2_winter
## 
## REML criterion at convergence: 237.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.70580 -0.58679 -0.02087  0.42401  2.08772 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  PlotID   (Intercept) 1.934    1.391   
##  Residual             2.167    1.472   
## Number of obs: 61, groups:  PlotID, 16
## 
## Fixed effects:
##                     Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)          3.99785    0.71296 29.97984   5.607 4.21e-06 ***
## FishTreatmentNoFish -0.97885    0.80178 13.97763  -1.221   0.2423    
## Pre_PostPre          1.14063    0.45109 43.24618   2.529   0.0152 *  
## scale(Redox)        -0.04671    0.32890 49.06398  -0.142   0.8877    
## scale(NO3)          -1.04947    0.97472 45.77495  -1.077   0.2873    
## scale(NH4)           0.24190    0.18802 50.24321   1.287   0.2041    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) FshTNF Pr_PsP scl(R) s(NO3)
## FshTrtmntNF -0.574                            
## Pre_PostPre -0.174  0.021                     
## scale(Redx) -0.350 -0.044 -0.466              
## scale(NO3)   0.548 -0.042  0.131 -0.671       
## scale(NH4)  -0.092  0.076 -0.114  0.156 -0.026

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(M1)
## W = 0.97544, p-value = 0.2571
##  contrast      estimate    SE df t.ratio p.value
##  Fish - NoFish    0.979 0.802 14   1.220  0.2426
## 
## Results are averaged over the levels of: Pre_Post 
## Degrees-of-freedom method: kenward-roger
##  FishTreatment Pre_Post emmean    SE   df lower.CL upper.CL
##  Fish          Post       4.29 0.617 19.2     3.00     5.58
##  NoFish        Post       3.31 0.613 18.8     2.03     4.60
##  Fish          Pre        5.43 0.598 17.1     4.17     6.69
##  NoFish        Pre        4.45 0.606 17.9     3.18     5.73
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

CH4 Cumulative emissions

To investigate cumulative emissions, we will follow the summation method described in Echeverría-Progulakis et al 2025 J. Enviro. Manag.). This method calculates cummulative greenhouse gas emissions using the following equation:

\(GHG_c = \sum_{i = 1}^{n} (\frac{GHG_{fi} * 𝛥time_{i+1}} {100})\)

Where \(GHG_c\) is cumulative gas emissions in kg/ha over the entirety of a given season or time period. \(GHG_{fi}\) is gas flux in mg m-2 h-1 , and \(𝜟time_{i+1}\) is the difference in time (in hours) between the \(ith\) and the the \((i + 1)th\) sampling events. Our sampling method, using a LI-COR trace gas analyzer, gave us real-time flux measurements in nm m-2 s-1, so we start by converting this to the correct units using the following conversion, by first converting nmol to molar mass of methane (16.04 g mol-1), and then converting seconds to hours. This results in nmol m-2 s-1 x 16.04 x 10-6 = mg m-2 h-1. After converting our flux values to mg m-2 h-1, then the cumulative GHG units result in kg/ha.

Year 1 Summation: Full year

Panel A) shows cumulative emissions between each sampling date (approximating area under the curve flux) by treatment, and panel B) shows total seasonal annual methane flux by treatment.

Modeling the cumulative full year data (and instantaneous full year flux data) is difficult, likely because of seasonal variance heterogeneity (i.e., different variance structures during the winter flooding season and rice growing season). Below is shown model comparison for the total season annual methane flux - utilizing simple linear model structure lm(total seasonal flux ~ treatment). Sample size is low here, with N = 8 for each treatment set (fish treatment and bird treatment). There is a marginal effect of fish treatment on annual cumulative emissions.

Response Model Comparison F statistic P-value
Seasonal cumulative CH4 flux (kg ha-1) M1: fish treatment * bird treatment M1 vs M2 0.203 0.660
M2: fish treatment + bird treatment M2 vs M3 0.077 0.786
M3: fish treatment NA NA NA
M4: bird treatment M2 vs M4 3.840 0.072 · 
## 
## Call:
## lm(formula = season_total ~ FishTreatment, data = season_totals)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1879.2  -544.2   185.3   502.7  1593.7 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           1998.0      321.5   6.215 2.26e-05 ***
## FishTreatmentNoFish    921.8      454.7   2.028   0.0621 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 909.3 on 14 degrees of freedom
## Multiple R-squared:  0.227,  Adjusted R-squared:  0.1718 
## F-statistic: 4.111 on 1 and 14 DF,  p-value: 0.06209
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(M1)
## W = 0.96715, p-value = 0.7906

Year 1 Summation: Winter flooding period

Panel A) shows cumulative emissions between each sampling date (approximating area under the curve flux) by treatment, and panel B) shows total seasonal methane flux by treatment for the winter flooded period of year 1 only.

Below is model comparison for the total season annual methane flux - utilizing simple linear model structure lm(total seasonal flux ~ treatment). Here, our response variable (total seasonal flux) is log transformed to ensure normal distribution of residuals. Sample size is low, with N = 8 for each treatment set (fish treatment and bird treatment). There is no evidence of an effect of treatment on seasonal flux.

Response Model Comparison F statistic P-value
Seasonal cumulative CH4 flux (kg ha-1) M1: fish treatment * bird treatment M1 vs M2 2.410 0.147
M2: fish treatment + bird treatment M2 vs M3 2.635 0.129
M3: fish treatment NA NA NA
M4: bird treatment M2 vs M4 0.874 0.367

Year 2 Summation: Full year

Panel A) shows cumulative emissions between each sampling date (approximating area under the curve flux) by treatment, and panel B) shows total annual seasonal methane flux by treatment for year 2.

Below is model comparison for the total season annual methane flux - utilizing simple linear model structure lm(total seasonal flux ~ treatment). Sample size is low, with N = 8 for each treatment set (fish treatment and bird treatment). There is no evidence of an effect of treatment on seasonal flux.

Response Model Comparison F statistic P-value
Seasonal cumulative CH4 flux (kg ha-1) M1: fish treatment * bird treatment M1 vs M2 0.049 0.828
M2: fish treatment + bird treatment M2 vs M3 0.111 0.745
M3: fish treatment NA NA NA
M4: bird treatment M2 vs M4 0.030 0.865

Year 2 Summation: Winter flooded period

Panel A) shows cumulative emissions between each sampling date (approximating area under the curve flux) by treatment, and panel B) shows total annual seasonal methane flux by treatment for year 2.

Below is model comparison for the total season methane flux for the winter flooded period in study year 2- utilizing simple linear model structure lm(total seasonal flux ~ treatment). Sample size is low, with N = 8 for each treatment set (fish treatment and bird treatment). Here, our response variable (seasonal flux) was log transformed to ensure normal distribution of residuals. There is no evidence of an effect of treatment on seasonal flux.

Response Model Comparison F statistic P-value
Seasonal cumulative CH4 flux (kg ha-1) M1: fish treatment * bird treatment M1 vs M2 0.526 0.482
M2: fish treatment + bird treatment M2 vs M3 0.026 0.875
M3: fish treatment NA NA NA
M4: bird treatment M2 vs M4 2.082 0.173