The data

Sampling sites (Fig. 1)

Phloem samples were taken from 12 trees, three from each of four forest sites in southern Sweden: Hälleskogsbrännan (HSB; 59°54’35.5“N 16°08’42.2”E), Ecopark Öjesjön (OJ; 59°51’04.2“N 16°14’57.8”E), Knutby (KB; 59°56’46.1“N 18°18’00.4”E and Skyttorp (SK; 60°06’32.4“N 17°49’34.0”E). In addition, 48 adult male and 48 female beetles were sampled from each site.

HSB and OJ were hit by wildfire during summer 2014, while KB and SK were not.

Sampling scheme

Beetles were paired together (due to such small amount of tissue), resulting in 24 female and 24 male beetle samples from each site.

Five trees from each site were felled in the spring, prior to potential bark beetle infestation, and a single phloem sample was taken from each tree (‘t0’/pre-colonization phloem). Three trees from each site were then chosen for inclusion in this study based on successful beetle colonization.

Phloem samples were taken from each tree at two time points: first when the bark beetle larvae where in the first instar stage of development (‘t1’/ 1st instar phloem), and again when they were in the 3rd instar of larval development (‘t2’/3rd instar phloem). At each of these time points, ten phloem samples were taken from each tree: two from uncolonized phloem tissue, and eight from the edges of bark beetle galleries.

In total, 192 beetle samples and 252 phloem samples were taken for inclusion in this study. DNA was sampled from a body wash and from the guts from each of the beetle samples and from all phloem samples, and sequenced using PacBio.

Carbon (C), nitrogen (N), and phosphorus (P) concentrations were measured for each phloem sample, and molar C:N, C:P, and N:P ratios were calculated.

Data files

Statistical analyses conducted in publication can be reproduced using the following files:

  • metadata file: ‘tp_metadata_20211202.txt’
  • OTU / sample matrix: ‘tp_initial_otu_table.txt’
  • SINTAX (Edgar, 2016) taxonomy predictions for each OTU sequence from which the ITS2 region could be successfully extracted: ‘tp_sintax_results.txt’

Additional data files associated with this publication:

  • OTU representative sequences for the OTUs included in the initial OTU table (‘tp_initial_otu_table.txt’): ‘tp_otu_seqs.fasta’

  • output files from OTU clustering (‘scata3997_JUVgMBeZyA.zip’) generated via the SCATA pipeline (https://scata.mykopat.slu.se/), which contains:
    - unfiltered OTU representative sequences (‘all_clusters_scata3997.fas’)
    - unfiltered OTU table (‘all_tag_by_cluster_counts.txt’)

  • tagged primer sequences (sample primer tag combinations can be found in the metadata file)

  • raw sequence data




Analysis Part I: fungal communities

Overview of community sequencing dataset

Initial OTU table

# otutab_in is file 'tp_initial_otu_table.txt'
initial_otutab <- otutab_in

# get list of fungi only OTUs
sintax_results_fungi <-
    sintax_results[startsWith(sintax_results$cutoff_0.8, "d:Fungi"),]

# subset the initial_otutab to include only the fungal OTUs
fungi_otus <- row.names(sintax_results_fungi)

# filter initial_otutab
otutab <- initial_otutab[, names(initial_otutab) %in% fungi_otus]

There are 412 samples, 362 OTUs, and 15854 sequence reads in OTU table (SCATA output, samples with primer tag combinations not included in this study plus singleton OTUs filtered out).

fungal OTUs only

There are 412 samples, 290 OTUs, and 13397 sequence reads in the OTU table after filtering non-fungal OTUs out.

Filtered OTU table

Remove positive control sequence

One OTU (scata3997_1; 1084 reads) that corresponds to positive control DNA used in PCR, so it should be removed (DNA extracted from Agaricus bisporus was used as a positive control, and scata3997_1 matches A. bisporus)

# remove OTU_1 from otutab
otutab <- otutab[, names(otutab) != "scata3997_1"]

# remove samples with 0 reads left
otutab <- otutab[rowSums(otutab) > 0,]  

There are 401 samples, 289 OTUs, and 12313 sequence reads in the fungal community dataset after removing scata3997_1.

Remove samples with less than 10 reads

Some of the samples have sequencing depths of less than 10 reads:

Total samples with less than 10 reads: 70
Minimum sequencing depth: 1
Maximum sequencing depth: 156 (sample HBTPFB11)

# remove samples with less than 10 reads for downstream analysis
otutab_min10 <- subset(otutab, rowSums(otutab) > 9)

# remove OTUs that have zero reads after removing samples with less than 10 reads
otutab_min10 <- otutab_min10[, colSums(otutab_min10) > 0]

After removing samples with less than 10 reads, the community dataset used for downstream analysis includes:

Total reads: 11962
Total samples: 331
Total OTUs: 286




Fig. 2 nMDS beetle samples

nMDS results

    
    Call:
    metaMDS(comm = otutab_beetle_ra, distance = "bray", k = 2, trymax = 200) 
    
    global Multidimensional Scaling using monoMDS
    
    Data:     otutab_beetle_ra 
    Distance: bray 
    
    Dimensions: 2 
    Stress:     0.2405601 
    Stress type 1, weak ties
    No convergent solutions - best solution after 200 tries
    Scaling: centring, PC rotation, halfchange scaling 
    Species: expanded scores based on 'otutab_beetle_ra'

nMDS plot

plot saved as ‘tp_Fig2_nMDS_beetles.pdf’

<b>Figure 2.</b> Non-metric multidimensional scaling (nMDS) ordination showing fungal community differentiation in bark beetle samples associated with fire history of sampling locality.  Filled symbols represent samples taken from plots that experienced forest fire, and open symbols represent samples from unburnt sites.  Plot excludes samples with less than 10 sequence reads. There are 167 samples, 258 OTUs and 7155 sequence reads represented in the plot.  Plot stress = 0.241.

Figure 2. Non-metric multidimensional scaling (nMDS) ordination showing fungal community differentiation in bark beetle samples associated with fire history of sampling locality. Filled symbols represent samples taken from plots that experienced forest fire, and open symbols represent samples from unburnt sites. Plot excludes samples with less than 10 sequence reads. There are 167 samples, 258 OTUs and 7155 sequence reads represented in the plot. Plot stress = 0.241.

This is the same plot as above, just highlighting sex and body region:

PERMANOVA

Fire history (burnt vs unburnt)

The most conspicuous pattern for beetle samples is a difference based on burnt vs unburnt collection locality. This is supported by PERMANOVA. A test for homogeneity of group dispersion tests for a significant difference in dispersion, which can lead to a false positive PERMANOVA test. However, a clear distinction between samples from burnt and unburnt sites can be observed.

# PERMANOVA test for significant difference between beetle sample fungal communities based on fire history (burnt vs unburnt)
adonis2(
    otutab_beetle_ra ~ fire_hist,
    data = meta_seq_beetle,
    by = 'terms',
    method = 'bray'
)
    Permutation test for adonis under reduced model
    Terms added sequentially (first to last)
    Permutation: free
    Number of permutations: 999
    
    adonis2(formula = otutab_beetle_ra ~ fire_hist, data = meta_seq_beetle, method = "bray", by = "terms")
               Df SumOfSqs      R2      F Pr(>F)    
    fire_hist   1    4.473 0.07355 13.099  0.001 ***
    Residual  165   56.346 0.92645                  
    Total     166   60.819 1.00000                  
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ANOVA testing for difference in dispersion of groups
anova(betadisper(vegdist(otutab_beetle_ra, method = "bray"), meta_seq_beetle$fire_hist))
    Analysis of Variance Table
    
    Response: Distances
               Df  Sum Sq  Mean Sq F value      Pr(>F)    
    Groups      1 0.27326 0.273258  22.387 0.000004764 ***
    Residuals 165 2.01402 0.012206                        
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sex (female vs male) and body region (body vs gut)

Grouping of samples based on beetle sex or body region is not apparent in ordination. PERMANOVA indicates that a small amount of variation can be explained by sex and body region, but this is a negligible amount. PERMANOVA indicates that a slightly larger amount of variation can be attributed to beetle body region, however, this may be due to the significantly different group dispersions.

sex

# PERMANOVA test for significant difference between beetle sample fungal communities based on sex (female vs male)
perm <- how(nperm = 199)
setBlocks(perm) <- with(meta_seq_beetle, site)

adonis2(
    otutab_beetle_ra ~ sex,
    data = meta_seq_beetle,
    by = 'terms',
    method = 'bray',
    permutations = perm
)
    Permutation test for adonis under reduced model
    Terms added sequentially (first to last)
    Blocks:  with(meta_seq_beetle, site) 
    Permutation: free
    Number of permutations: 199
    
    adonis2(formula = otutab_beetle_ra ~ sex, data = meta_seq_beetle, permutations = perm, method = "bray", by = "terms")
              Df SumOfSqs      R2      F Pr(>F)  
    sex        1    0.580 0.00954 1.5885  0.025 *
    Residual 165   60.239 0.99046                
    Total    166   60.819 1.00000                
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ANOVA testing for difference in dispersion of groups
anova(betadisper(vegdist(otutab_beetle_ra, method = "bray"), meta_seq_beetle$sex))
    Analysis of Variance Table
    
    Response: Distances
               Df  Sum Sq  Mean Sq F value Pr(>F)
    Groups      1 0.02047 0.020472  1.6313 0.2033
    Residuals 165 2.07062 0.012549

body region

# PERMANOVA test for significant difference between beetle sample fungal communities based on body region (body vs gut)
perm <- how(nperm = 199)
setBlocks(perm) <- with(meta_seq_beetle, site)

adonis2(
    otutab_beetle_ra ~ region,
    data = meta_seq_beetle,
    by = 'terms',
    method = 'bray',
    permutations = perm
)
    Permutation test for adonis under reduced model
    Terms added sequentially (first to last)
    Blocks:  with(meta_seq_beetle, site) 
    Permutation: free
    Number of permutations: 199
    
    adonis2(formula = otutab_beetle_ra ~ region, data = meta_seq_beetle, permutations = perm, method = "bray", by = "terms")
              Df SumOfSqs      R2      F Pr(>F)   
    region     1    0.953 0.01567 2.6264  0.005 **
    Residual 165   59.866 0.98433                 
    Total    166   60.819 1.00000                 
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ANOVA testing for difference in dispersion of groups
anova(betadisper(vegdist(otutab_beetle_ra, method = "bray"), meta_seq_beetle$region))
    Analysis of Variance Table
    
    Response: Distances
               Df  Sum Sq  Mean Sq F value  Pr(>F)  
    Groups      1 0.04267 0.042667  3.5917 0.05982 .
    Residuals 165 1.96007 0.011879                  
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



Fig. S1 displaying body and gut samples from the same individuals

plot saved as ‘tp_FigS1_nMDS_beetles_fig2_with_lines.pdf’

<b>Figure S1.</b> Plot from Fig. 2 with lines added to connect body and gut samples that originate from the same individual beetle sample, demonstrating that there is not a pattern driven by the non-independence of the samples.

Figure S1. Plot from Fig. 2 with lines added to connect body and gut samples that originate from the same individual beetle sample, demonstrating that there is not a pattern driven by the non-independence of the samples.

IDs for beetle samples from which both body and gut samples are included in analysis:

    $beetle_individual
     [1] "HBTPF10" "HBTPF13" "HBTPF8"  "HBTPM4"  "KBTPF1"  "KBTPF10" "KBTPF11" "KBTPF13" "KBTPF14" "KBTPF16" "KBTPF2"  "KBTPF21" "KBTPF23" "KBTPF5"  "KBTPM1" 
    [16] "KBTPM23" "OJTPF14" "OJTPF16" "OJTPF20" "OJTPF21" "OJTPF22" "OJTPF24" "OJTPF25" "OJTPF26" "OJTPM27" "OJTPM30" "OJTPM31" "OJTPM32" "OJTPM33" "OJTPM35"
    [31] "SKTPF1"  "SKTPF13" "SKTPF16" "SKTPF18" "SKTPF20" "SKTPF21" "SKTPF22" "SKTPF23" "SKTPF5"  "SKTPF7"  "SKTPM1"  "SKTPM2"  "SKTPM3"  "SKTPM4"  "SKTPM8" 
    [46] "SKTPM9"



Fig. 3 - nMDS phloem samples

nMDS results

    
    Call:
    metaMDS(comm = otutab_phloem_ra, distance = "bray", k = 2, trymax = 200) 
    
    global Multidimensional Scaling using monoMDS
    
    Data:     otutab_phloem_ra 
    Distance: bray 
    
    Dimensions: 2 
    Stress:     0.1859736 
    Stress type 1, weak ties
    No convergent solutions - best solution after 200 tries
    Scaling: centring, PC rotation, halfchange scaling 
    Species: expanded scores based on 'otutab_phloem_ra'

nMDS plot

plot saved as ‘tp_Fig3_nMDS_phloem.pdf’

<b>Figure 3.</b> Non-metric multidimensional scaling (nMDS) ordination displaying differentiation of fungal communities in phloem samples based on fire history of sampling locality.  Filled symbols represent samples taken from plots that experienced forest fire, and open symbols represent samples from unburnt sites.  Plot excludes samples with less than 10 sequence reads.  There are 164 samples, 96 OTUs and 4807 sequence reads represented in the plot.  Plot stress = 0.186.

Figure 3. Non-metric multidimensional scaling (nMDS) ordination displaying differentiation of fungal communities in phloem samples based on fire history of sampling locality. Filled symbols represent samples taken from plots that experienced forest fire, and open symbols represent samples from unburnt sites. Plot excludes samples with less than 10 sequence reads. There are 164 samples, 96 OTUs and 4807 sequence reads represented in the plot. Plot stress = 0.186.

PERMANOVA

Fire history (burnt vs unburnt)
# PERMANOVA test for significant difference between phloem sample fungal communities based fire history (burnt vs unburnt)
adonis2(
    otutab_phloem_ra ~ fire_hist,
    data = meta_seq_phloem,
    by = 'terms',
    method = 'bray'
)
    Permutation test for adonis under reduced model
    Terms added sequentially (first to last)
    Permutation: free
    Number of permutations: 999
    
    adonis2(formula = otutab_phloem_ra ~ fire_hist, data = meta_seq_phloem, method = "bray", by = "terms")
               Df SumOfSqs      R2      F Pr(>F)    
    fire_hist   1    2.033 0.05439 9.3182  0.001 ***
    Residual  162   35.346 0.94561                  
    Total     163   37.379 1.00000                  
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ANOVA testing for difference in dispersion of groups
anova(betadisper(vegdist(otutab_phloem_ra, method = "bray"), meta_seq_phloem$fire_hist))
    Analysis of Variance Table
    
    Response: Distances
               Df Sum Sq   Mean Sq F value Pr(>F)
    Groups      1 0.0003 0.0003371  0.0149 0.9029
    Residuals 162 3.6593 0.0225883
Uncolonized phloem vs colonized first instar vs colonized 3rd instar
# PERMAVONA test of significance for samples grouped (into 3 groups) as 1) uncolonized samples (pre-col, 1st instar, and 3rd instar collapsed), 2) 1st instar colonized, and 3) 3rd instar colonized) 

# 'hm_grp' is grouping in Fig 3 heatmap and nMDS below
perm <- how(nperm = 199)
setBlocks(perm) <- with(meta_seq_phloem, site)

adonis2(
    otutab_phloem_ra ~ hm_grp,
    data = meta_seq_phloem,
    by = 'terms',
    method = 'bray',
    permutations = perm
)
    Permutation test for adonis under reduced model
    Terms added sequentially (first to last)
    Blocks:  with(meta_seq_phloem, site) 
    Permutation: free
    Number of permutations: 199
    
    adonis2(formula = otutab_phloem_ra ~ hm_grp, data = meta_seq_phloem, permutations = perm, method = "bray", by = "terms")
              Df SumOfSqs      R2      F Pr(>F)   
    hm_grp     5    7.211 0.19293 7.5539  0.005 **
    Residual 158   30.167 0.80707                 
    Total    163   37.379 1.00000                 
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise PERMAVONA

Uncolonized phloem vs colonized first instar, and colonized first instar vs colonized 3rd instar

test of significance for uncolonized samples (pre-colonization, 1st instar, and 3rd instar combined) vs 1st instar colonized samples, and for 1st instar vs 3rd instar colonized samples)

# PERMANOVA for uncolonized vs 1st instar colonized samples
perm <- how(nperm = 199)
setBlocks(perm) <- with(meta_seq_phloem_uncol_v_1st, site)

permanova_uv1 <- adonis2(
    otutab_phloem_uncol_vs_1st_ra ~ hm_grp,
    data = meta_seq_phloem_uncol_v_1st,
    by = 'terms',
    method = 'bray',
    permutations = perm
)
permanova_uv1
    Permutation test for adonis under reduced model
    Terms added sequentially (first to last)
    Blocks:  with(meta_seq_phloem_uncol_v_1st, site) 
    Permutation: free
    Number of permutations: 199
    
    adonis2(formula = otutab_phloem_uncol_vs_1st_ra ~ hm_grp, data = meta_seq_phloem_uncol_v_1st, permutations = perm, method = "bray", by = "terms")
             Df SumOfSqs      R2      F Pr(>F)   
    hm_grp    3   2.8119 0.14808 4.5772  0.005 **
    Residual 79  16.1776 0.85192                 
    Total    82  18.9895 1.00000                 
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
perm <- how(nperm = 199)
setBlocks(perm) <- with(meta_seq_phloem_1st_vs_3rd, site)

# PERMANOVA for 1st vs 3rd instar colonized samples
permanova_1v3 <- adonis2(
    otutab_phloem_1st_vs_3rd_ra ~ hm_grp,
    data = meta_seq_phloem_1st_vs_3rd,
    by = 'terms',
    method = 'bray',
    permutations = perm
)
permanova_1v3
    Permutation test for adonis under reduced model
    Terms added sequentially (first to last)
    Blocks:  with(meta_seq_phloem_1st_vs_3rd, site) 
    Permutation: free
    Number of permutations: 199
    
    adonis2(formula = otutab_phloem_1st_vs_3rd_ra ~ hm_grp, data = meta_seq_phloem_1st_vs_3rd, permutations = perm, method = "bray", by = "terms")
              Df SumOfSqs      R2      F Pr(>F)   
    hm_grp     3    4.963 0.15488 8.6137  0.005 **
    Residual 141   27.081 0.84512                 
    Total    144   32.044 1.00000                 
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Bonferroni correction of p-values for multiple comparisons (pairwise PERMANOVA tests)

    ADJUSTED P-VALUES:  0.01 0.01



Fig. 4 - nMDS all samples

In order to investigate the OTUs that may be the highest contributors to community structure, we collapsed the samples to into group based on the results above.

  • beetle samples were collapsed, as visual inspection of nMDS ordination, PERMANOVA and tests of homogeneity of dispersion suggest that sex (male/female) and region (body/gut) do not explain a large proportion of the observed community variation. We also chose to retain both samples (body and gut) from individual beetle samples when available, as this aspect does not appear to contribute to variation in community composition, either. nMDS does, however indicate a difference between samples from burnt and unburnt sites along axis 1, which is supported by PERMANOVA.

  • uncolonized phloem samples taken at pre-colonization, 1st instar, and 3rd instar sampling time points were collapsed into a single group because they are so few and appear to cluster together in ordination, and the most notable pattern revealed is a general shift along axis 1 as uncolonized – colonized 1st instar – colonized, 3rd instar, and burnt – unburnt along axis 2, a pattern supported by PERMANOVA. Fire history is also supported by PERMANOVA as an explanatory factor of variation in fungal communities.

Therefore, we analyzed the OTUs that were both present in 30% of samples and had an average relative read abundance of 1% or greater in any of the following eight groups:

beetles, burnt
beetles, unburnt

phloem, no larvae, burnt
phloem, no larvae, unburnt

phloem, 1st instar, burnt
phloem, 1st instar, unburnt

phloem, 3rd instar, burnt
phloem, 3rd instar, unburnt

It is important to note that the beetle and phloem samples were run on two separate sequencing plates, so this could contribute to the differences between them. However, nMDS indicates that the 3rd instar samples are most like the beetle samples, followed by 3rd instar, and then uncolonized samples, which makes sense biologically. For this reason, we believe that the signal here is reliable.

nMDS results

    
    Call:
    metaMDS(comm = otutab_min10_ra, distance = "bray", k = 2, trymax = 200) 
    
    global Multidimensional Scaling using monoMDS
    
    Data:     otutab_min10_ra 
    Distance: bray 
    
    Dimensions: 2 
    Stress:     0.2230584 
    Stress type 1, weak ties
    No convergent solutions - best solution after 200 tries
    Scaling: centring, PC rotation, halfchange scaling 
    Species: expanded scores based on 'otutab_min10_ra'

nMDS plot

plot saved as ‘tp_Fig4_nMDS_all_samples.pdf’

<b>Figure 4.</b> Non-metric multidimensional scaling (nMDS) ordination displaying fungal com-munity shifts associated with forest fire, sampling type, and substrate (phloem vs. beetle samples).  Filled symbols represent samples taken from plots that experienced forest fire, and open symbols represent samples from unburnt sites.  Plot excludes samples with less than 10 sequence reads.  There are 331 samples, 286 OTUs, and 11962sequence reads represented in the plot.  Plot stress = 0.223.

Figure 4. Non-metric multidimensional scaling (nMDS) ordination displaying fungal com-munity shifts associated with forest fire, sampling type, and substrate (phloem vs. beetle samples). Filled symbols represent samples taken from plots that experienced forest fire, and open symbols represent samples from unburnt sites. Plot excludes samples with less than 10 sequence reads. There are 331 samples, 286 OTUs, and 11962sequence reads represented in the plot. Plot stress = 0.223.

This is the same plot, just simplified to show 1) only burnt vs unburnt, and 2) only beetle vs. phloem:

PERMANOVA

Phloem vs beetle samples
# PERMAVONA test of significant difference between fungal communities in phloem vs beetle samples (from all four sites)
perm <- how(nperm = 199)
setBlocks(perm) <- with(meta_seq_min10, site)

adonis2(
    otutab_min10_ra ~ seq_plate,
    data = meta_seq_min10,
    by = 'terms',
    method = 'bray',
    permutations = perm
)
    Permutation test for adonis under reduced model
    Terms added sequentially (first to last)
    Blocks:  with(meta_seq_min10, site) 
    Permutation: free
    Number of permutations: 199
    
    adonis2(formula = otutab_min10_ra ~ seq_plate, data = meta_seq_min10, permutations = perm, method = "bray", by = "terms")
               Df SumOfSqs      R2      F Pr(>F)   
    seq_plate   1   16.327 0.14256 54.701  0.005 **
    Residual  329   98.198 0.85744                 
    Total     330  114.525 1.00000                 
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Fire history (burnt vs unburnt)
# PERMAVONA test of significant difference between fungal communities in burnt vs unburnt sites (from all fungal communities, samples from both phloem and beetles)
adonis2(
    otutab_min10_ra ~ fire_hist,
    data = meta_seq_min10,
    by = 'terms',
    method = 'bray'
)
    Permutation test for adonis under reduced model
    Terms added sequentially (first to last)
    Permutation: free
    Number of permutations: 999
    
    adonis2(formula = otutab_min10_ra ~ fire_hist, data = meta_seq_min10, method = "bray", by = "terms")
               Df SumOfSqs      R2     F Pr(>F)    
    fire_hist   1    6.137 0.05359 18.63  0.001 ***
    Residual  329  108.388 0.94641                 
    Total     330  114.525 1.00000                 
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Phloem vs beetle samples plus fire history (marginal variation)
# marginal PERMANOVA for the two terms (amount of variation explained by each term after other term(s) accounted for)
adonis2(
    otutab_min10_ra ~ fire_hist + seq_plate,
    data = meta_seq_min10,
    by = 'margin',
    method = 'bray'
)
    Permutation test for adonis under reduced model
    Marginal effects of terms
    Permutation: free
    Number of permutations: 999
    
    adonis2(formula = otutab_min10_ra ~ fire_hist + seq_plate, data = meta_seq_min10, method = "bray", by = "margin")
               Df SumOfSqs      R2      F Pr(>F)    
    fire_hist   1    5.128 0.04478 18.072  0.001 ***
    seq_plate   1   15.317 0.13375 53.982  0.001 ***
    Residual  328   93.070 0.81266                  
    Total     330  114.525 1.00000                  
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

pairwise PERMANOVA

Fire history for uncolonized, colonized first, and colonized third separately

pairwise PERMANOVA test of significance between burnt and unburnt phloem samples as grouped in heatmap below

# pairwise PERMANOVA, uncolonized burnt vs uncolonized unburnt

permanova_nolarvae <- adonis2(
    otutab_phloem_nolarvae_ra ~ fire_hist,
    data = meta_seq_phloem_nolarvae,
    by = 'terms',
    method = 'bray')

permanova_nolarvae
    Permutation test for adonis under reduced model
    Terms added sequentially (first to last)
    Permutation: free
    Number of permutations: 999
    
    adonis2(formula = otutab_phloem_nolarvae_ra ~ fire_hist, data = meta_seq_phloem_nolarvae, method = "bray", by = "terms")
              Df SumOfSqs      R2     F Pr(>F)
    fire_hist  1   0.1392 0.04317 0.767  0.675
    Residual  17   3.0863 0.95683             
    Total     18   3.2255 1.00000
# pairwise PERMANOVA, 1st instar burnt vs 1st instar unburnt
meta_seq_phloem_1 <- subset(meta_seq_phloem,
                            substrate_collapsed_4 == "phloem, 1st instar")
meta_seq_phloem_1 <- droplevels(meta_seq_phloem_1)
otutab_phloem_1 <- otutab_phloem[row.names(meta_seq_phloem_1),]
otutab_phloem_1 <- otutab_phloem_1[, colSums(otutab_phloem_1) > 0]
otutab_phloem_1_ra <- decostand(otutab_phloem_1, method = "total")

permanova_1st <- adonis2(
    otutab_phloem_1_ra ~ fire_hist,
    data = meta_seq_phloem_1,
    by = 'terms',
    method = 'bray')

permanova_1st
    Permutation test for adonis under reduced model
    Terms added sequentially (first to last)
    Permutation: free
    Number of permutations: 999
    
    adonis2(formula = otutab_phloem_1_ra ~ fire_hist, data = meta_seq_phloem_1, method = "bray", by = "terms")
              Df SumOfSqs      R2      F Pr(>F)    
    fire_hist  1   1.4026 0.09677 6.6425  0.001 ***
    Residual  62  13.0913 0.90323                  
    Total     63  14.4939 1.00000                  
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# pairwise PERMANOVA, 3rd instar burnt vs 3rd instar unburnt
meta_seq_phloem_3 <- subset(meta_seq_phloem,
                            substrate_collapsed_4 == "phloem, 3rd instar")
meta_seq_phloem_3 <- droplevels(meta_seq_phloem_3)
otutab_phloem_3 <- otutab_phloem[row.names(meta_seq_phloem_3),]
otutab_phloem_3 <- otutab_phloem_3[, colSums(otutab_phloem_3) > 0]
otutab_phloem_3_ra <- decostand(otutab_phloem_3, method = "total")

permanova_3rd <- adonis2(
    otutab_phloem_3_ra ~ fire_hist,
    data = meta_seq_phloem_3,
    by = 'terms',
    method = 'bray')

permanova_3rd
    Permutation test for adonis under reduced model
    Terms added sequentially (first to last)
    Permutation: free
    Number of permutations: 999
    
    adonis2(formula = otutab_phloem_3_ra ~ fire_hist, data = meta_seq_phloem_3, method = "bray", by = "terms")
              Df SumOfSqs      R2      F Pr(>F)    
    fire_hist  1   1.8826 0.11861 10.631  0.001 ***
    Residual  79  13.9898 0.88139                  
    Total     80  15.8724 1.00000                  
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Bonferroni correction of p-values for multiple comparisons

    ADJUSTED P-VALUES:  1 0.003 0.003

Fig. 5 - heatmap of core community OTUs

These are the core OTUs:

     [1] "scata3997_0"  "scata3997_6"  "scata3997_5"  "scata3997_7"  "scata3997_11" "scata3997_8"  "scata3997_28" "scata3997_36" "scata3997_2"  "scata3997_3" 
    [11] "scata3997_13" "scata3997_22" "scata3997_15" "scata3997_17" "scata3997_19"

These OTUs account for 72.86 percent of the filtered OTU table (otutab_min10).

heatmap plot

plot saved as ‘tp_Fig5_heatmap.pdf’

<b>Figure 5.</b> Heatmap of the average relative sequence read abundance across samples for 15 OTUs identified as the core community in this study.  To the right, OTU number and assigned taxonomy.

Figure 5. Heatmap of the average relative sequence read abundance across samples for 15 OTUs identified as the core community in this study. To the right, OTU number and assigned taxonomy.

Number of samples included in each category:

site tree substrate n
KB K2 1st instar, uncolonized 1
KB K2 1st instar, colonized 5
KB K2 3rd instar, uncolonized 2
KB K2 3rd instar, colonized 5
KB K3 pre-colonization 1
KB K3 1st instar, uncolonized 1
KB K3 1st instar, colonized 6
KB K3 3rd instar, uncolonized 2
KB K3 3rd instar, colonized 6
KB K5 pre-colonization 1
KB K5 1st instar, uncolonized 1
KB K5 1st instar, colonized 2
KB K5 3rd instar, colonized 7
KB female, body 21
KB female, gut 12
KB male, body 22
KB male, gut 3
SK S2 1st instar, colonized 4
SK S2 3rd instar, colonized 8
SK S3 1st instar, colonized 6
SK S3 3rd instar, uncolonized 1
SK S3 3rd instar, colonized 7
SK S5 1st instar, colonized 7
SK S5 3rd instar, colonized 7
SK female, body 21
SK female, gut 11
SK male, body 12
SK male, gut 10
HSB H1 1st instar, colonized 8
HSB H1 3rd instar, uncolonized 1
HSB H1 3rd instar, colonized 6
HSB H2 1st instar, colonized 7
HSB H2 3rd instar, colonized 8
HSB H3 1st instar, uncolonized 1
HSB H3 1st instar, colonized 5
HSB H3 3rd instar, colonized 6
HSB female, body 7
HSB female, gut 3
HSB male, body 5
HSB male, gut 2
OJ O1 1st instar, uncolonized 1
OJ O1 1st instar, colonized 1
OJ O1 3rd instar, uncolonized 2
OJ O1 3rd instar, colonized 7
OJ O2 1st instar, colonized 5
OJ O2 3rd instar, uncolonized 1
OJ O2 3rd instar, colonized 6
OJ O4 1st instar, uncolonized 1
OJ O4 1st instar, colonized 8
OJ O4 3rd instar, uncolonized 2
OJ O4 3rd instar, colonized 8
OJ female, body 12
OJ female, gut 9
OJ male, body 10
OJ male, gut 7



Analysis Part II: phloem chemistry

C, N, and P concentrations were measured for each phloem sample. Three samples—SK5TPWC1, SK3TPWC3, and K3TPWC3—were excluded from analysis because either N or C was measured to be zero. (N and C were measured simultaneously for each sample, so the sample had to be removed entirely from analysis.) For the remaining 249 samples, molar ratios (C:N, C:P, N:P) were calculated.

Quick look at the phloem chemistry data

In general, more tree-to-tree variation is visible than site-to-site in this quick overview of the data.

C:N

C:P

N:P

Burnt vs unburnt C:N, C:P and N:P pre-colonization

Here are the pre-col values again:

<br><br><br>




One-way ANOVA tests were used to test for difference in nutrient ratios in trees from burnt vs unburnt sites prior to bark beetle colonization, with site treated as a random effects factor.

Based on the results, we fail to reject \(H_0\) that there is no difference in C:P, C:N and N:P ratios of the pre-colonization samples between the burnt and the unburnt sites. These results should be interpreted with caution, however, as our sample size may result in low power of detection.

ANOVA, burnt vs unburnt C:N pre-colonization

# ANOVA, burnt vs unburnt C:N pre-colonization
aov(cn_ratio ~ fire_hist + Error(site),
    data = subset(meta_chem, sampling_time == "t0"))
    
    Call:
    aov(formula = cn_ratio ~ fire_hist + Error(site), data = subset(meta_chem, 
        sampling_time == "t0"))
    
    Grand Mean: 95.71526
    
    Stratum 1: site
    
    Terms:
                    fire_hist Residuals
    Sum of Squares   133.5371  215.4130
    Deg. of Freedom         1         2
    
    Residual standard error: 10.37817
    Estimated effects are balanced
    
    Stratum 2: Within
    
    Terms:
                    Residuals
    Sum of Squares   892.5092
    Deg. of Freedom         8
    
    Residual standard error: 10.56237
summary(aov(cn_ratio ~ fire_hist + Error(site),
            data = subset(meta_chem, sampling_time == "t0")))
    
    Error: site
              Df Sum Sq Mean Sq F value Pr(>F)
    fire_hist  1  133.5   133.5    1.24  0.381
    Residuals  2  215.4   107.7               
    
    Error: Within
              Df Sum Sq Mean Sq F value Pr(>F)
    Residuals  8  892.5   111.6

ANOVA, burnt vs unburnt C:P pre-colonization

# ANOVA, burnt vs unburnt C:N pre-colonization
aov(cp_ratio ~ fire_hist + Error(site),
    data = subset(meta_chem, sampling_time == "t0"))
    
    Call:
    aov(formula = cp_ratio ~ fire_hist + Error(site), data = subset(meta_chem, 
        sampling_time == "t0"))
    
    Grand Mean: 1993.308
    
    Stratum 1: site
    
    Terms:
                    fire_hist Residuals
    Sum of Squares  1703021.7  945897.1
    Deg. of Freedom         1         2
    
    Residual standard error: 687.7125
    Estimated effects are balanced
    
    Stratum 2: Within
    
    Terms:
                    Residuals
    Sum of Squares    2262637
    Deg. of Freedom         8
    
    Residual standard error: 531.8172
summary(aov(cp_ratio ~ fire_hist + Error(site),
            data = subset(meta_chem, sampling_time == "t0")))
    
    Error: site
              Df  Sum Sq Mean Sq F value Pr(>F)
    fire_hist  1 1703022 1703022   3.601  0.198
    Residuals  2  945897  472949               
    
    Error: Within
              Df  Sum Sq Mean Sq F value Pr(>F)
    Residuals  8 2262637  282830

ANOVA, burnt vs unburnt N:P pre-colonization

# ANOVA, burnt vs unburnt C:N pre-colonization
aov(np_ratio ~ fire_hist + Error(site),
    data = subset(meta_chem, sampling_time == "t0"))
    
    Call:
    aov(formula = np_ratio ~ fire_hist + Error(site), data = subset(meta_chem, 
        sampling_time == "t0"))
    
    Grand Mean: 20.81932
    
    Stratum 1: site
    
    Terms:
                    fire_hist Residuals
    Sum of Squares   135.1451  118.5604
    Deg. of Freedom         1         2
    
    Residual standard error: 7.699363
    Estimated effects are balanced
    
    Stratum 2: Within
    
    Terms:
                    Residuals
    Sum of Squares   243.4461
    Deg. of Freedom         8
    
    Residual standard error: 5.516409
summary(aov(np_ratio ~ fire_hist + Error(site),
            data = subset(meta_chem, sampling_time == "t0")))
    
    Error: site
              Df Sum Sq Mean Sq F value Pr(>F)
    fire_hist  1  135.2  135.15    2.28   0.27
    Residuals  2  118.6   59.28               
    
    Error: Within
              Df Sum Sq Mean Sq F value Pr(>F)
    Residuals  8  243.4   30.43



Fig. 6 - nutrient ratio values for all samples

plot saved as ‘tp_Fig6_chem_by_fire_history_tree_means.pdf’

<b>Figure 6.</b> C:N, C:P and N:P ratios (molar) for all samples.  Lines indicate mean ratio value for colonized (solid line) and uncolonized (dashed line) phloem samples for each tree.  Faint line connects the pre-colonization value to both the colonized (closed symbols) and uncolonized (open symbols) value for each tree.  See Fig. S2 to examine values for samples from each tree separately, and Fig. S3 for model-predicted values plotted.

Figure 6. C:N, C:P and N:P ratios (molar) for all samples. Lines indicate mean ratio value for colonized (solid line) and uncolonized (dashed line) phloem samples for each tree. Faint line connects the pre-colonization value to both the colonized (closed symbols) and uncolonized (open symbols) value for each tree. See Fig. S2 to examine values for samples from each tree separately, and Fig. S3 for model-predicted values plotted.




Fig. S2 - nutrient ratio values for all samples organized by tree within burnt vs unburnt sites

plot saved as ‘tp_FigS2_chem_by_tree.pdf’

<b>Figure S2.</b> C:N, C:P, and N:P ratios for each tree separately.  See Fig. 6.

Figure S2. C:N, C:P, and N:P ratios for each tree separately. See Fig. 6.




Linear mixed-effects models, partially Bayesian implementation

Here, the relationships between each nutrient ratio (cn_ratio, cp_ratio, and np_ratio) and sampling time (sampling_time; 1st vs 3rd instar), fire history (fire_hist; burnt vs unburnt), and beetle infestation (bb_colonization; colonized vs uncolonized phloem) are examined. The nutrient ratio measured for each tree pre-colonization are treated as a nuisance covariate.

# center (but not scale) pre-col values and code fixed effects factors as sum-to-zero contrasts
meta_chem_t1t2_cent <- meta_chem_t1t2
contrasts(meta_chem_t1t2_cent$sampling_time) <- contr.sum
contrasts(meta_chem_t1t2_cent$bb_colonization) <- contr.sum
contrasts(meta_chem_t1t2_cent$fire_hist) <- contr.sum
meta_chem_t1t2_cent$precol_cn_c <- scale(meta_chem_t1t2_cent$precol_cn, scale = FALSE)
meta_chem_t1t2_cent$precol_cp_c <- scale(meta_chem_t1t2_cent$precol_cp, scale = FALSE)
meta_chem_t1t2_cent$precol_np_c <- scale(meta_chem_t1t2_cent$precol_np, scale = FALSE)
# test the random effects structure using intercept-only model

# interaction term [(1|site:tree)]
cn_fit_err1 <- blme::blmer(cn_ratio ~ 1 + (1|site:tree), data = meta_chem_t1t2_cent)
cp_fit_err1 <- blme::blmer(cp_ratio ~ 1 + (1|site:tree), data = meta_chem_t1t2_cent)
np_fit_err1 <- blme::blmer(np_ratio ~ 1 + (1|site:tree), data = meta_chem_t1t2_cent)

# nested term [(1|site/tree)]
# equivalent in this case to [(1|site) + (1|tree)]
cn_fit_err2 <- blme::blmer(cn_ratio ~ 1 + (1|site/tree), data = meta_chem_t1t2_cent)
cp_fit_err2 <- blme::blmer(cp_ratio ~ 1 + (1|site/tree), data = meta_chem_t1t2_cent)
np_fit_err2 <- blme::blmer(np_ratio ~ 1 + (1|site/tree), data = meta_chem_t1t2_cent) # this model fails to converge
    Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.600704 (tol = 0.002, component 1)
# compare models
anova(cn_fit_err1, cn_fit_err2, refit = FALSE)
    Data: meta_chem_t1t2_cent
    Models:
    cn_fit_err1: cn_ratio ~ 1 + (1 | site:tree)
    cn_fit_err2: cn_ratio ~ 1 + (1 | site/tree)
                npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
    cn_fit_err1    3 2024.7 2035.1 -1009.4   2018.7                    
    cn_fit_err2    4 2029.0 2042.8 -1010.5   2021.0     0  1          1
anova(cp_fit_err1, cp_fit_err2, refit = FALSE)
    Data: meta_chem_t1t2_cent
    Models:
    cp_fit_err1: cp_ratio ~ 1 + (1 | site:tree)
    cp_fit_err2: cp_ratio ~ 1 + (1 | site/tree)
                npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
    cp_fit_err1    3 3467.1 3477.5 -1730.5   3461.1                     
    cp_fit_err2    4 3468.7 3482.6 -1730.4   3460.7 0.3301  1     0.5656
anova(np_fit_err1, np_fit_err2, refit = FALSE)
    Data: meta_chem_t1t2_cent
    Models:
    np_fit_err1: np_ratio ~ 1 + (1 | site:tree)
    np_fit_err2: np_ratio ~ 1 + (1 | site/tree)
                npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
    np_fit_err1    3 1183.1 1193.5 -588.57   1177.1                    
    np_fit_err2    4 1185.1 1199.0 -588.56   1177.1 0.011  1     0.9166

fit the model of interest and test the assumption of equal slopes

In order to treat precol_xx as a nuisance covariate, we must first check that the assumption of equal slopes is met by comparing the models above to models that include an interaction between the pre-colonization value and each of the three explanatory variables (fire_hist, bb_colonization, and sampling_time).

# fit a model with variables of interest plus nuisance covariate pre-col value
cn_fit_nus <- blme::blmer(cn_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cn_c +
                              (1|site:tree), data = meta_chem_t1t2_cent)

cp_fit_nus <- blme::blmer(cp_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cp_c +
                              (1|site:tree), data = meta_chem_t1t2_cent)

np_fit_nus <- blme::blmer(np_ratio ~ fire_hist * bb_colonization * sampling_time + precol_np_c +
                              (1|site:tree), data = meta_chem_t1t2_cent)

# fit a model that includes an interaction between the covariate and each of the three fixed effects variables of interest
# see if the model with the interactions is significantly better to test assumption of equal slopes
cn_fit_int <- blme::blmer(cn_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cn_c +
                              precol_cn_c:fire_hist + precol_cn_c:bb_colonization + precol_cn_c:sampling_time +
                              (1|site:tree), data = meta_chem_t1t2_cent)

cp_fit_int <- blme::blmer(cp_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cp_c +
                              precol_cp_c:fire_hist + precol_cp_c:bb_colonization + precol_cp_c:sampling_time +
                              (1|site:tree), data = meta_chem_t1t2_cent)

np_fit_int <- blme::blmer(np_ratio ~ fire_hist * bb_colonization * sampling_time + precol_np_c +
                              precol_np_c:fire_hist + precol_np_c:bb_colonization + precol_np_c:sampling_time +
                              (1|site:tree), data = meta_chem_t1t2_cent)

anova(cn_fit_nus, cn_fit_int)
    refitting model(s) with ML (instead of REML)
    Data: meta_chem_t1t2_cent
    Models:
    cn_fit_nus: cn_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cn_c + (1 | site:tree)
    cn_fit_int: cn_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cn_c + precol_cn_c:fire_hist + precol_cn_c:bb_colonization + precol_cn_c:sampling_time + (1 | site:tree)
               npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
    cn_fit_nus   11 1853.0 1891.1 -915.48   1831.0                       
    cn_fit_int   14 1849.2 1897.8 -910.61   1821.2 9.7412  3     0.0209 *
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(cp_fit_nus, cp_fit_int)
    refitting model(s) with ML (instead of REML)
    Data: meta_chem_t1t2_cent
    Models:
    cp_fit_nus: cp_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cp_c + (1 | site:tree)
    cp_fit_int: cp_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cp_c + precol_cp_c:fire_hist + precol_cp_c:bb_colonization + precol_cp_c:sampling_time + (1 | site:tree)
               npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
    cp_fit_nus   11 3378.7 3416.8 -1678.3   3356.7                     
    cp_fit_int   14 3378.5 3427.0 -1675.2   3350.5 6.2103  3     0.1018
anova(np_fit_nus, np_fit_int)
    refitting model(s) with ML (instead of REML)
    Data: meta_chem_t1t2_cent
    Models:
    np_fit_nus: np_ratio ~ fire_hist * bb_colonization * sampling_time + precol_np_c + (1 | site:tree)
    np_fit_int: np_ratio ~ fire_hist * bb_colonization * sampling_time + precol_np_c + precol_np_c:fire_hist + precol_np_c:bb_colonization + precol_np_c:sampling_time + (1 | site:tree)
               npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
    np_fit_nus   11 1097.4 1135.6 -537.71   1075.4                         
    np_fit_int   14 1085.5 1134.0 -528.73   1057.5 17.974  3  0.0004454 ***
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The model with the interactions is significantly better in the model comparisons for C:N and N:P, indicating that the ANCOVA assumption of equal slopes is violated for these. For the C:P model, both alternatives have nearly the same AIC, but with different number of parameters. Therefore, a mixed-effects model is retained with the interactions included in all three cases. However, it seems that it is mostly driven by the interaction between precol_xx and sampling_time, so models are fitted with only the precol_xx:sampling_time interaction included to see if the model with all three interactions is still significantly better.

# fit a model with only the precol_xx:sampling_time interaction and see if the model with all three interactions is still significantly better
cn_fit_int_red <- blme::blmer(cn_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cn_c +
                                  precol_cn_c:sampling_time + (1|site:tree), data = meta_chem_t1t2_cent)

cp_fit_int_red <- blme::blmer(cp_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cp_c +
                                  precol_cp_c:sampling_time + (1|site:tree), data = meta_chem_t1t2_cent)

np_fit_int_red <- blme::blmer(np_ratio ~ fire_hist * bb_colonization * sampling_time + precol_np_c +
                                  precol_np_c:sampling_time + (1|site:tree), data = meta_chem_t1t2_cent)

anova(cn_fit_int, cn_fit_int_red)
    refitting model(s) with ML (instead of REML)
    Data: meta_chem_t1t2_cent
    Models:
    cn_fit_int_red: cn_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cn_c + precol_cn_c:sampling_time + (1 | site:tree)
    cn_fit_int: cn_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cn_c + precol_cn_c:fire_hist + precol_cn_c:bb_colonization + precol_cn_c:sampling_time + (1 | site:tree)
                   npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
    cn_fit_int_red   12 1845.5 1887.1 -910.76   1821.5                     
    cn_fit_int       14 1849.2 1897.8 -910.61   1821.2 0.3075  2     0.8575
anova(cp_fit_int, cp_fit_int_red)
    refitting model(s) with ML (instead of REML)
    Data: meta_chem_t1t2_cent
    Models:
    cp_fit_int_red: cp_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cp_c + precol_cp_c:sampling_time + (1 | site:tree)
    cp_fit_int: cp_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cp_c + precol_cp_c:fire_hist + precol_cp_c:bb_colonization + precol_cp_c:sampling_time + (1 | site:tree)
                   npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
    cp_fit_int_red   12 3380.3 3421.9 -1678.2   3356.3                       
    cp_fit_int       14 3378.5 3427.0 -1675.2   3350.5 5.8374  2      0.054 .
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(np_fit_int, np_fit_int_red)
    refitting model(s) with ML (instead of REML)
    Data: meta_chem_t1t2_cent
    Models:
    np_fit_int_red: np_ratio ~ fire_hist * bb_colonization * sampling_time + precol_np_c + precol_np_c:sampling_time + (1 | site:tree)
    np_fit_int: np_ratio ~ fire_hist * bb_colonization * sampling_time + precol_np_c + precol_np_c:fire_hist + precol_np_c:bb_colonization + precol_np_c:sampling_time + (1 | site:tree)
                   npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
    np_fit_int_red   12 1083.9 1125.5 -529.93   1059.9                     
    np_fit_int       14 1085.5 1134.0 -528.73   1057.5 2.4065  2     0.3002
# assign the final models (model assumptions tested below)
cn_fit <- cn_fit_int_red
cp_fit <- cp_fit_int_red
np_fit <- np_fit_int_red



Fig. S3 - model-predicted values

plot saved as ‘tp_FigS3_chem_by_fire_history_predicted_values.pdf’

<b>Figure S3.</b> Predicted C:N, C:P and N:P ratios based on partially Bayesian linear mixed-effects models.  Measured C:N, C:P and N:P ratio values for each sample are plotted in Fig. 6 and Fig. S2.

Figure S3. Predicted C:N, C:P and N:P ratios based on partially Bayesian linear mixed-effects models. Measured C:N, C:P and N:P ratio values for each sample are plotted in Fig. 6 and Fig. S2.




Comparison of model-predicted vs mean nutrient ratio values (subsets of Fig. 5 and Fig. S3)




Checking model assumptions

C:N

1. linearity of predictors

visual inspection, plot of residuals against response values:

2. homogeneity of variance

visual inspection, plot of residuals against fitted values:

statistical test, Levene’s test:

The assumption of equal variance of residuals is violated for several grouping factors, however this should not be a problem, since large sample sizes almost always result in a significant difference. The robustness of LMERs against violations like this has also been documented (Schielzeth et al 2020, DOI: 10.1111/2041-210X.13434).

leveneTest(residuals(cn_fit) ~ meta_w_preds_and_resids$tree)
    Levene's Test for Homogeneity of Variance (center = median)
           Df F value     Pr(>F)    
    group  11  3.9057 0.00003454 ***
          225                       
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
leveneTest(residuals(cn_fit) ~ meta_w_preds_and_resids$bb_colonization)
    Levene's Test for Homogeneity of Variance (center = median)
           Df F value Pr(>F)
    group   1  0.7226 0.3962
          235
leveneTest(residuals(cn_fit) ~ meta_w_preds_and_resids$sampling_time)
    Levene's Test for Homogeneity of Variance (center = median)
           Df F value  Pr(>F)  
    group   1  6.7237 0.01011 *
          235                  
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
leveneTest(residuals(cn_fit) ~ meta_w_preds_and_resids$site)
    Levene's Test for Homogeneity of Variance (center = median)
           Df F value  Pr(>F)  
    group   3  2.7622 0.04282 *
          233                  
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3. normality of errors

visual inspection, histogram of residuals:

visual inspection, plot of residuals for each explanatory factor:

visual inspection, plot of measured value against residuals:

statistical test, Shapiro-Wilk test:

As with, homogeneity of variance, assumption violations occur here, but this should not be a problem for the same reasons as homogeneity of variance.

shapiro.test(meta_w_preds_and_resids$residuals_cn)
    
        Shapiro-Wilk normality test
    
    data:  meta_w_preds_and_resids$residuals_cn
    W = 0.98414, p-value = 0.009555
shapiro.test(scale(meta_w_preds_and_resids$residuals_cn))
    
        Shapiro-Wilk normality test
    
    data:  scale(meta_w_preds_and_resids$residuals_cn)
    W = 0.98414, p-value = 0.009555

Also, check the normality of random effects:

qqnorm(ranef(cn_fit)$`site:tree`$`(Intercept)`)
qqline(ranef(cn_fit)$`site:tree`$`(Intercept)`)
shapiro.test(ranef(cn_fit)$`site:tree`$`(Intercept)`)
    
        Shapiro-Wilk normality test
    
    data:  ranef(cn_fit)$`site:tree`$`(Intercept)`
    W = 0.90692, p-value = 0.1948

C:P

1. linearity of predictors

visual inspection, plot of residuals against response values:

2. homogeneity of variance

visual inspection, plot of residuals against fitted values:

statistical test, Levene’s test:

The assumption of equal variance of residuals is violated for several grouping factors, however this should not be a problem, since large sample sizes almost always result in a significant difference. The robustness of LMERs against violations like this has also been documented (Schielzeth et al 2020, DOI: 10.1111/2041-210X.13434).

leveneTest(residuals(cp_fit) ~ meta_w_preds_and_resids$tree)
    Levene's Test for Homogeneity of Variance (center = median)
           Df F value Pr(>F)
    group  11  1.1438 0.3283
          225
leveneTest(residuals(cp_fit) ~ meta_w_preds_and_resids$bb_colonization)
    Levene's Test for Homogeneity of Variance (center = median)
           Df F value      Pr(>F)    
    group   1  20.811 0.000008167 ***
          235                        
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
leveneTest(residuals(cp_fit) ~ meta_w_preds_and_resids$sampling_time)
    Levene's Test for Homogeneity of Variance (center = median)
           Df F value Pr(>F)
    group   1       0  0.998
          235
leveneTest(residuals(cp_fit) ~ meta_w_preds_and_resids$site)
    Levene's Test for Homogeneity of Variance (center = median)
           Df F value   Pr(>F)   
    group   3  4.1986 0.006443 **
          233                    
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3. normality of errors

visual inspection, histogram of residuals:

visual inspection, plot of residuals for each explanatory factor:

visual inspection, plot of measured value against residuals:

statistical test, Shapiro-Wilk test:

As with, homogeneity of variance, assumption violations occur here, but this should not be a problem for the same reasons as homogeneity of variance.

shapiro.test(meta_w_preds_and_resids$residuals_cp)
    
        Shapiro-Wilk normality test
    
    data:  meta_w_preds_and_resids$residuals_cp
    W = 0.96513, p-value = 0.00001518
shapiro.test(scale(meta_w_preds_and_resids$residuals_cp))
    
        Shapiro-Wilk normality test
    
    data:  scale(meta_w_preds_and_resids$residuals_cp)
    W = 0.96513, p-value = 0.00001518

Also, check the normality of random effects:

qqnorm(ranef(cp_fit)$`site:tree`$`(Intercept)`)
qqline(ranef(cp_fit)$`site:tree`$`(Intercept)`)
shapiro.test(ranef(cp_fit)$`site:tree`$`(Intercept)`)
    
        Shapiro-Wilk normality test
    
    data:  ranef(cp_fit)$`site:tree`$`(Intercept)`
    W = 0.97804, p-value = 0.9746

N:P

1. linearity of predictors

visual inspection, plot of residuals against response values:

2. homogeneity of variance

visual inspection, plot of residuals against fitted values:

statistical test, Levene’s test:

The assumption of equal variance of residuals is violated for several grouping factors, however this should not be a problem, since large sample sizes almost always result in a significant difference. The robustness of LMERs against violations like this has also been documented (Schielzeth et al 2020, DOI: 10.1111/2041-210X.13434).

leveneTest(residuals(np_fit) ~ meta_w_preds_and_resids$tree)
    Levene's Test for Homogeneity of Variance (center = median)
           Df F value Pr(>F)
    group  11  1.4346 0.1585
          225
leveneTest(residuals(np_fit) ~ meta_w_preds_and_resids$bb_colonization)
    Levene's Test for Homogeneity of Variance (center = median)
           Df F value  Pr(>F)  
    group   1  4.0556 0.04516 *
          235                  
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
leveneTest(residuals(np_fit) ~ meta_w_preds_and_resids$sampling_time)
    Levene's Test for Homogeneity of Variance (center = median)
           Df F value Pr(>F)
    group   1   0.132 0.7167
          235
leveneTest(residuals(np_fit) ~ meta_w_preds_and_resids$site)
    Levene's Test for Homogeneity of Variance (center = median)
           Df F value   Pr(>F)   
    group   3  4.6263 0.003654 **
          233                    
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3. normality of errors

visual inspection, histogram of residuals:

visual inspection, plot of residuals for each explanatory factor:

visual inspection, plot of measured value against residuals:

statistical test, Shapiro-Wilk test:

As with, homogeneity of variance, assumption violations occur here, but this should not be a problem for the same reasons as homogeneity of variance.

shapiro.test(meta_w_preds_and_resids$residuals_np)
    
        Shapiro-Wilk normality test
    
    data:  meta_w_preds_and_resids$residuals_np
    W = 0.98206, p-value = 0.004281
shapiro.test(scale(meta_w_preds_and_resids$residuals_np))
    
        Shapiro-Wilk normality test
    
    data:  scale(meta_w_preds_and_resids$residuals_np)
    W = 0.98206, p-value = 0.004281

Also, check the normality of random effects:

qqnorm(ranef(np_fit)$`site:tree`$`(Intercept)`)
qqline(ranef(np_fit)$`site:tree`$`(Intercept)`)
shapiro.test(ranef(np_fit)$`site:tree`$`(Intercept)`)
    
        Shapiro-Wilk normality test
    
    data:  ranef(np_fit)$`site:tree`$`(Intercept)`
    W = 0.95388, p-value = 0.6942




\(R^{2}\), ICC, Wald \(\chi^{2}\) test, 95% CI, and summary for each model

C:N

# R^2
performance::r2(cn_fit)
    # R2 for Mixed Models
    
      Conditional R2: 0.731
         Marginal R2: 0.469
# Interclass Correlation Coefficient (ICC)
performance::icc(cn_fit)
    # Intraclass Correlation Coefficient
    
         Adjusted ICC: 0.494
      Conditional ICC: 0.262
# Type II Wald Chi^2 test (Analysis of Deviance table produced)
Anova(cn_fit)
    Analysis of Deviance Table (Type II Wald chisquare tests)
    
    Response: cn_ratio
                                               Chisq Df            Pr(>Chisq)    
    fire_hist                                 0.1576  1              0.691353    
    bb_colonization                          51.5438  1    0.0000000000007002 ***
    sampling_time                           160.1106  1 < 0.00000000000000022 ***
    precol_cn_c                               5.2685  1              0.021714 *  
    fire_hist:bb_colonization                21.3155  1    0.0000038956150564 ***
    fire_hist:sampling_time                  44.9590  1    0.0000000000201202 ***
    bb_colonization:sampling_time            22.4943  1    0.0000021076720233 ***
    sampling_time:precol_cn_c                 9.4124  1              0.002155 ** 
    fire_hist:bb_colonization:sampling_time   0.7451  1              0.388047    
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 95% confident interval
confint(cn_fit)
    Computing profile confidence intervals ...
                                                     2.5 %       97.5 %
    .sig01                                       5.5040364  13.32601345
    .sigma                                       9.6762915  11.64216078
    (Intercept)                                103.3788808 114.05217704
    fire_hist1                                  -6.8616237   4.41613764
    bb_colonization1                            -7.7790426  -4.31815700
    sampling_time1                              -7.9413362  -4.47807442
    precol_cn_c                                  0.2093577   1.30118611
    fire_hist1:bb_colonization1                  2.3967124   5.85759801
    fire_hist1:sampling_time1                   -6.2650854  -2.67321686
    bb_colonization1:sampling_time1             -5.9932267  -2.52999415
    sampling_time1:precol_cn_c                  -0.3662111  -0.08199132
    fire_hist1:bb_colonization1:sampling_time1  -2.5007257   0.96250693
# model summary
summary(cn_fit)
    Cov prior  : site:tree ~ wishart(df = 3.5, scale = Inf, posterior.scale = cov, common.scale = TRUE)
    Prior dev  : 0.0377
    
    Linear mixed model fit by REML ['blmerMod']
    Formula: cn_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cn_c +      precol_cn_c:sampling_time + (1 | site:tree)
       Data: meta_chem_t1t2_cent
    
    REML criterion at convergence: 1809.7
    
    Scaled residuals: 
        Min      1Q  Median      3Q     Max 
    -3.3942 -0.5249  0.0247  0.5969  2.8124 
    
    Random effects:
     Groups    Name        Variance Std.Dev.
     site:tree (Intercept) 111.9    10.58   
     Residual              114.8    10.71   
    Number of obs: 237, groups:  site:tree, 12
    
    Fixed effects:
                                                Estimate Std. Error t value
    (Intercept)                                108.71630    3.18115  34.175
    fire_hist1                                  -1.22186    3.36210  -0.363
    bb_colonization1                            -6.05156    0.88995  -6.800
    sampling_time1                              -6.21011    0.89058  -6.973
    precol_cn_c                                  0.75530    0.32614   2.316
    fire_hist1:bb_colonization1                  4.12419    0.88995   4.634
    fire_hist1:sampling_time1                   -4.46995    0.92365  -4.839
    bb_colonization1:sampling_time1             -4.26122    0.89058  -4.785
    sampling_time1:precol_cn_c                  -0.22422    0.07308  -3.068
    fire_hist1:bb_colonization1:sampling_time1  -0.76871    0.89058  -0.863
    
    Correlation of Fixed Effects:
                (Intr) fr_hs1 bb_cl1 smpl_1 prcl__ fr_hst1:b_1 fr_hst1:s_1 b_1:_1 s_1:__
    fire_hist1   0.005                                                                  
    bb_colnztn1 -0.174 -0.015                                                           
    samplng_tm1 -0.006 -0.006  0.021                                                    
    precol_cn_c  0.004  0.324  0.000  0.000                                             
    fr_hst1:b_1 -0.016 -0.165  0.056  0.021  0.000                                      
    fr_hst1:s_1 -0.006 -0.006  0.020  0.057 -0.001  0.020                               
    bb_clnz1:_1  0.006  0.006 -0.021 -0.623  0.000 -0.021      -0.056                   
    smplng_1:__  0.000 -0.002 -0.002  0.004 -0.007 -0.002       0.265      -0.001       
    fr_h1:_1:_1  0.006  0.006 -0.021 -0.058  0.000 -0.021      -0.601       0.058 -0.001

C:P

# R^2
performance::r2(cp_fit)
    # R2 for Mixed Models
    
      Conditional R2: 0.774
         Marginal R2: 0.352
# Interclass Correlation Coefficient (ICC)
performance::icc(cp_fit)
    # Intraclass Correlation Coefficient
    
         Adjusted ICC: 0.651
      Conditional ICC: 0.422
# Type II Wald Chi^2 test (Analysis of Deviance table produced)
Anova(cp_fit)
    Analysis of Deviance Table (Type II Wald chisquare tests)
    
    Response: cp_ratio
                                              Chisq Df            Pr(>Chisq)    
    fire_hist                                0.3713  1             0.5422934    
    bb_colonization                          5.6866  1             0.0170945 *  
    sampling_time                           78.4676  1 < 0.00000000000000022 ***
    precol_cp_c                              2.5643  1             0.1093029    
    fire_hist:bb_colonization                0.0402  1             0.8410526    
    fire_hist:sampling_time                  3.2554  1             0.0711892 .  
    bb_colonization:sampling_time           33.0173  1        0.000000009134 ***
    sampling_time:precol_cp_c                0.3703  1             0.5428501    
    fire_hist:bb_colonization:sampling_time 14.7234  1             0.0001245 ***
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 95% confident interval
confint(cp_fit)
    Computing profile confidence intervals ...
                                                       2.5 %        97.5 %
    .sig01                                      197.14654597  461.47439522
    .sigma                                      242.66215053  291.96579135
    (Intercept)                                1704.17086516 2067.98868988
    fire_hist1                                 -303.35346749  145.34857058
    bb_colonization1                             13.59968590  100.43771280
    sampling_time1                             -113.89826848  -27.03452531
    precol_cp_c                                  -0.01600956    0.68101675
    fire_hist1:bb_colonization1                 -46.18613658   40.65189032
    fire_hist1:sampling_time1                   -35.44100908   64.82133279
    bb_colonization1:sampling_time1            -176.63560881  -89.77510786
    sampling_time1:precol_cp_c                   -0.08646784    0.04552046
    fire_hist1:bb_colonization1:sampling_time1 -129.24737792  -42.38687697
# model summary
summary(cp_fit)
    Cov prior  : site:tree ~ wishart(df = 3.5, scale = Inf, posterior.scale = cov, common.scale = TRUE)
    Prior dev  : -0.9366
    
    Linear mixed model fit by REML ['blmerMod']
    Formula: cp_ratio ~ fire_hist * bb_colonization * sampling_time + precol_cp_c +      precol_cp_c:sampling_time + (1 | site:tree)
       Data: meta_chem_t1t2_cent
    
    REML criterion at convergence: 3294.1
    
    Scaled residuals: 
        Min      1Q  Median      3Q     Max 
    -2.7138 -0.6336 -0.0479  0.5273  4.2803 
    
    Random effects:
     Groups    Name        Variance Std.Dev.
     site:tree (Intercept) 134799   367.1   
     Residual               72198   268.7   
    Number of obs: 237, groups:  site:tree, 12
    
    Fixed effects:
                                                 Estimate Std. Error t value
    (Intercept)                                1886.05435  108.31653  17.412
    fire_hist1                                  -79.03938  133.64113  -0.591
    bb_colonization1                             57.12225   22.33091   2.558
    sampling_time1                              -70.57560   22.33759  -3.159
    precol_cp_c                                   0.33247    0.20776   1.600
    fire_hist1:bb_colonization1                  -2.66357   22.33091  -0.119
    fire_hist1:sampling_time1                    14.51387   25.78271   0.563
    bb_colonization1:sampling_time1            -133.09680   22.33676  -5.959
    sampling_time1:precol_cp_c                   -0.02065    0.03394  -0.609
    fire_hist1:bb_colonization1:sampling_time1  -85.70857   22.33676  -3.837
    
    Correlation of Fixed Effects:
                (Intr) fr_hs1 bb_cl1 smpl_1 prcl__ fr_hst1:b_1 fr_hst1:s_1 b_1:_1 s_1:__
    fire_hist1   0.006                                                                  
    bb_colnztn1 -0.128 -0.010                                                           
    samplng_tm1 -0.004  0.000  0.021                                                    
    precol_cp_c  0.007  0.586 -0.001  0.005                                             
    fr_hst1:b_1 -0.012 -0.104  0.057  0.021 -0.001                                      
    fr_hst1:s_1 -0.001  0.003  0.003  0.055  0.005  0.003                               
    bb_clnz1:_1  0.004  0.001 -0.021 -0.623 -0.005 -0.021      -0.052                   
    smplng_1:__  0.006  0.006 -0.031  0.010  0.002 -0.031       0.499      -0.004       
    fr_h1:_1:_1  0.004  0.001 -0.021 -0.058 -0.005 -0.021      -0.542       0.058 -0.004

N:P

# R^2
performance::r2(np_fit)
    # R2 for Mixed Models
    
      Conditional R2: 0.777
         Marginal R2: 0.320
# Interclass Correlation Coefficient (ICC)
performance::icc(np_fit)
    # Intraclass Correlation Coefficient
    
         Adjusted ICC: 0.672
      Conditional ICC: 0.457
# Type II Wald Chi^2 test (Analysis of Deviance table produced)
Anova(np_fit)
    Analysis of Deviance Table (Type II Wald chisquare tests)
    
    Response: np_ratio
                                              Chisq Df            Pr(>Chisq)    
    fire_hist                                1.6250  1             0.2023915    
    bb_colonization                         72.6535  1 < 0.00000000000000022 ***
    sampling_time                            0.2238  1             0.6361339    
    precol_np_c                              0.8075  1             0.3688732    
    fire_hist:bb_colonization               23.4658  1           0.000001272 ***
    fire_hist:sampling_time                 14.2598  1             0.0001592 ***
    bb_colonization:sampling_time            7.5743  1             0.0059206 ** 
    sampling_time:precol_np_c               15.6855  1           0.000074797 ***
    fire_hist:bb_colonization:sampling_time 18.7486  1           0.000014913 ***
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 95% confident interval
confint(np_fit)
    Computing profile confidence intervals ...
                                                     2.5 %     97.5 %
    .sig01                                      1.62591532  3.7927453
    .sigma                                      1.90514682  2.2922253
    (Intercept)                                16.03280506 19.0184422
    fire_hist1                                 -2.52923167  0.9616255
    bb_colonization1                            1.12649924  1.8082652
    sampling_time1                              0.07541901  0.7573734
    precol_np_c                                -0.12398286  0.4149505
    fire_hist1:bb_colonization1                -1.17498193 -0.4932160
    fire_hist1:sampling_time1                   0.70462091  1.4620188
    bb_colonization1:sampling_time1            -0.86713001 -0.1851935
    sampling_time1:precol_np_c                  0.05063376  0.1480716
    fire_hist1:bb_colonization1:sampling_time1 -1.10066365 -0.4187272
# model summary
summary(np_fit)
    Cov prior  : site:tree ~ wishart(df = 3.5, scale = Inf, posterior.scale = cov, common.scale = TRUE)
    Prior dev  : -1.0744
    
    Linear mixed model fit by REML ['blmerMod']
    Formula: np_ratio ~ fire_hist * bb_colonization * sampling_time + precol_np_c +      precol_np_c:sampling_time + (1 | site:tree)
       Data: meta_chem_t1t2_cent
    
    REML criterion at convergence: 1076.1
    
    Scaled residuals: 
        Min      1Q  Median      3Q     Max 
    -3.9613 -0.6221 -0.0219  0.5938  3.2688 
    
    Random effects:
     Groups    Name        Variance Std.Dev.
     site:tree (Intercept) 9.109    3.018   
     Residual              4.450    2.110   
    Number of obs: 237, groups:  site:tree, 12
    
    Fixed effects:
                                               Estimate Std. Error t value
    (Intercept)                                17.52536    0.88872  19.720
    fire_hist1                                 -0.78413    1.03940  -0.754
    bb_colonization1                            1.46847    0.17532   8.376
    sampling_time1                              0.41604    0.17537   2.372
    precol_np_c                                 0.14547    0.16059   0.906
    fire_hist1:bb_colonization1                -0.83301    0.17532  -4.751
    fire_hist1:sampling_time1                   1.08258    0.19477   5.558
    bb_colonization1:sampling_time1            -0.52581    0.17537  -2.998
    sampling_time1:precol_np_c                  0.09924    0.02506   3.960
    fire_hist1:bb_colonization1:sampling_time1 -0.75934    0.17537  -4.330
    
    Correlation of Fixed Effects:
                (Intr) fr_hs1 bb_cl1 smpl_1 prcl__ fr_hst1:b_1 fr_hst1:s_1 b_1:_1 s_1:__
    fire_hist1   0.005                                                                  
    bb_colnztn1 -0.123 -0.010                                                           
    samplng_tm1 -0.004 -0.001  0.021                                                    
    precol_np_c  0.006  0.519 -0.001  0.005                                             
    fr_hst1:b_1 -0.011 -0.105  0.057  0.021 -0.001                                      
    fr_hst1:s_1 -0.001  0.002  0.006  0.056  0.005  0.006                               
    bb_clnz1:_1  0.004  0.001 -0.021 -0.623 -0.005 -0.021      -0.054                   
    smplng_1:__  0.006  0.006 -0.031  0.008  0.002 -0.031       0.435      -0.004       
    fr_h1:_1:_1  0.004  0.001 -0.021 -0.058 -0.005 -0.021      -0.563       0.058 -0.004



Linear mixed-effects models, partially Bayesian implementation, burnt and unburnt subsets

Fire history does not have a significant main effect in any of the three models, but there are significant interactions that include fire history, in different combinations for each model. To begin teasing apart the relationships between fire history and the remaining two explanatory variables, bark beetle colonization status and sampling time, we can analyze the data for samples from burnt and unburnt sites separately. Fitted models follow the same formula as the models above, only with fire_hist terms excluded.

fit models

# subset by fire history, center precol_XX variables, and apply sum-to-zero contrasts

meta_chem_t1t2_cent_burnt <- subset(meta_chem_t1t2, fire_hist == "burnt") %>% droplevels()
contrasts(meta_chem_t1t2_cent_burnt$sampling_time) <- contr.sum
contrasts(meta_chem_t1t2_cent_burnt$bb_colonization) <- contr.sum
meta_chem_t1t2_cent_burnt$precol_cn_c <- scale(meta_chem_t1t2_cent_burnt$precol_cn, scale = FALSE)
meta_chem_t1t2_cent_burnt$precol_cp_c <- scale(meta_chem_t1t2_cent_burnt$precol_cp, scale = FALSE)
meta_chem_t1t2_cent_burnt$precol_np_c <- scale(meta_chem_t1t2_cent_burnt$precol_np, scale = FALSE)

meta_chem_t1t2_cent_unburnt <- subset(meta_chem_t1t2, fire_hist == "unburnt") %>% droplevels()
contrasts(meta_chem_t1t2_cent_unburnt$sampling_time) <- contr.sum
contrasts(meta_chem_t1t2_cent_unburnt$bb_colonization) <- contr.sum
meta_chem_t1t2_cent_unburnt$precol_cn_c <- scale(meta_chem_t1t2_cent_unburnt$precol_cn, scale = FALSE)
meta_chem_t1t2_cent_unburnt$precol_cp_c <- scale(meta_chem_t1t2_cent_unburnt$precol_cp, scale = FALSE)
meta_chem_t1t2_cent_unburnt$precol_np_c <- scale(meta_chem_t1t2_cent_unburnt$precol_np, scale = FALSE)
# fit models following structure from full model (model including samples from both burnt and unburnt sites)

cn_fit_burnt <- blme::blmer(cn_ratio ~ bb_colonization * sampling_time + precol_cn_c +
                                precol_cn_c:sampling_time + (1|site:tree), data = meta_chem_t1t2_cent_burnt)

cp_fit_burnt <- blme::blmer(cp_ratio ~ bb_colonization * sampling_time + precol_cp_c +
                                precol_cp_c:sampling_time + (1|site:tree), data = meta_chem_t1t2_cent_burnt)

np_fit_burnt <- blme::blmer(np_ratio ~ bb_colonization * sampling_time + precol_np_c +
                                precol_np_c:sampling_time + (1|site:tree), data = meta_chem_t1t2_cent_burnt)

cn_fit_unburnt <- blme::blmer(cn_ratio ~ bb_colonization * sampling_time + precol_cn_c +
                                  precol_cn_c:sampling_time + (1|site:tree), data = meta_chem_t1t2_cent_unburnt)

cp_fit_unburnt <- blme::blmer(cp_ratio ~ bb_colonization * sampling_time + precol_cp_c +
                                  precol_cp_c:sampling_time + (1|site:tree), data = meta_chem_t1t2_cent_unburnt)

np_fit_unburnt <- blme::blmer(np_ratio ~ bb_colonization * sampling_time + precol_np_c +
                                  precol_np_c:sampling_time + (1|site:tree), data = meta_chem_t1t2_cent_unburnt)

Model results

C:N burnt
# R^2
performance::r2(cn_fit_burnt)
    # R2 for Mixed Models
    
      Conditional R2: 0.762
         Marginal R2: 0.209
# Interclass Correlation Coefficient (ICC)
performance::icc(cn_fit_burnt)
    # Intraclass Correlation Coefficient
    
         Adjusted ICC: 0.699
      Conditional ICC: 0.553
# Type II Wald Chi^2 test (Analysis of Deviance table produced)
Anova(cn_fit_burnt)
    Analysis of Deviance Table (Type II Wald chisquare tests)
    
    Response: cn_ratio
                                    Chisq Df          Pr(>Chisq)    
    bb_colonization               59.7424  1 0.00000000000001081 ***
    sampling_time                 18.8727  1 0.00001397384680382 ***
    precol_cn_c                    0.3335  1            0.563587    
    bb_colonization:sampling_time  7.0376  1            0.007982 ** 
    sampling_time:precol_cn_c      2.7833  1            0.095251 .  
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 95% confident interval
confint(cn_fit_burnt)
    Computing profile confidence intervals ...
                                           2.5 %       97.5 %
    .sig01                            6.59502112  22.69988827
    .sigma                           10.07671203  13.06995771
    (Intercept)                     101.56126974 123.29901915
    bb_colonization1                -12.74964396  -7.60186605
    sampling_time1                   -5.05382053   0.09395739
    precol_cn_c                      -0.88052444   1.99789667
    bb_colonization1:sampling_time1  -6.06639004  -0.91861212
    sampling_time1:precol_cn_c       -0.04028855   0.50908560
# model summary
summary(cn_fit_burnt)
    Cov prior  : site:tree ~ wishart(df = 3.5, scale = Inf, posterior.scale = cov, common.scale = TRUE)
    Prior dev  : -1.2626
    
    Linear mixed model fit by REML ['blmerMod']
    Formula: cn_ratio ~ bb_colonization * sampling_time + precol_cn_c + precol_cn_c:sampling_time +      (1 | site:tree)
       Data: meta_chem_t1t2_cent_burnt
    
    REML criterion at convergence: 932.4
    
    Scaled residuals: 
         Min       1Q   Median       3Q      Max 
    -2.64507 -0.55870 -0.08217  0.63479  2.49597 
    
    Random effects:
     Groups    Name        Variance Std.Dev.
     site:tree (Intercept) 308.9    17.57   
     Residual              133.1    11.54   
    Number of obs: 120, groups:  site:tree, 6
    
    Fixed effects:
                                    Estimate Std. Error t value
    (Intercept)                     112.4301     7.2946  15.413
    bb_colonization1                -10.1758     1.3165  -7.729
    sampling_time1                   -2.4799     1.3165  -1.884
    precol_cn_c                       0.5587     0.9674   0.578
    bb_colonization1:sampling_time1  -3.4925     1.3165  -2.653
    sampling_time1:precol_cn_c        0.2344     0.1405   1.668
    
    Correlation of Fixed Effects:
                (Intr) bb_cl1 smpl_1 prcl__ b_1:_1
    bb_colnztn1 -0.108                            
    samplng_tm1  0.000  0.000                     
    precol_cn_c  0.000  0.000  0.000              
    bb_clnz1:_1  0.000  0.000 -0.600  0.000       
    smplng_1:__  0.000  0.000  0.000  0.000  0.000
C:N unburnt
# R^2
performance::r2(cn_fit_unburnt)
    # R2 for Mixed Models
    
      Conditional R2: 0.814
         Marginal R2: 0.742
# Interclass Correlation Coefficient (ICC)
performance::icc(cn_fit_unburnt)
    # Intraclass Correlation Coefficient
    
         Adjusted ICC: 0.277
      Conditional ICC: 0.072
# Type II Wald Chi^2 test (Analysis of Deviance table produced)
Anova(cn_fit_unburnt)
    Analysis of Deviance Table (Type II Wald chisquare tests)
    
    Response: cn_ratio
                                    Chisq Df            Pr(>Chisq)    
    bb_colonization                 4.003  1               0.04542 *  
    sampling_time                 265.206  1 < 0.00000000000000022 ***
    precol_cn_c                    16.026  1         0.00006246412 ***
    bb_colonization:sampling_time  22.664  1         0.00000192930 ***
    sampling_time:precol_cn_c      36.421  1         0.00000000159 ***
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 95% confident interval
confint(cn_fit_unburnt)
    Computing profile confidence intervals ...
                                          2.5 %      97.5 %
    .sig01                            0.9878653   6.8873900
    .sigma                            7.6181318   9.9154532
    (Intercept)                     101.3352645 108.4777333
    bb_colonization1                 -3.9318497   0.1943767
    sampling_time1                  -11.9814505  -7.8535201
    precol_cn_c                       0.5513502   1.1491094
    bb_colonization1:sampling_time1  -7.0953268  -2.9673907
    sampling_time1:precol_cn_c       -0.5742674  -0.2929042
# model summary
summary(cn_fit_unburnt)
    Cov prior  : site:tree ~ wishart(df = 3.5, scale = Inf, posterior.scale = cov, common.scale = TRUE)
    Prior dev  : 1.4356
    
    Linear mixed model fit by REML ['blmerMod']
    Formula: cn_ratio ~ bb_colonization * sampling_time + precol_cn_c + precol_cn_c:sampling_time +      (1 | site:tree)
       Data: meta_chem_t1t2_cent_unburnt
    
    REML criterion at convergence: 841.8
    
    Scaled residuals: 
        Min      1Q  Median      3Q     Max 
    -4.1720 -0.5802 -0.0368  0.5781  2.1912 
    
    Random effects:
     Groups    Name        Variance Std.Dev.
     site:tree (Intercept) 29.29    5.412   
     Residual              76.26    8.733   
    Number of obs: 117, groups:  site:tree, 6
    
    Fixed effects:
                                     Estimate Std. Error t value
    (Intercept)                     104.91608    2.44795  42.859
    bb_colonization1                 -1.90597    1.05431  -1.808
    sampling_time1                   -9.92338    1.05561  -9.401
    precol_cn_c                       0.85035    0.20777   4.093
    bb_colonization1:sampling_time1  -5.02546    1.05562  -4.761
    sampling_time1:precol_cn_c       -0.43399    0.07191  -6.035
    
    Correlation of Fixed Effects:
                (Intr) bb_cl1 smpl_1 prcl__ b_1:_1
    bb_colnztn1 -0.277                            
    samplng_tm1 -0.017  0.040                     
    precol_cn_c -0.001 -0.001  0.001              
    bb_clnz1:_1  0.017 -0.040 -0.644 -0.001       
    smplng_1:__  0.001 -0.003  0.001 -0.015 -0.002






C:P burnt
# R^2
performance::r2(cp_fit_burnt)
    # R2 for Mixed Models
    
      Conditional R2: 0.724
         Marginal R2: 0.310
# Interclass Correlation Coefficient (ICC)
performance::icc(cp_fit_burnt)
    # Intraclass Correlation Coefficient
    
         Adjusted ICC: 0.601
      Conditional ICC: 0.415
# Type II Wald Chi^2 test (Analysis of Deviance table produced)
Anova(cp_fit_burnt)
    Analysis of Deviance Table (Type II Wald chisquare tests)
    
    Response: cp_ratio
                                    Chisq Df Pr(>Chisq)    
    bb_colonization                3.0299  1    0.08174 .  
    sampling_time                 19.4656  1 0.00001024 ***
    precol_cp_c                    3.5148  1    0.06082 .  
    bb_colonization:sampling_time  1.9036  1    0.16768    
    sampling_time:precol_cp_c      0.1096  1    0.74056    
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 95% confident interval
confint(cp_fit_burnt)
    Computing profile confidence intervals ...
                                            2.5 %       97.5 %
    .sig01                           133.09046863  476.3141237
    .sigma                           262.91147842  341.0082471
    (Intercept)                     1859.40513592 2318.8910835
    bb_colonization1                  -7.36951045  126.9411554
    sampling_time1                  -159.95116186  -25.6404960
    precol_cp_c                        0.11671905    0.9987833
    bb_colonization1:sampling_time1 -114.54356382   19.7671020
    sampling_time1:precol_cp_c        -0.08664678    0.1219758
# model summary
summary(cp_fit_burnt)
    Cov prior  : site:tree ~ wishart(df = 3.5, scale = Inf, posterior.scale = cov, common.scale = TRUE)
    Prior dev  : -0.613
    
    Linear mixed model fit by REML ['blmerMod']
    Formula: cp_ratio ~ bb_colonization * sampling_time + precol_cp_c + precol_cp_c:sampling_time +      (1 | site:tree)
       Data: meta_chem_t1t2_cent_burnt
    
    REML criterion at convergence: 1691.3
    
    Scaled residuals: 
        Min      1Q  Median      3Q     Max 
    -2.3517 -0.6346 -0.0307  0.5378  3.7755 
    
    Random effects:
     Groups    Name        Variance Std.Dev.
     site:tree (Intercept) 136338   369.2   
     Residual               90599   301.0   
    Number of obs: 120, groups:  site:tree, 6
    
    Fixed effects:
                                      Estimate Std. Error t value
    (Intercept)                     2089.14801  154.60499  13.513
    bb_colonization1                  59.78582   34.34644   1.741
    sampling_time1                   -92.79583   34.34644  -2.702
    precol_cp_c                        0.55775    0.29750   1.875
    bb_colonization1:sampling_time1  -47.38823   34.34644  -1.380
    sampling_time1:precol_cp_c         0.01766    0.05335   0.331
    
    Correlation of Fixed Effects:
                (Intr) bb_cl1 smpl_1 prcl__ b_1:_1
    bb_colnztn1 -0.133                            
    samplng_tm1  0.000  0.000                     
    precol_cp_c  0.000  0.000  0.000              
    bb_clnz1:_1  0.000  0.000 -0.600  0.000       
    smplng_1:__  0.000  0.000  0.000  0.000  0.000
C:P unburnt
# R^2
performance::r2(cp_fit_unburnt)
    # R2 for Mixed Models
    
      Conditional R2: 0.843
         Marginal R2: 0.209
# Interclass Correlation Coefficient (ICC)
performance::icc(cp_fit_unburnt)
    # Intraclass Correlation Coefficient
    
         Adjusted ICC: 0.802
      Conditional ICC: 0.634
# Type II Wald Chi^2 test (Analysis of Deviance table produced)
Anova(cp_fit_unburnt)
    Analysis of Deviance Table (Type II Wald chisquare tests)
    
    Response: cp_ratio
                                    Chisq Df            Pr(>Chisq)    
    bb_colonization                2.9804  1               0.08428 .  
    sampling_time                 80.2984  1 < 0.00000000000000022 ***
    precol_cp_c                    0.0930  1               0.76045    
    bb_colonization:sampling_time 62.6113  1  0.000000000000002518 ***
    sampling_time:precol_cp_c      2.1756  1               0.14021    
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 95% confident interval
confint(cp_fit_unburnt)
    Computing profile confidence intervals ...
                                           2.5 %        97.5 %
    .sig01                           177.2625408  593.65591054
    .sigma                           198.8528551  258.82711679
    (Intercept)                     1395.5402688 1961.49638884
    bb_colonization1                   2.1180279  110.01652103
    sampling_time1                  -102.4289405    5.41979579
    precol_cp_c                       -0.4326445    0.65172278
    bb_colonization1:sampling_time1 -272.2899292 -164.43880385
    sampling_time1:precol_cp_c        -0.1405676    0.02004044
# model summary
summary(cp_fit_unburnt)
    Cov prior  : site:tree ~ wishart(df = 3.5, scale = Inf, posterior.scale = cov, common.scale = TRUE)
    Prior dev  : -2.0988
    
    Linear mixed model fit by REML ['blmerMod']
    Formula: cp_ratio ~ bb_colonization * sampling_time + precol_cp_c + precol_cp_c:sampling_time +      (1 | site:tree)
       Data: meta_chem_t1t2_cent_unburnt
    
    REML criterion at convergence: 1590.5
    
    Scaled residuals: 
         Min       1Q   Median       3Q      Max 
    -2.83998 -0.65983 -0.09913  0.62621  2.39668 
    
    Random effects:
     Groups    Name        Variance Std.Dev.
     site:tree (Intercept) 210862   459.2   
     Residual               52040   228.1   
    Number of obs: 117, groups:  site:tree, 6
    
    Fixed effects:
                                      Estimate Std. Error t value
    (Intercept)                     1678.47939  189.48902   8.858
    bb_colonization1                  56.28936   27.59790   2.040
    sampling_time1                   -48.58482   27.58574  -1.761
    precol_cp_c                        0.10952    0.36345   0.301
    bb_colonization1:sampling_time1 -218.28334   27.58634  -7.913
    sampling_time1:precol_cp_c        -0.06059    0.04108  -1.475
    
    Correlation of Fixed Effects:
                (Intr) bb_cl1 smpl_1 prcl__ b_1:_1
    bb_colnztn1 -0.094                            
    samplng_tm1 -0.006  0.040                     
    precol_cp_c -0.004 -0.001  0.007              
    bb_clnz1:_1  0.006 -0.040 -0.644 -0.007       
    smplng_1:__  0.009 -0.061  0.004  0.002 -0.008






N:P burnt
# R^2
performance::r2(np_fit_burnt)
    # R2 for Mixed Models
    
      Conditional R2: 0.733
         Marginal R2: 0.348
# Interclass Correlation Coefficient (ICC)
performance::icc(np_fit_burnt)
    # Intraclass Correlation Coefficient
    
         Adjusted ICC: 0.591
      Conditional ICC: 0.385
# Type II Wald Chi^2 test (Analysis of Deviance table produced)
Anova(np_fit_burnt)
    Analysis of Deviance Table (Type II Wald chisquare tests)
    
    Response: np_ratio
                                    Chisq Df            Pr(>Chisq)    
    bb_colonization               78.4014  1 < 0.00000000000000022 ***
    sampling_time                  0.8924  1                0.3448    
    precol_np_c                    1.7991  1                0.1798    
    bb_colonization:sampling_time  0.8073  1                0.3689    
    sampling_time:precol_np_c     21.2492  1           0.000004033 ***
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 95% confident interval
confint(np_fit_burnt)
    Computing profile confidence intervals ...
                                          2.5 %     97.5 %
    .sig01                           0.98100045  3.5286402
    .sigma                           1.98965206  2.5806700
    (Intercept)                     17.08966961 20.4967022
    bb_colonization1                 1.79326545  2.8096969
    sampling_time1                  -0.84477429  0.1716571
    precol_np_c                     -0.02755255  0.5535212
    bb_colonization1:sampling_time1 -0.27468208  0.7417494
    sampling_time1:precol_np_c       0.09525869  0.2355934
# model summary
summary(np_fit_burnt)
    Cov prior  : site:tree ~ wishart(df = 3.5, scale = Inf, posterior.scale = cov, common.scale = TRUE)
    Prior dev  : -0.5496
    
    Linear mixed model fit by REML ['blmerMod']
    Formula: np_ratio ~ bb_colonization * sampling_time + precol_np_c + precol_np_c:sampling_time +      (1 | site:tree)
       Data: meta_chem_t1t2_cent_burnt
    
    REML criterion at convergence: 559.7
    
    Scaled residuals: 
        Min      1Q  Median      3Q     Max 
    -3.0286 -0.6498 -0.1108  0.6674  2.9046 
    
    Random effects:
     Groups    Name        Variance Std.Dev.
     site:tree (Intercept) 7.485    2.736   
     Residual              5.189    2.278   
    Number of obs: 120, groups:  site:tree, 6
    
    Fixed effects:
                                    Estimate Std. Error t value
    (Intercept)                     18.79319    1.14673  16.389
    bb_colonization1                 2.30148    0.25992   8.854
    sampling_time1                  -0.33656    0.25992  -1.295
    precol_np_c                      0.26298    0.19607   1.341
    bb_colonization1:sampling_time1  0.23353    0.25992   0.898
    sampling_time1:precol_np_c       0.16543    0.03589   4.610
    
    Correlation of Fixed Effects:
                (Intr) bb_cl1 smpl_1 prcl__ b_1:_1
    bb_colnztn1 -0.136                            
    samplng_tm1  0.000  0.000                     
    precol_np_c  0.000  0.000  0.000              
    bb_clnz1:_1  0.000  0.000 -0.600  0.000       
    smplng_1:__  0.000  0.000  0.000  0.000  0.000
N:P unburnt
# R^2
performance::r2(np_fit_unburnt)
    # R2 for Mixed Models
    
      Conditional R2: 0.854
         Marginal R2: 0.059
# Interclass Correlation Coefficient (ICC)
performance::icc(np_fit_unburnt)
    # Intraclass Correlation Coefficient
    
         Adjusted ICC: 0.844
      Conditional ICC: 0.795
# Type II Wald Chi^2 test (Analysis of Deviance table produced)
Anova(np_fit_unburnt)
    Analysis of Deviance Table (Type II Wald chisquare tests)
    
    Response: np_ratio
                                    Chisq Df     Pr(>Chisq)    
    bb_colonization                8.1073  1       0.004409 ** 
    sampling_time                  3.9313  1       0.047395 *  
    precol_np_c                    0.0002  1       0.988471    
    bb_colonization:sampling_time 34.0067  1 0.000000005492 ***
    sampling_time:precol_np_c      0.1284  1       0.720080    
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 95% confident interval
confint(np_fit_unburnt)
    Computing profile confidence intervals ...
                                          2.5 %      97.5 %
    .sig01                           1.64842715  5.46636294
    .sigma                           1.58115782  2.05802842
    (Intercept)                     13.60947431 18.80887503
    bb_colonization1                 0.24428682  1.10241383
    sampling_time1                   0.72730751  1.58486209
    precol_np_c                     -0.50505491  0.49601447
    bb_colonization1:sampling_time1 -1.70776804 -0.85019469
    sampling_time1:precol_np_c      -0.05216145  0.07605762
# model summary
summary(np_fit_unburnt)
    Cov prior  : site:tree ~ wishart(df = 3.5, scale = Inf, posterior.scale = cov, common.scale = TRUE)
    Prior dev  : -2.5378
    
    Linear mixed model fit by REML ['blmerMod']
    Formula: np_ratio ~ bb_colonization * sampling_time + precol_np_c + precol_np_c:sampling_time +      (1 | site:tree)
       Data: meta_chem_t1t2_cent_unburnt
    
    REML criterion at convergence: 499.9
    
    Scaled residuals: 
        Min      1Q  Median      3Q     Max 
    -4.0005 -0.5508  0.0111  0.5687  3.3721 
    
    Random effects:
     Groups    Name        Variance Std.Dev.
     site:tree (Intercept) 17.87    4.227   
     Residual               3.29    1.814   
    Number of obs: 117, groups:  site:tree, 6
    
    Fixed effects:
                                     Estimate Std. Error t value
    (Intercept)                     16.208832   1.739468   9.318
    bb_colonization1                 0.675333   0.219492   3.077
    sampling_time1                   1.156293   0.219352   5.271
    precol_np_c                     -0.004533   0.335174  -0.014
    bb_colonization1:sampling_time1 -1.279185   0.219357  -5.832
    sampling_time1:precol_np_c       0.011753   0.032796   0.358
    
    Correlation of Fixed Effects:
                (Intr) bb_cl1 smpl_1 prcl__ b_1:_1
    bb_colnztn1 -0.081                            
    samplng_tm1 -0.005  0.040                     
    precol_np_c -0.004 -0.001  0.006              
    bb_clnz1:_1  0.005 -0.040 -0.644 -0.006       
    smplng_1:__  0.008 -0.064  0.004  0.003 -0.008



Session Info

devtools::session_info()
    ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
     setting  value                       
     version  R version 3.6.2 (2019-12-12)
     os       macOS High Sierra 10.13.6   
     system   x86_64, darwin15.6.0        
     ui       X11                         
     language (EN)                        
     collate  en_US.UTF-8                 
     ctype    en_US.UTF-8                 
     tz       Europe/Stockholm            
     date     2021-12-02                  
    
    ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
     package     * version date       lib source        
     abind         1.4-5   2016-07-21 [1] CRAN (R 3.6.0)
     assertthat    0.2.1   2019-03-21 [1] CRAN (R 3.6.0)
     backports     1.2.1   2020-12-09 [1] CRAN (R 3.6.2)
     blme        * 1.0-5   2021-01-05 [1] CRAN (R 3.6.2)
     boot          1.3-28  2021-05-03 [1] CRAN (R 3.6.2)
     broom         0.7.7   2021-06-13 [1] CRAN (R 3.6.2)
     bslib         0.2.5.1 2021-05-18 [1] CRAN (R 3.6.2)
     cachem        1.0.5   2021-05-15 [1] CRAN (R 3.6.2)
     callr         3.7.0   2021-04-20 [1] CRAN (R 3.6.2)
     car         * 3.0-10  2020-09-29 [1] CRAN (R 3.6.2)
     carData     * 3.0-4   2020-05-22 [1] CRAN (R 3.6.2)
     cellranger    1.1.0   2016-07-27 [1] CRAN (R 3.6.0)
     cli           2.5.0   2021-04-26 [1] CRAN (R 3.6.2)
     cluster       2.1.2   2021-04-17 [1] CRAN (R 3.6.2)
     colorspace    2.0-1   2021-05-04 [1] CRAN (R 3.6.2)
     cowplot       1.1.1   2020-12-30 [1] CRAN (R 3.6.2)
     crayon        1.4.1   2021-02-08 [1] CRAN (R 3.6.2)
     curl          4.3.1   2021-04-30 [1] CRAN (R 3.6.2)
     data.table    1.14.0  2021-02-21 [1] CRAN (R 3.6.2)
     DBI           1.1.1   2021-01-15 [1] CRAN (R 3.6.2)
     desc          1.3.0   2021-03-05 [1] CRAN (R 3.6.2)
     devtools      2.4.2   2021-06-07 [1] CRAN (R 3.6.2)
     digest        0.6.27  2020-10-24 [1] CRAN (R 3.6.2)
     dplyr       * 1.0.6   2021-05-05 [1] CRAN (R 3.6.2)
     ellipsis      0.3.2   2021-04-29 [1] CRAN (R 3.6.2)
     evaluate      0.14    2019-05-28 [1] CRAN (R 3.6.0)
     fansi         0.5.0   2021-05-25 [1] CRAN (R 3.6.2)
     farver        2.1.0   2021-02-28 [1] CRAN (R 3.6.2)
     fastmap       1.1.0   2021-01-25 [1] CRAN (R 3.6.2)
     forcats       0.5.1   2021-01-27 [1] CRAN (R 3.6.2)
     foreign       0.8-76  2020-03-03 [1] CRAN (R 3.6.0)
     fs            1.5.0   2020-07-31 [1] CRAN (R 3.6.2)
     generics      0.1.0   2020-10-31 [1] CRAN (R 3.6.2)
     ggplot2     * 3.3.3   2020-12-30 [1] CRAN (R 3.6.2)
     ggpubr        0.4.0   2020-06-27 [1] CRAN (R 3.6.2)
     ggsignif      0.6.2   2021-06-14 [1] CRAN (R 3.6.2)
     glue          1.4.2   2020-08-27 [1] CRAN (R 3.6.2)
     gridExtra     2.3     2017-09-09 [1] CRAN (R 3.6.0)
     gtable        0.3.0   2019-03-25 [1] CRAN (R 3.6.0)
     haven         2.4.1   2021-04-23 [1] CRAN (R 3.6.2)
     highr         0.9     2021-04-16 [1] CRAN (R 3.6.2)
     hms           1.1.0   2021-05-17 [1] CRAN (R 3.6.2)
     htmltools     0.5.1.1 2021-01-22 [1] CRAN (R 3.6.2)
     insight       0.14.1  2021-05-28 [1] CRAN (R 3.6.2)
     jquerylib     0.1.4   2021-04-26 [1] CRAN (R 3.6.2)
     jsonlite      1.7.2   2020-12-09 [1] CRAN (R 3.6.2)
     knitr         1.33    2021-04-24 [1] CRAN (R 3.6.2)
     labeling      0.4.2   2020-10-20 [1] CRAN (R 3.6.2)
     lattice     * 0.20-44 2021-05-02 [1] CRAN (R 3.6.2)
     lifecycle     1.0.0   2021-02-15 [1] CRAN (R 3.6.2)
     lme4        * 1.1-27  2021-05-15 [1] CRAN (R 3.6.2)
     magrittr      2.0.1   2020-11-17 [1] CRAN (R 3.6.2)
     maps          3.3.0   2018-04-03 [1] CRAN (R 3.6.0)
     MASS          7.3-54  2021-05-03 [1] CRAN (R 3.6.2)
     Matrix      * 1.3-4   2021-06-01 [1] CRAN (R 3.6.2)
     memoise       2.0.0   2021-01-26 [1] CRAN (R 3.6.2)
     mgcv          1.8-36  2021-06-01 [1] CRAN (R 3.6.2)
     minqa         1.2.4   2014-10-09 [1] CRAN (R 3.6.0)
     modelr      * 0.1.8   2020-05-19 [1] CRAN (R 3.6.2)
     munsell       0.5.0   2018-06-12 [1] CRAN (R 3.6.0)
     nlme          3.1-152 2021-02-04 [1] CRAN (R 3.6.2)
     nloptr        1.2.2.2 2020-07-02 [1] CRAN (R 3.6.2)
     openxlsx      4.2.3   2020-10-27 [1] CRAN (R 3.6.2)
     performance   0.7.2   2021-05-17 [1] CRAN (R 3.6.2)
     permute     * 0.9-5   2019-03-12 [1] CRAN (R 3.6.0)
     pillar        1.6.1   2021-05-16 [1] CRAN (R 3.6.2)
     pkgbuild      1.2.0   2020-12-15 [1] CRAN (R 3.6.2)
     pkgconfig     2.0.3   2019-09-22 [1] CRAN (R 3.6.0)
     pkgload       1.2.1   2021-04-06 [1] CRAN (R 3.6.2)
     prettyunits   1.1.1   2020-01-24 [1] CRAN (R 3.6.0)
     processx      3.5.2   2021-04-30 [1] CRAN (R 3.6.2)
     ps            1.6.0   2021-02-28 [1] CRAN (R 3.6.2)
     purrr         0.3.4   2020-04-17 [1] CRAN (R 3.6.2)
     R6            2.5.0   2020-10-28 [1] CRAN (R 3.6.2)
     Rcpp          1.0.6   2021-01-15 [1] CRAN (R 3.6.2)
     readxl        1.3.1   2019-03-13 [1] CRAN (R 3.6.0)
     remotes       2.4.0   2021-06-02 [1] CRAN (R 3.6.2)
     rio           0.5.26  2021-03-01 [1] CRAN (R 3.6.2)
     rlang         0.4.11  2021-04-30 [1] CRAN (R 3.6.2)
     rmarkdown     2.8     2021-05-07 [1] CRAN (R 3.6.2)
     rprojroot     2.0.2   2020-11-15 [1] CRAN (R 3.6.2)
     rstatix       0.7.0   2021-02-13 [1] CRAN (R 3.6.2)
     sass          0.4.0   2021-05-12 [1] CRAN (R 3.6.2)
     scales        1.1.1   2020-05-11 [1] CRAN (R 3.6.2)
     sessioninfo   1.1.1   2018-11-05 [1] CRAN (R 3.6.0)
     stringi       1.6.2   2021-05-17 [1] CRAN (R 3.6.2)
     stringr       1.4.0   2019-02-10 [1] CRAN (R 3.6.0)
     testthat      3.0.2   2021-02-14 [1] CRAN (R 3.6.2)
     tibble        3.1.2   2021-05-16 [1] CRAN (R 3.6.2)
     tidyr       * 1.1.3   2021-03-03 [1] CRAN (R 3.6.2)
     tidyselect    1.1.1   2021-04-30 [1] CRAN (R 3.6.2)
     usethis       2.0.1   2021-02-10 [1] CRAN (R 3.6.2)
     utf8          1.2.1   2021-03-12 [1] CRAN (R 3.6.2)
     vctrs         0.3.8   2021-04-29 [1] CRAN (R 3.6.2)
     vegan       * 2.5-7   2020-11-28 [1] CRAN (R 3.6.2)
     viridis     * 0.6.1   2021-05-11 [1] CRAN (R 3.6.2)
     viridisLite * 0.4.0   2021-04-13 [1] CRAN (R 3.6.2)
     withr         2.4.2   2021-04-18 [1] CRAN (R 3.6.2)
     xfun          0.23    2021-05-15 [1] CRAN (R 3.6.2)
     yaml          2.2.1   2020-02-01 [1] CRAN (R 3.6.0)
     zip           2.2.0   2021-05-31 [1] CRAN (R 3.6.2)
    
    [1] /Library/Frameworks/R.framework/Versions/3.6/Resources/library