library("ggplot2")
library("GGally")
library("Hmisc")
library("readr")
library("DESeq2")
library("phyloseq")
library('ggrepel')
library('ggalluvial')
library("metacoder")
library("microeco")
library("dplyr")
library("iNEXT")
library("networkD3")
library("tidyr")
library('tibble')
library('pals')
library("rstatix")
library("cowplot")
library("vegan")

sessionInfo()

## R version 4.5.0 (2025-04-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=English_United Kingdom.utf8 
## [2] LC_CTYPE=English_United Kingdom.utf8   
## [3] LC_MONETARY=English_United Kingdom.utf8
## [4] LC_NUMERIC=C                           
## [5] LC_TIME=English_United Kingdom.utf8    
## 
## time zone: Europe/Berlin
## tzcode source: internal
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] vegan_2.7-2                 permute_0.9-10             
##  [3] cowplot_1.2.0               rstatix_1.1.0              
##  [5] pals_1.10                   tibble_3.3.1               
##  [7] tidyr_1.3.2                 networkD3_0.4.1            
##  [9] iNEXT_3.0.2                 dplyr_1.2.1                
## [11] microeco_2.0.0              metacoder_0.3.8            
## [13] ggalluvial_0.12.5           ggrepel_0.9.6              
## [15] phyloseq_1.52.0             DESeq2_1.48.2              
## [17] SummarizedExperiment_1.38.1 Biobase_2.68.0             
## [19] MatrixGenerics_1.20.0       matrixStats_1.5.0          
## [21] GenomicRanges_1.60.0        GenomeInfoDb_1.44.3        
## [23] IRanges_2.42.0              S4Vectors_0.46.0           
## [25] BiocGenerics_0.54.1         generics_0.1.4             
## [27] readr_2.1.6                 Hmisc_5.2-6                
## [29] GGally_2.4.0                ggplot2_4.0.2              
## 
## loaded via a namespace (and not attached):
##  [1] gridExtra_2.3           rlang_1.2.0             magrittr_2.0.3         
##  [4] ade4_1.7-23             otel_0.2.0              compiler_4.5.0         
##  [7] mgcv_1.9-4              maps_3.4.3              vctrs_0.7.3            
## [10] reshape2_1.4.5          stringr_1.6.0           pkgconfig_2.0.3        
## [13] crayon_1.5.3            fastmap_1.2.0           backports_1.5.0        
## [16] XVector_0.48.0          rmarkdown_2.30          tzdb_0.5.0             
## [19] UCSC.utils_1.4.0        purrr_1.2.1             xfun_0.61              
## [22] cachem_1.1.0            jsonlite_2.0.0          biomformat_1.36.0      
## [25] rhdf5filters_1.20.0     DelayedArray_0.34.1     Rhdf5lib_1.30.0        
## [28] BiocParallel_1.42.2     broom_1.0.12            data.tree_1.2.0        
## [31] parallel_4.5.0          cluster_2.1.8.2         R6_2.6.1               
## [34] bslib_0.10.0            stringi_1.8.7           RColorBrewer_1.1-3     
## [37] car_3.1-5               rpart_4.1.24            jquerylib_0.1.4        
## [40] Rcpp_1.1.0              iterators_1.0.14        knitr_1.51             
## [43] base64enc_0.1-6         splines_4.5.0           Matrix_1.7-4           
## [46] nnet_7.3-20             igraph_2.2.2            tidyselect_1.2.1       
## [49] rstudioapi_0.18.0       dichromat_2.0-0.1       abind_1.4-8            
## [52] yaml_2.3.10             codetools_0.2-20        lattice_0.22-9         
## [55] plyr_1.8.9              withr_3.0.2             S7_0.2.1               
## [58] evaluate_1.0.5          foreign_0.8-91          survival_3.8-6         
## [61] ggstats_0.12.0          Biostrings_2.76.0       pillar_1.11.1          
## [64] carData_3.0-6           checkmate_2.3.4         foreach_1.5.2          
## [67] hms_1.1.4               scales_1.4.0            glue_1.8.0             
## [70] mapproj_1.2.12          tools_4.5.0             data.table_1.18.2.1    
## [73] locfit_1.5-9.12         rhdf5_2.52.1            grid_4.5.0             
## [76] ape_5.8-1               colorspace_2.1-2        nlme_3.1-168           
## [79] GenomeInfoDbData_1.2.14 htmlTable_2.4.3         Formula_1.2-5          
## [82] cli_3.6.5               S4Arrays_1.8.1          gtable_0.3.6           
## [85] sass_0.4.10             digest_0.6.37           SparseArray_1.8.1      
## [88] htmlwidgets_1.6.4       farver_2.1.2            htmltools_0.5.8.1      
## [91] multtest_2.64.0         lifecycle_1.0.5         httr_1.4.8             
## [94] MASS_7.3-65

sample_info <-read.table("Alaska_info.txt", row.names=1, header=T, check.names=F)
sample_info_no_R <- read.table("Alaska_info_no_R.txt", row.names=1, header=T, check.names=F)

tax_tab <- read_delim("Alaska_taxonomy.txt", delim = "\t", escape_double = FALSE, trim_ws = TRUE)

## New names:
## Rows: 467 Columns: 18
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "\t" chr
## (16): OTU, Species...4, ID1, ID2, reps/refs, Kingdom, Phylum, Class, Ord... dbl
## (2): size, match
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `Species` -> `Species...4`
## • `Species` -> `Species...14`

count_tab <-read.table("Alaska_counts.txt", header=T, row.names=1, check.names=F)
count_tab <- count_tab[rownames(count_tab) %in% tax_tab$OTU, ]

Treatment of NTC

NTC <-rowSums(count_tab[, 1:8]) # sum of NTC counts for each OTU
sample_counts <- rowSums(count_tab[, 9:143]) # sum of sample counts for each OTU

# normalize sample counts by dividing the samples' total by 17 – as there are 17 as many samples (135) as there are NTCs (8) (135/8=17)
norm_sample_OTU_counts <- sample_counts/17

# which OTUs are deemed likely contaminants based on the threshold noted above:
blank_OTUs <- names(NTC[NTC *10 > norm_sample_OTU_counts]) # OTUs which have more counts than 1/10 of the sample counts
length(blank_OTUs)

## [1] 4

# How many counts are kept
colSums(count_tab[!rownames(count_tab) %in% blank_OTUs, ]) / colSums(count_tab) * 100

##        NC3        NC1        NC2       NEG1       NEG2        NC4       NTC1 
##  100.00000   10.64372   77.06767        NaN        NaN  100.00000  100.00000 
##       NTC2       BT38      45207    DF46144      BF193      PE048       BT51 
##  100.00000   99.32096   94.24296   99.57122   99.47561   99.75702   98.92461 
##     BT25-1 DF45570-r4    DF45955      BF187    PE168-1      BF234   BT118-r4 
##   97.84747   98.95379   98.68708   99.93803   99.49115   97.23105   97.45719 
##      45261    DF46114    DF46025      45263       BT20   45252-r2      PE455 
##   90.18482   98.56276   99.83654   91.27261   92.07628   94.80401   98.85126 
##   45252-r5      PE184      PE062    DF46110      45224      PE043      PE475 
##   94.69995   99.79756   99.92551   99.83511   96.55938   99.81288   99.55835 
##      45375    DF46112   BT118-r1    DF45989 DF45570-r2      BF188    PE168-3 
##   98.09209   99.51626   99.99434   99.60050   97.65897   96.58367   99.99696 
##    DF46090       BT35   45252-r3      PE411      BT104      BT114    BF278-2 
##   99.80990   99.44001   90.69389   99.95844   98.69873   99.99720   99.24072 
##    DF46099       BT33       BT70      BF182    DF46159      45418      45311 
##   99.73984   99.98684   99.63897   97.73173   99.97253   98.33936   96.71749 
##      BF255       BT30 DF45570-r3   BT118-r2      BF252      BF340    PE168-2 
##   98.78447   98.72714   75.00000   99.42968   97.12428   94.50588  100.00000 
##      45366       BT62      45368      BF153      45504      45351   45252-r4 
##   95.74179   99.31377   97.09630   99.51104   85.48662   90.26309   98.01217 
##    DF45701      45267      45293      45257      45252      45234      45302 
##   99.32073   97.91928   97.00096   96.33471   89.66172   97.69718   93.29878 
##      BF165      BT121      PE106      BT118      PE116      BF149    DF45698 
##   83.02010   99.88608   99.85454   97.69966   99.98451   99.84160   97.43421 
##    45209-1    DF45570    DF46001      BF190   BT118-r3    45209-2       BT31 
##   92.75857   98.58597   99.15343   99.83561   99.68251   92.86548   99.50164 
##      PE176      45533      BF324   BT118-r5      PE134      45354      45430 
##   99.81732   94.52863   97.88887   94.82207   99.52614   89.88541   94.85348 
##      45381      PE130      45465      PE058      45313      BF346      PE162 
##   88.89717   99.99321   94.83919  100.00000   96.85421   98.28296  100.00000 
##      BT100      BT105      45329        BT5    BF278-3      BF277      45393 
##   99.01163   98.30174   95.72949   95.94683   99.16105   99.72295   96.72852 
##      PE160    DF45979     BT25-2       BT57      PE042 DF45570-r5    DF45944 
##   99.28195   98.57165   97.17837   98.32061   99.52242   99.76227   99.23664 
##   45252-r1    45209-3    DF45567      BF283      45516       BT43      45346 
##   94.08105   92.36757   99.98362   99.83531   97.56185   97.46956   89.95767 
##      PE038      PE044    DF46039        BT4    DF46162      BF191      45442 
##   99.92056   99.60134   98.99022   97.78467   95.29840   96.24360   96.45190 
##      PE164      BF300      PE035 DF45570-r1     BT25-3      45201      BF202 
##   99.99452   93.14790   99.98656   92.93257   97.68558   90.31334   99.40860 
##      PE045      BF280    BF278-1 
##  100.00000   84.89532   99.23334

filt_count_tab <- count_tab[!rownames(count_tab) %in% blank_OTUs, -c(1:8)] # removing blank_OTUs and the blank samples from further analysis
filt_count_tab <- filt_count_tab[rowSums(filt_count_tab == 0) != ncol(filt_count_tab), ] # delete OTUs with 0 counts

Normalizing for sampling depth (variance stabilizing transformation)

sample_info<-sample_info[order(rownames(sample_info)), ]
filt_count_tab1<- filt_count_tab[,order(names(filt_count_tab))]

deseq_counts <- DESeqDataSetFromMatrix(filt_count_tab1, colData = sample_info, design = ~ Site)

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

deseq_counts <- estimateSizeFactors(deseq_counts,type = "poscounts")
deseq_counts_vst <- varianceStabilizingTransformation(deseq_counts,blind = TRUE)
vst_trans_count_tab <- assay(deseq_counts_vst) #save as matrix
vst_trans_count_tab_df <- as.data.frame(vst_trans_count_tab)

PCoA (Plot)

palette <- c("#EE7733","#CC6677", "#009988","#337538","#33BBEE", "#5E4FA2")

# create phyloseq object with normalized data
vst_count_phy <- otu_table(vst_trans_count_tab, taxa_are_rows=T)
sample_info_tab_phy <- sample_data(sample_info)
vst_physeq <- phyloseq(vst_count_phy, sample_info_tab_phy)

#generating and visualizing the PCoA with phyloseq
vst_pcoa <- ordinate(vst_physeq, method="PCoA", distance="euclidean")

eigen_vals <- vst_pcoa$values$Eigenvalues # to scale the axes
total_var <- sum(eigen_vals)
axis1 <- 100 * eigen_vals[1] / total_var
axis2 <- 100 * eigen_vals[2] / total_var

plot_ordination(vst_physeq, vst_pcoa, color = "Plot", shape = "Habitat") +
  geom_point(size = 2) +
  geom_text_repel(aes(label = rownames(sample_info)), size = 2, force = 1, max.overlaps = 100) +
  stat_ellipse(aes(fill = Plot), geom = "polygon", color = NA, alpha = 0.1) +
  coord_fixed(sqrt(eigen_vals[2] / eigen_vals[1])) +
    labs(x = paste0("Principal Coordinate Axis 1 (", round(axis1, 2), "%)"),
    y = paste0("Principal Coordinate Axis 2 (", round(axis2, 2), "%)")) +
  theme_classic() +
  scale_shape_manual(values = c(20, 17)) +
  scale_color_manual(values = palette) +
  scale_fill_manual(values = palette)

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the phyloseq package.
##   Please report the issue at <https://github.com/joey711/phyloseq/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Without technical replicates

Normalizing for sampling depth (variance stabilizing transformation)

#reorder the row and column names
sample_info_no_R<-sample_info_no_R[order(rownames(sample_info_no_R)), ]
filt_count_tab<- filt_count_tab[,order(names(filt_count_tab))]
filt_count_tab_no_r <- filt_count_tab[colnames(filt_count_tab) %in% rownames(sample_info_no_R) ]

# create a DESeq2 object
deseq_counts <- DESeqDataSetFromMatrix(filt_count_tab_no_r, colData = sample_info_no_R, design = ~ Habitat)

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

deseq_counts <- estimateSizeFactors(deseq_counts,type = "poscounts")
deseq_counts_vst <- varianceStabilizingTransformation(deseq_counts,blind = TRUE)
vst_trans_count_tab <- assay(deseq_counts_vst) #save as matrix

PCoA (Plot, All OTUs)

# create phyloseq object with normalized data
vst_count_phy <- otu_table(vst_trans_count_tab, taxa_are_rows=T)
sample_info_tab_phy <- sample_data(sample_info_no_R)
vst_physeq <- phyloseq(vst_count_phy, sample_info_tab_phy)

#generating and visualizing the PCoA with phyloseq
vst_pcoa <- ordinate(vst_physeq, method="PCoA", distance="euclidean")

eigen_vals <- vst_pcoa$values$Eigenvalues # ato scale the axes
total_var <- sum(eigen_vals)
axis1 <- 100 * eigen_vals[1] / total_var
axis2 <- 100 * eigen_vals[2] / total_var

plot_ordination(vst_physeq, vst_pcoa, color = "Plot", shape = "Habitat") +
  geom_point(size = 2, alpha = 0.5) +
  stat_ellipse(aes(fill = Plot), geom = "polygon", color = NA, alpha = 0.2) +
  coord_fixed(sqrt(eigen_vals[2] / eigen_vals[1])) +
    labs(x = paste0("Principal Coordinate Axis 1 (", round(axis1, 2), "%)"),
    y = paste0("Principal Coordinate Axis 2 (", round(axis2, 2), "%)")) +
  theme_classic() +
  scale_shape_manual(values = c(20, 17)) +
  scale_color_manual(values = palette) +
  scale_fill_manual(values = palette)

PCoA (Plot, only ECM)

ECM_OTUs <- tax_tab$OTU[tax_tab$Guild %in% c("Ectomycorrhizal")]
length(ECM_OTUs)

## [1] 247

Sapro_OTU<-tax_tab$OTU[tax_tab$Guild %in% c("Wood Saprotroph","Plant Saprotroph","Undefined Saprotroph", "Plant Pathogen")]
length(Sapro_OTU)

## [1] 147

ECM_OTUs_present <- intersect(ECM_OTUs, rownames(vst_trans_count_tab))

# create phyloseq object with normalized data
vst_count_phy_ECM <- otu_table(vst_trans_count_tab[ECM_OTUs_present, ], taxa_are_rows=T)
sample_info_tab_phy_ECM <- sample_data(sample_info_no_R)
vst_physeq_ECM <- phyloseq(vst_count_phy_ECM, sample_info_tab_phy)

#generating and visualizing the PCoA with phyloseq
vst_pcoa_ECM <- ordinate(vst_physeq_ECM, method="PCoA", distance="euclidean")

eigen_vals_ECM <- vst_pcoa_ECM$values$Eigenvalues # to scale the axes
total_var <- sum(eigen_vals_ECM)
axis1 <- 100 * eigen_vals_ECM[1] / total_var
axis2 <- 100 * eigen_vals_ECM[2] / total_var

plot_ordination(vst_physeq_ECM, vst_pcoa_ECM, color = "Plot", shape = "Habitat") +
  geom_point(size = 2, alpha = 0.5) +
  stat_ellipse(aes(fill = Plot), geom = "polygon", color = NA, alpha = 0.2) +
  coord_fixed(sqrt(eigen_vals_ECM[2] / eigen_vals_ECM[1])) +
    labs(x = paste0("Principal Coordinate Axis 1 (", round(axis1, 2), "%)"),
    y = paste0("Principal Coordinate Axis 2 (", round(axis2, 2), "%)")) +
  theme_classic() +
  scale_shape_manual(values = c(20, 17)) +
  scale_color_manual(values = palette) +
  scale_fill_manual(values = palette)

## Warning in MASS::cov.trob(data[, vars], wt = weight * nrow(data)): Probable
## convergence failure

Relative Abundancy

tax_table<-as.data.frame(tax_tab)
row.names(tax_table) <- tax_table$OTU

tax_table <- tax_table[, -1]
tax_table<- tax_table[order(rownames(tax_table)),]

# create microeco file
meco_fungi <- microtable$new(sample_table = sample_info_no_R, otu_table = filt_count_tab, tax_table = tax_table)
 
mycol <- c("#A6CEE3" ,"#1F78B4" ,"#B2DF8A" ,"#33A02C" ,"#FB9A99", "#E31A1C" ,"#FDBF6F" ,"#FF7F00" ,"#CAB2D6" ,
           "#6A3D9A" ,"#FFFF99" ,"#B15928" , "#6fdc8c", "green4", "#004144")

t1 <- trans_abund$new(dataset = meco_fungi, taxrank = "Phylum", ntaxa = 6)

## No taxa_abund list found. Calculate relative abundance with cal_abund function ...

## The result is stored in object$taxa_abund ...

## The transformed abundance data is stored in object$data_abund ...

t2 <- trans_abund$new(dataset = meco_fungi, taxrank = "Class", ntaxa = 11)

## The transformed abundance data is stored in object$data_abund ...

t3 <- trans_abund$new(dataset = meco_fungi, taxrank = "Order", ntaxa = 11)

## The transformed abundance data is stored in object$data_abund ...

t4 <- trans_abund$new(dataset = meco_fungi, taxrank = "Family", ntaxa = 11)

## The transformed abundance data is stored in object$data_abund ...

t5 <- trans_abund$new(dataset = meco_fungi, taxrank = "Genus", ntaxa = 12)

## The transformed abundance data is stored in object$data_abund ...

p1<-t1$plot_bar(xtext_keep = TRUE, xtext_angle = 90,facet = "Plot",xtext_size = 6, barwidth = 1, color_values = mycol)
p2<-t2$plot_bar(xtext_keep = TRUE, xtext_angle = 90,facet = "Plot",xtext_size = 6, barwidth = 1,color_values = mycol)
p3<-t3$plot_bar(xtext_keep = TRUE, xtext_angle = 90,facet = "Plot",xtext_size = 6, barwidth = 1,color_values = mycol)
p4<-t4$plot_bar(xtext_keep = TRUE, xtext_angle = 90,facet = "Plot",xtext_size = 6, barwidth = 1,color_values = mycol)

library(gridExtra)

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

## The following object is masked from 'package:Biobase':
## 
##     combine

## The following object is masked from 'package:BiocGenerics':
## 
##     combine

grid.arrange(p1, p2, p3, p4, ncol = 1)

Relative abundane - Genus

t5$plot_bar(xtext_keep = TRUE, xtext_angle = 90,facet = "Plot",xtext_size = 6, barwidth = 1,color_values = mycol)

Prediction of the RAF composition (all OTUs): Variance partitioning

meta <- as(sample_data(vst_physeq), "data.frame")
keep <- complete.cases(meta[, c("Age","pH","Habitat","Site")])
vst_mat <- as.data.frame(t(otu_table(vst_physeq)))
vst_mat_clean <- vst_mat[keep, ]
meta_clean <- meta[keep, ]

# List of predictors
env <- meta_clean[, c("pH", "Habitat")]
env$Habitat <- as.factor(env$Habitat)
space <- data.frame(Site = as.factor(meta_clean$Site))
growth <- meta_clean$Age

# Variance partitioning

vp <- varpart(vst_mat_clean, ~ pH, ~ Habitat, ~ Site, ~ Age, data = meta_clean)
vals <- vp$part$indfract$Adj.R.square[1:4]

plot(vp, Xnames = c("pH", "Habitat","Site",  "Age"), bg = c("#EE7733","#5E4FA2","#337538" ,"#33BBEE"), cex = 1, digits = 2)

Redundance analyses - All OTUs

rda_env <- rda(vst_mat_clean ~ pH + Habitat + Condition(Site), data = meta_clean)
anova.cca(rda_env, by = "term", permutations = 9999)

## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
## 
## Model: rda(formula = vst_mat_clean ~ pH + Habitat + Condition(Site), data = meta_clean)
##           Df Variance      F Pr(>F)    
## pH         1    159.7 5.0912  1e-04 ***
## Habitat    1     74.5 2.3754  4e-04 ***
## Residual 103   3230.2                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

rda_space <- rda(vst_mat_clean ~ Site, data = meta_clean)
anova.cca(rda_space, permutations = 9999)

## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 9999
## 
## Model: rda(formula = vst_mat_clean ~ Site, data = meta_clean)
##           Df Variance      F Pr(>F)    
## Model      2    867.0 13.139  1e-04 ***
## Residual 105   3464.4                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

rda_age <- rda(vst_mat_clean ~ Age + Condition(pH + Habitat + Site), data = meta_clean)
anova.cca(rda_age, by = "term", permutations = 9999)

## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
## 
## Model: rda(formula = vst_mat_clean ~ Age + Condition(pH + Habitat + Site), data = meta_clean)
##           Df Variance      F Pr(>F)   
## Age        1     59.3 1.9072 0.0011 **
## Residual 102   3170.9                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Prediction of the RAF composition (only ECM OTUs): Variance partitioning

meta <- as(sample_data(vst_physeq), "data.frame")
keep <- complete.cases(meta[, c("Age","pH","Habitat","Site")])
vst_mat <- as.data.frame(t(otu_table(vst_physeq)))
vst_mat_clean <- vst_mat[keep, ]
meta_clean <- meta[keep, ]

tax_data <- tax_tab
tax_data_ECM <- tax_data[tax_data$Guild == "Ectomycorrhizal", ]
vst_mat_clean_ECM <- vst_mat_clean[, colnames(vst_mat_clean) %in% tax_data_ECM$OTU]


# List of predictors
env <- meta_clean[, c("pH", "Habitat")]
env$Habitat <- as.factor(env$Habitat)
space <- data.frame(Site = as.factor(meta_clean$Site))
growth <- meta_clean$Age

# Variance partitioning

vp <- varpart(vst_mat_clean_ECM, ~ pH, ~ Habitat, ~ Site, ~ Age, data = meta_clean)
vals <- vp$part$indfract$Adj.R.square[1:4]

plot(vp, Xnames = c("pH", "Habitat","Site",  "Age"), bg = c("#EE7733","#5E4FA2","#337538" ,"#33BBEE"), cex = 1, digits = 2)

Redundance analyses - Only ECM

rda_env <- rda(vst_mat_clean_ECM ~ pH + Habitat + Condition(Site), data = meta_clean)
anova.cca(rda_env, by = "term", permutations = 9999)

## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
## 
## Model: rda(formula = vst_mat_clean_ECM ~ pH + Habitat + Condition(Site), data = meta_clean)
##           Df Variance      F Pr(>F)    
## pH         1    68.88 3.8858 0.0001 ***
## Habitat    1    31.64 1.7853 0.0047 ** 
## Residual 103  1825.67                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

rda_space <- rda(vst_mat_clean_ECM ~ Site, data = meta_clean)
anova.cca(rda_space, permutations = 9999)

## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 9999
## 
## Model: rda(formula = vst_mat_clean_ECM ~ Site, data = meta_clean)
##           Df Variance      F Pr(>F)    
## Model      2   367.73 10.023  1e-04 ***
## Residual 105  1926.19                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

rda_age <- rda(vst_mat_clean_ECM ~ Age + Condition(pH + Habitat + Site), data = meta_clean)
anova.cca(rda_age, by = "term", permutations = 9999)

## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
## 
## Model: rda(formula = vst_mat_clean_ECM ~ Age + Condition(pH + Habitat + Site), data = meta_clean)
##           Df Variance     F Pr(>F)  
## Age        1    25.28 1.432 0.0386 *
## Residual 102  1800.40               
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Mantel test (all OTUs)

# test whether RAF community dissimilarity increases with the distance between
# trees inside a plot (distance decay), and whether such a pattern is explained
# by differences in soil pH between the trees

plots <- sort(unique(as.character(meta_clean$Plot)))   # AF AT BF BT IF IT

set.seed(50)

# ---------- AF
sel <- which(as.character(meta_clean$Plot) == "AF" &
               !is.na(meta_clean$Latitude) & !is.na(meta_clean$Longitude))
meta_plot <- meta_clean[sel, ]
lat_mean  <- mean(meta_plot$Latitude)
x <- (meta_plot$Longitude - mean(meta_plot$Longitude)) * 111320 * cos(lat_mean * pi / 180)
y <- (meta_plot$Latitude - lat_mean) * 110540
dist_geo <- dist(cbind(x, y))
dist_com <- dist(vst_mat_clean[sel, ])
dist_pH  <- dist(as.matrix(meta_plot$pH))
mt_geo    <- mantel(dist_com, dist_geo, method = "pearson", permutations = 9999)
mt_geo_pH <- mantel.partial(dist_com, dist_geo, dist_pH, method = "pearson", permutations = 9999)
mt_pH_geo <- mantel(dist_pH, dist_geo, method = "pearson", permutations = 9999)
m_af <- data.frame(Plot = "AF", n = length(sel),
                r_geo    = mt_geo$statistic,    p_geo    = mt_geo$signif,
                r_geo_pH = mt_geo_pH$statistic, p_geo_pH = mt_geo_pH$signif,
                r_pH_geo = mt_pH_geo$statistic, p_pH_geo = mt_pH_geo$signif)

# ---------- AT
sel <- which(as.character(meta_clean$Plot) == "AT" &
               !is.na(meta_clean$Latitude) & !is.na(meta_clean$Longitude))
meta_plot <- meta_clean[sel, ]
lat_mean  <- mean(meta_plot$Latitude)
x <- (meta_plot$Longitude - mean(meta_plot$Longitude)) * 111320 * cos(lat_mean * pi / 180)
y <- (meta_plot$Latitude - lat_mean) * 110540
dist_geo <- dist(cbind(x, y))
dist_com <- dist(vst_mat_clean[sel, ])
dist_pH  <- dist(as.matrix(meta_plot$pH))
mt_geo    <- mantel(dist_com, dist_geo, method = "pearson", permutations = 9999)
mt_geo_pH <- mantel.partial(dist_com, dist_geo, dist_pH, method = "pearson", permutations = 9999)
mt_pH_geo <- mantel(dist_pH, dist_geo, method = "pearson", permutations = 9999)
m_at <- data.frame(Plot = "AT", n = length(sel),
                r_geo    = mt_geo$statistic,    p_geo    = mt_geo$signif,
                r_geo_pH = mt_geo_pH$statistic, p_geo_pH = mt_geo_pH$signif,
                r_pH_geo = mt_pH_geo$statistic, p_pH_geo = mt_pH_geo$signif)

# ---------- BF
sel <- which(as.character(meta_clean$Plot) == "BF" &
               !is.na(meta_clean$Latitude) & !is.na(meta_clean$Longitude))
meta_plot <- meta_clean[sel, ]
lat_mean  <- mean(meta_plot$Latitude)
x <- (meta_plot$Longitude - mean(meta_plot$Longitude)) * 111320 * cos(lat_mean * pi / 180)
y <- (meta_plot$Latitude - lat_mean) * 110540
dist_geo <- dist(cbind(x, y))
dist_com <- dist(vst_mat_clean[sel, ])
dist_pH  <- dist(as.matrix(meta_plot$pH))
mt_geo    <- mantel(dist_com, dist_geo, method = "pearson", permutations = 9999)
mt_geo_pH <- mantel.partial(dist_com, dist_geo, dist_pH, method = "pearson", permutations = 9999)
mt_pH_geo <- mantel(dist_pH, dist_geo, method = "pearson", permutations = 9999)
m_bf <- data.frame(Plot = "BF", n = length(sel),
                r_geo    = mt_geo$statistic,    p_geo    = mt_geo$signif,
                r_geo_pH = mt_geo_pH$statistic, p_geo_pH = mt_geo_pH$signif,
                r_pH_geo = mt_pH_geo$statistic, p_pH_geo = mt_pH_geo$signif)

# ---------- BT
sel <- which(as.character(meta_clean$Plot) == "BT" &
               !is.na(meta_clean$Latitude) & !is.na(meta_clean$Longitude))
meta_plot <- meta_clean[sel, ]
lat_mean  <- mean(meta_plot$Latitude)
x <- (meta_plot$Longitude - mean(meta_plot$Longitude)) * 111320 * cos(lat_mean * pi / 180)
y <- (meta_plot$Latitude - lat_mean) * 110540
dist_geo <- dist(cbind(x, y))
dist_com <- dist(vst_mat_clean[sel, ])
dist_pH  <- dist(as.matrix(meta_plot$pH))
mt_geo    <- mantel(dist_com, dist_geo, method = "pearson", permutations = 9999)
mt_geo_pH <- mantel.partial(dist_com, dist_geo, dist_pH, method = "pearson", permutations = 9999)
mt_pH_geo <- mantel(dist_pH, dist_geo, method = "pearson", permutations = 9999)
m_bt <- data.frame(Plot = "BT", n = length(sel),
                r_geo    = mt_geo$statistic,    p_geo    = mt_geo$signif,
                r_geo_pH = mt_geo_pH$statistic, p_geo_pH = mt_geo_pH$signif,
                r_pH_geo = mt_pH_geo$statistic, p_pH_geo = mt_pH_geo$signif)

# ---------- IF
sel <- which(as.character(meta_clean$Plot) == "IF" &
               !is.na(meta_clean$Latitude) & !is.na(meta_clean$Longitude))
meta_plot <- meta_clean[sel, ]
lat_mean  <- mean(meta_plot$Latitude)
x <- (meta_plot$Longitude - mean(meta_plot$Longitude)) * 111320 * cos(lat_mean * pi / 180)
y <- (meta_plot$Latitude - lat_mean) * 110540
dist_geo <- dist(cbind(x, y))
dist_com <- dist(vst_mat_clean[sel, ])
dist_pH  <- dist(as.matrix(meta_plot$pH))
mt_geo    <- mantel(dist_com, dist_geo, method = "pearson", permutations = 9999)
mt_geo_pH <- mantel.partial(dist_com, dist_geo, dist_pH, method = "pearson", permutations = 9999)
mt_pH_geo <- mantel(dist_pH, dist_geo, method = "pearson", permutations = 9999)
m_if <- data.frame(Plot = "IF", n = length(sel),
                r_geo    = mt_geo$statistic,    p_geo    = mt_geo$signif,
                r_geo_pH = mt_geo_pH$statistic, p_geo_pH = mt_geo_pH$signif,
                r_pH_geo = mt_pH_geo$statistic, p_pH_geo = mt_pH_geo$signif)

# ---------- IT
sel <- which(as.character(meta_clean$Plot) == "IT" &
               !is.na(meta_clean$Latitude) & !is.na(meta_clean$Longitude))
meta_plot <- meta_clean[sel, ]
lat_mean  <- mean(meta_plot$Latitude)
x <- (meta_plot$Longitude - mean(meta_plot$Longitude)) * 111320 * cos(lat_mean * pi / 180)
y <- (meta_plot$Latitude - lat_mean) * 110540
dist_geo <- dist(cbind(x, y))
dist_com <- dist(vst_mat_clean[sel, ])
dist_pH  <- dist(as.matrix(meta_plot$pH))
mt_geo    <- mantel(dist_com, dist_geo, method = "pearson", permutations = 9999)
mt_geo_pH <- mantel.partial(dist_com, dist_geo, dist_pH, method = "pearson", permutations = 9999)
mt_pH_geo <- mantel(dist_pH, dist_geo, method = "pearson", permutations = 9999)
m_it <- data.frame(Plot = "IT", n = length(sel),
                r_geo    = mt_geo$statistic,    p_geo    = mt_geo$signif,
                r_geo_pH = mt_geo_pH$statistic, p_geo_pH = mt_geo_pH$signif,
                r_pH_geo = mt_pH_geo$statistic, p_pH_geo = mt_pH_geo$signif)

mantel_res <- rbind(m_af, m_at, m_bf, m_bt, m_if, m_it)

# adjust the p-values across the six plots
mantel_res$p_geo_adj    <- p.adjust(mantel_res$p_geo,    method = "BH")
mantel_res$p_geo_pH_adj <- p.adjust(mantel_res$p_geo_pH, method = "BH")
mantel_res$p_pH_geo_adj <- p.adjust(mantel_res$p_pH_geo, method = "BH")

mantel_res

##   Plot  n       r_geo  p_geo    r_geo_pH p_geo_pH     r_pH_geo p_pH_geo
## 1   AF 19 -0.10913090 0.8051 -0.11086896   0.7990 -0.051944461   0.6309
## 2   AT 18  0.01586464 0.4168  0.01544612   0.4300  0.007359699   0.4622
## 3   BF 21  0.01627938 0.3815  0.03014043   0.3383 -0.024435790   0.5790
## 4   BT 20 -0.17570682 0.9821 -0.18561676   0.9859 -0.034464525   0.5789
## 5   IF 19 -0.12709225 0.8341 -0.13691800   0.8586  0.052928550   0.2925
## 6   IT 11 -0.06342222 0.6414 -0.06976820   0.6512 -0.200568892   0.9408
##   p_geo_adj p_geo_pH_adj p_pH_geo_adj
## 1    0.9821       0.9859      0.75708
## 2    0.9821       0.9859      0.75708
## 3    0.9821       0.9859      0.75708
## 4    0.9821       0.9859      0.75708
## 5    0.9821       0.9859      0.75708
## 6    0.9821       0.9859      0.94080

Mantel test (ECM OTUs)

# test whether RAF community dissimilarity increases with the distance between
# trees inside a plot (distance decay), and whether such a pattern is explained
# by differences in soil pH between the trees

plots <- sort(unique(as.character(meta_clean$Plot)))   # AF AT BF BT IF IT

set.seed(50)

# ---------- AF
sel <- which(as.character(meta_clean$Plot) == "AF" &
               !is.na(meta_clean$Latitude) & !is.na(meta_clean$Longitude))
meta_plot <- meta_clean[sel, ]
lat_mean  <- mean(meta_plot$Latitude)
x <- (meta_plot$Longitude - mean(meta_plot$Longitude)) * 111320 * cos(lat_mean * pi / 180)
y <- (meta_plot$Latitude - lat_mean) * 110540
dist_geo <- dist(cbind(x, y))
dist_com <- dist(vst_mat_clean_ECM[sel, ])
dist_pH  <- dist(as.matrix(meta_plot$pH))
mt_geo    <- mantel(dist_com, dist_geo, method = "pearson", permutations = 9999)
mt_geo_pH <- mantel.partial(dist_com, dist_geo, dist_pH, method = "pearson", permutations = 9999)
mt_pH_geo <- mantel(dist_pH, dist_geo, method = "pearson", permutations = 9999)
m_af <- data.frame(Plot = "AF", n = length(sel),
                r_geo    = mt_geo$statistic,    p_geo    = mt_geo$signif,
                r_geo_pH = mt_geo_pH$statistic, p_geo_pH = mt_geo_pH$signif,
                r_pH_geo = mt_pH_geo$statistic, p_pH_geo = mt_pH_geo$signif)

# ---------- AT
sel <- which(as.character(meta_clean$Plot) == "AT" &
               !is.na(meta_clean$Latitude) & !is.na(meta_clean$Longitude))
meta_plot <- meta_clean[sel, ]
lat_mean  <- mean(meta_plot$Latitude)
x <- (meta_plot$Longitude - mean(meta_plot$Longitude)) * 111320 * cos(lat_mean * pi / 180)
y <- (meta_plot$Latitude - lat_mean) * 110540
dist_geo <- dist(cbind(x, y))
dist_com <- dist(vst_mat_clean_ECM[sel, ])
dist_pH  <- dist(as.matrix(meta_plot$pH))
mt_geo    <- mantel(dist_com, dist_geo, method = "pearson", permutations = 9999)
mt_geo_pH <- mantel.partial(dist_com, dist_geo, dist_pH, method = "pearson", permutations = 9999)
mt_pH_geo <- mantel(dist_pH, dist_geo, method = "pearson", permutations = 9999)
m_at <- data.frame(Plot = "AT", n = length(sel),
                r_geo    = mt_geo$statistic,    p_geo    = mt_geo$signif,
                r_geo_pH = mt_geo_pH$statistic, p_geo_pH = mt_geo_pH$signif,
                r_pH_geo = mt_pH_geo$statistic, p_pH_geo = mt_pH_geo$signif)

# ---------- BF
sel <- which(as.character(meta_clean$Plot) == "BF" &
               !is.na(meta_clean$Latitude) & !is.na(meta_clean$Longitude))
meta_plot <- meta_clean[sel, ]
lat_mean  <- mean(meta_plot$Latitude)
x <- (meta_plot$Longitude - mean(meta_plot$Longitude)) * 111320 * cos(lat_mean * pi / 180)
y <- (meta_plot$Latitude - lat_mean) * 110540
dist_geo <- dist(cbind(x, y))
dist_com <- dist(vst_mat_clean_ECM[sel, ])
dist_pH  <- dist(as.matrix(meta_plot$pH))
mt_geo    <- mantel(dist_com, dist_geo, method = "pearson", permutations = 9999)
mt_geo_pH <- mantel.partial(dist_com, dist_geo, dist_pH, method = "pearson", permutations = 9999)
mt_pH_geo <- mantel(dist_pH, dist_geo, method = "pearson", permutations = 9999)
m_bf <- data.frame(Plot = "BF", n = length(sel),
                r_geo    = mt_geo$statistic,    p_geo    = mt_geo$signif,
                r_geo_pH = mt_geo_pH$statistic, p_geo_pH = mt_geo_pH$signif,
                r_pH_geo = mt_pH_geo$statistic, p_pH_geo = mt_pH_geo$signif)

# ---------- BT
sel <- which(as.character(meta_clean$Plot) == "BT" &
               !is.na(meta_clean$Latitude) & !is.na(meta_clean$Longitude))
meta_plot <- meta_clean[sel, ]
lat_mean  <- mean(meta_plot$Latitude)
x <- (meta_plot$Longitude - mean(meta_plot$Longitude)) * 111320 * cos(lat_mean * pi / 180)
y <- (meta_plot$Latitude - lat_mean) * 110540
dist_geo <- dist(cbind(x, y))
dist_com <- dist(vst_mat_clean_ECM[sel, ])
dist_pH  <- dist(as.matrix(meta_plot$pH))
mt_geo    <- mantel(dist_com, dist_geo, method = "pearson", permutations = 9999)
mt_geo_pH <- mantel.partial(dist_com, dist_geo, dist_pH, method = "pearson", permutations = 9999)
mt_pH_geo <- mantel(dist_pH, dist_geo, method = "pearson", permutations = 9999)
m_bt <- data.frame(Plot = "BT", n = length(sel),
                r_geo    = mt_geo$statistic,    p_geo    = mt_geo$signif,
                r_geo_pH = mt_geo_pH$statistic, p_geo_pH = mt_geo_pH$signif,
                r_pH_geo = mt_pH_geo$statistic, p_pH_geo = mt_pH_geo$signif)

# ---------- IF
sel <- which(as.character(meta_clean$Plot) == "IF" &
               !is.na(meta_clean$Latitude) & !is.na(meta_clean$Longitude))
meta_plot <- meta_clean[sel, ]
lat_mean  <- mean(meta_plot$Latitude)
x <- (meta_plot$Longitude - mean(meta_plot$Longitude)) * 111320 * cos(lat_mean * pi / 180)
y <- (meta_plot$Latitude - lat_mean) * 110540
dist_geo <- dist(cbind(x, y))
dist_com <- dist(vst_mat_clean_ECM[sel, ])
dist_pH  <- dist(as.matrix(meta_plot$pH))
mt_geo    <- mantel(dist_com, dist_geo, method = "pearson", permutations = 9999)
mt_geo_pH <- mantel.partial(dist_com, dist_geo, dist_pH, method = "pearson", permutations = 9999)
mt_pH_geo <- mantel(dist_pH, dist_geo, method = "pearson", permutations = 9999)
m_if <- data.frame(Plot = "IF", n = length(sel),
                r_geo    = mt_geo$statistic,    p_geo    = mt_geo$signif,
                r_geo_pH = mt_geo_pH$statistic, p_geo_pH = mt_geo_pH$signif,
                r_pH_geo = mt_pH_geo$statistic, p_pH_geo = mt_pH_geo$signif)

# ---------- IT
sel <- which(as.character(meta_clean$Plot) == "IT" &
               !is.na(meta_clean$Latitude) & !is.na(meta_clean$Longitude))
meta_plot <- meta_clean[sel, ]
lat_mean  <- mean(meta_plot$Latitude)
x <- (meta_plot$Longitude - mean(meta_plot$Longitude)) * 111320 * cos(lat_mean * pi / 180)
y <- (meta_plot$Latitude - lat_mean) * 110540
dist_geo <- dist(cbind(x, y))
dist_com <- dist(vst_mat_clean_ECM[sel, ])
dist_pH  <- dist(as.matrix(meta_plot$pH))
mt_geo    <- mantel(dist_com, dist_geo, method = "pearson", permutations = 9999)
mt_geo_pH <- mantel.partial(dist_com, dist_geo, dist_pH, method = "pearson", permutations = 9999)
mt_pH_geo <- mantel(dist_pH, dist_geo, method = "pearson", permutations = 9999)
m_it <- data.frame(Plot = "IT", n = length(sel),
                r_geo    = mt_geo$statistic,    p_geo    = mt_geo$signif,
                r_geo_pH = mt_geo_pH$statistic, p_geo_pH = mt_geo_pH$signif,
                r_pH_geo = mt_pH_geo$statistic, p_pH_geo = mt_pH_geo$signif)

mantel_res_ECM <- rbind(m_af, m_at, m_bf, m_bt, m_if, m_it)

# adjust the p-values across the six plots
mantel_res_ECM$p_geo_adj    <- p.adjust(mantel_res_ECM$p_geo,    method = "BH")
mantel_res_ECM$p_geo_pH_adj <- p.adjust(mantel_res_ECM$p_geo_pH, method = "BH")
mantel_res_ECM$p_pH_geo_adj <- p.adjust(mantel_res_ECM$p_pH_geo, method = "BH")

mantel_res_ECM

##   Plot  n       r_geo  p_geo    r_geo_pH p_geo_pH     r_pH_geo p_pH_geo
## 1   AF 19 -0.15294587 0.9032 -0.15711127   0.9028 -0.051944461   0.6309
## 2   AT 18 -0.01924579 0.5235 -0.01970546   0.5334  0.007359699   0.4622
## 3   BF 21 -0.05018487 0.6531 -0.04453159   0.6406 -0.024435790   0.5790
## 4   BT 20 -0.16719714 0.9766 -0.17561797   0.9776 -0.034464525   0.5789
## 5   IF 19 -0.08476312 0.7312 -0.09241273   0.7518  0.052928550   0.2925
## 6   IT 11 -0.04301710 0.5891 -0.04022095   0.5735 -0.200568892   0.9408
##   p_geo_adj p_geo_pH_adj p_pH_geo_adj
## 1    0.9766       0.9776      0.75708
## 2    0.9766       0.9776      0.75708
## 3    0.9766       0.9776      0.75708
## 4    0.9766       0.9776      0.75708
## 5    0.9766       0.9776      0.75708
## 6    0.9766       0.9776      0.94080

Prediction of Tree growth (rcsBAI) - all OTUs

meta2 <- as(sample_data(vst_physeq), "data.frame")
keep2 <- complete.cases(meta2[, c("Age","pH","Habitat","Site", "rcsBAI_5y", "rcsBAI_10y","rcsBAI_15y","rcsBAI_30y")])
vst_mat <- as.data.frame(t(otu_table(vst_physeq)))
vst_mat_clean2 <- vst_mat[keep2, ]
meta_clean2 <- meta2[keep2, ]

meta_clean2$rcsBAI_5y_trans <- log(meta_clean2$rcsBAI_5y)
meta_clean2$rcsBAI_10y_trans <- log(meta_clean2$rcsBAI_10y)
meta_clean2$rcsBAI_15y_trans <- log(meta_clean2$rcsBAI_15y)
meta_clean2$rcsBAI_30y_trans <- log(meta_clean2$rcsBAI_30y)

# PCA on OTU table
pca <- prcomp(vst_mat_clean2, scale. = TRUE)
pcs_df <- as.data.frame(pca$x[, 1:5]) 

# Add BAI and environmental predictors
pcs_df$rcsBAI_5y_trans <- meta_clean2$rcsBAI_5y_trans
pcs_df$rcsBAI_10y_trans <- meta_clean2$rcsBAI_10y_trans
pcs_df$rcsBAI_15y_trans <- meta_clean2$rcsBAI_15y_trans
pcs_df$rcsBAI_30y_trans <- meta_clean2$rcsBAI_30y_trans
pcs_df$pH <- meta_clean2$pH
pcs_df$Habitat <- meta_clean2$Habitat
pcs_df$Site <- meta_clean2$Site

#  RDA using PCs + environmental variables
rda_comm <-  rda(pcs_df$rcsBAI_5y_trans ~ PC1 + PC2 + PC3 + PC4 + PC5 + Condition(Site + Habitat + pH),data = pcs_df)
anova(rda_comm, permutations = 9999)

## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 9999
## 
## Model: rda(formula = pcs_df$rcsBAI_5y_trans ~ PC1 + PC2 + PC3 + PC4 + PC5 + Condition(Site + Habitat + pH), data = pcs_df)
##          Df Variance      F Pr(>F)
## Model     5  0.00514 0.1529 0.9787
## Residual 96  0.64578

RsquareAdj(rda_comm)

## $r.squared
## [1] 0.006007931
## 
## $adj.r.squared
## [1] -0.03461367

rda_comm <- rda(pcs_df$rcsBAI_10y_trans ~ PC1 + PC2 + PC3 + PC4 + PC5 + Condition(Site + Habitat + pH),data = pcs_df)
anova(rda_comm, permutations = 9999)

## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 9999
## 
## Model: rda(formula = pcs_df$rcsBAI_10y_trans ~ PC1 + PC2 + PC3 + PC4 + PC5 + Condition(Site + Habitat + pH), data = pcs_df)
##          Df Variance      F Pr(>F)
## Model     5  0.00413 0.1241 0.9887
## Residual 96  0.63915

RsquareAdj(rda_comm)

## $r.squared
## [1] 0.004992691
## 
## $adj.r.squared
## [1] -0.03663512

rda_comm <- rda(pcs_df$rcsBAI_15y_trans ~ PC1 + PC2 + PC3 + PC4 + PC5 + Condition(Site + Habitat + pH),data = pcs_df)
anova(rda_comm, permutations = 9999)

## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 9999
## 
## Model: rda(formula = pcs_df$rcsBAI_15y_trans ~ PC1 + PC2 + PC3 + PC4 + PC5 + Condition(Site + Habitat + pH), data = pcs_df)
##          Df Variance     F Pr(>F)
## Model     5  0.00413 0.127 0.9871
## Residual 96  0.62398

RsquareAdj(rda_comm)

## $r.squared
## [1] 0.005170612
## 
## $adj.r.squared
## [1] -0.03696668

rda_comm <- rda(pcs_df$rcsBAI_30y_trans ~ PC1 + PC2 + PC3 + PC4 + PC5 + Condition(Site + Habitat + pH),data = pcs_df)
anova(rda_comm, permutations = 9999)

## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 9999
## 
## Model: rda(formula = pcs_df$rcsBAI_30y_trans ~ PC1 + PC2 + PC3 + PC4 + PC5 + Condition(Site + Habitat + pH), data = pcs_df)
##          Df Variance      F Pr(>F)
## Model     5  0.00754 0.2712 0.9279
## Residual 96  0.53365

RsquareAdj(rda_comm)

## $r.squared
## [1] 0.01094423
## 
## $adj.r.squared
## [1] -0.03057621

Prediction of Tree growth (rcsBAI) - only ECM

tax_data <- tax_tab
tax_data_ECM <- tax_data[tax_data$Guild == "Ectomycorrhizal", ]
vst_mat_clean2_ECM <- vst_mat_clean2[, colnames(vst_mat_clean2) %in% tax_data_ECM$OTU]

# PCA on OTU table
pca_ECM <- prcomp(vst_mat_clean2_ECM, scale. = TRUE)
pcs_df_ECM <- as.data.frame(pca_ECM$x[, 1:5])

# Add BAI and environmental predictors
pcs_df_ECM$rcsBAI_5y_trans <- meta_clean2$rcsBAI_5y_trans
pcs_df_ECM$rcsBAI_10y_trans <- meta_clean2$rcsBAI_10y_trans
pcs_df_ECM$rcsBAI_15y_trans <- meta_clean2$rcsBAI_15y_trans
pcs_df_ECM$rcsBAI_30y_trans <- meta_clean2$rcsBAI_30y_trans
pcs_df_ECM$pH <- meta_clean2$pH
pcs_df_ECM$Habitat <- meta_clean2$Habitat
pcs_df_ECM$Site <- meta_clean2$Site

#  RDA using PCs + environmental variables
rda_comm <- rda(pcs_df_ECM$rcsBAI_5y_trans ~ PC1 + PC2 + PC3 + PC4 + PC5 + Condition(Site + Habitat + pH),data = pcs_df_ECM)
anova(rda_comm, permutations = 9999)

## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 9999
## 
## Model: rda(formula = pcs_df_ECM$rcsBAI_5y_trans ~ PC1 + PC2 + PC3 + PC4 + PC5 + Condition(Site + Habitat + pH), data = pcs_df_ECM)
##          Df Variance      F Pr(>F)
## Model     5  0.02170 0.6621 0.6656
## Residual 96  0.62922

RsquareAdj(rda_comm)

## $r.squared
## [1] 0.02535657
## 
## $adj.r.squared
## [1] -0.01345109

rda_comm <- rda(pcs_df_ECM$rcsBAI_10y_trans ~ PC1 + PC2 + PC3 + PC4 + PC5 + Condition(Site + Habitat + pH),data = pcs_df_ECM)
anova(rda_comm, permutations = 9999)

## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 9999
## 
## Model: rda(formula = pcs_df_ECM$rcsBAI_10y_trans ~ PC1 + PC2 + PC3 + PC4 + PC5 + Condition(Site + Habitat + pH), data = pcs_df_ECM)
##          Df Variance      F Pr(>F)
## Model     5  0.01974 0.6078 0.6986
## Residual 96  0.62355

RsquareAdj(rda_comm)

## $r.squared
## [1] 0.02385566
## 
## $adj.r.squared
## [1] -0.01600374

rda_comm <- rda(pcs_df_ECM$rcsBAI_15y_trans ~ PC1 + PC2 + PC3 + PC4 + PC5 + Condition(Site + Habitat + pH),data = pcs_df_ECM)
anova(rda_comm, permutations = 9999)

## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 9999
## 
## Model: rda(formula = pcs_df_ECM$rcsBAI_15y_trans ~ PC1 + PC2 + PC3 + PC4 + PC5 + Condition(Site + Habitat + pH), data = pcs_df_ECM)
##          Df Variance      F Pr(>F)
## Model     5  0.01937 0.6108 0.6907
## Residual 96  0.60874

RsquareAdj(rda_comm)

## $r.squared
## [1] 0.02426958
## 
## $adj.r.squared
## [1] -0.01607719

rda_comm <- rda(pcs_df_ECM$rcsBAI_30y_trans ~ PC1 + PC2 + PC3 + PC4 + PC5 + Condition(Site + Habitat + pH),data = pcs_df_ECM)
anova(rda_comm, permutations = 9999)

## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 9999
## 
## Model: rda(formula = pcs_df_ECM$rcsBAI_30y_trans ~ PC1 + PC2 + PC3 + PC4 + PC5 + Condition(Site + Habitat + pH), data = pcs_df_ECM)
##          Df Variance      F Pr(>F)
## Model     5  0.01004 0.3628  0.877
## Residual 96  0.53115

RsquareAdj(rda_comm)

## $r.squared
## [1] 0.01457127
## 
## $adj.r.squared
## [1] -0.02660913

OTU abundancy per Taxa per Site - treeline vs forest (age-adjusted): Genus

# genus-level count matrix
rank <- "Genus"
tax_map <- as.data.frame(tax_tab[, c("OTU", rank)])
names(tax_map) <- c("OTU", "taxon")
tax_map <- tax_map[tax_map$taxon != "-" & !is.na(tax_map$taxon), ]
tax_map <- tax_map[tax_map$OTU %in% rownames(filt_count_tab_no_r), ]
taxon_counts <- rowsum(as.matrix(filt_count_tab_no_r[tax_map$OTU, ]),
                       group = tax_map$taxon)

# ---------- Alaska_Range
samples <- rownames(sample_info_no_R[sample_info_no_R$Site == "Alaska_Range", ])
samples <- samples[!is.na(sample_info_no_R[samples, "Age"])]
info <- sample_info_no_R[samples, ]
info$Habitat <- factor(info$Habitat, levels = c("forest", "treeline"))
info$Age_s   <- as.numeric(scale(info$Age))

mat    <- taxon_counts[, samples]
mat    <- mat[rowSums(mat > 0) >= 3, ]                          # taxa in at least 3 trees
mat_ab <- mat[rowSums(mat > 0) >= ceiling(0.30 * ncol(mat)), ]  # taxa in at least 30% of trees

dds <- DESeqDataSetFromMatrix(mat_ab, colData = info, design = ~ Age_s + Habitat)
dds <- estimateSizeFactors(dds, type = "poscounts")
res_hab <- results(DESeq(dds, test = "LRT", reduced = ~ Age_s,   fitType = "local", quiet = TRUE),
                   pAdjustMethod = "BH", independentFiltering = FALSE, cooksCutoff = FALSE)
res_age <- results(DESeq(dds, test = "LRT", reduced = ~ Habitat, fitType = "local", quiet = TRUE),
                   pAdjustMethod = "BH", independentFiltering = FALSE, cooksCutoff = FALSE)
da_ar <- data.frame(Site = "Alaska_Range", taxon = rownames(res_hab),
                  log2FC = res_hab$log2FoldChange, p_adj = res_hab$padj,
                  log2FC_age = res_age$log2FoldChange, p_adj_age = res_age$padj)

occ_ar <- data.frame(Site = "Alaska_Range", taxon = rownames(mat),
                   n_forest = NA, n_treeline = NA, p = NA)

# one Fisher test per taxon
for (i in 1:nrow(mat)) {   
  tb <- table(factor(mat[i, ] > 0, levels = c(FALSE, TRUE)), info$Habitat)
  occ_ar$n_forest[i]   <- tb["TRUE", "forest"]
  occ_ar$n_treeline[i] <- tb["TRUE", "treeline"]
  occ_ar$p[i]          <- fisher.test(tb)$p.value}

occ_ar$p_adj <- p.adjust(occ_ar$p, method = "BH")
occ_ar$p     <- NULL

# ---------- Brooks_Range
samples <- rownames(sample_info_no_R[sample_info_no_R$Site == "Brooks_Range", ])
samples <- samples[!is.na(sample_info_no_R[samples, "Age"])]
info <- sample_info_no_R[samples, ]
info$Habitat <- factor(info$Habitat, levels = c("forest", "treeline"))
info$Age_s   <- as.numeric(scale(info$Age))

mat    <- taxon_counts[, samples]
mat    <- mat[rowSums(mat > 0) >= 3, ]                          # taxa in at least 3 trees
mat_ab <- mat[rowSums(mat > 0) >= ceiling(0.30 * ncol(mat)), ]  # taxa in at least 30% of trees

dds <- DESeqDataSetFromMatrix(mat_ab, colData = info, design = ~ Age_s + Habitat)
dds <- estimateSizeFactors(dds, type = "poscounts")
res_hab <- results(DESeq(dds, test = "LRT", reduced = ~ Age_s,   fitType = "local", quiet = TRUE),
                   pAdjustMethod = "BH", independentFiltering = FALSE, cooksCutoff = FALSE)
res_age <- results(DESeq(dds, test = "LRT", reduced = ~ Habitat, fitType = "local", quiet = TRUE),
                   pAdjustMethod = "BH", independentFiltering = FALSE, cooksCutoff = FALSE)
da_br <- data.frame(Site = "Brooks_Range", taxon = rownames(res_hab),
                  log2FC = res_hab$log2FoldChange, p_adj = res_hab$padj,
                  log2FC_age = res_age$log2FoldChange, p_adj_age = res_age$padj)

occ_br <- data.frame(Site = "Brooks_Range", taxon = rownames(mat),
                   n_forest = NA, n_treeline = NA, p = NA)

# one Fisher test per taxon
for (i in 1:nrow(mat)) {   
  tb <- table(factor(mat[i, ] > 0, levels = c(FALSE, TRUE)), info$Habitat)
  occ_br$n_forest[i]   <- tb["TRUE", "forest"]
  occ_br$n_treeline[i] <- tb["TRUE", "treeline"]
  occ_br$p[i]          <- fisher.test(tb)$p.value}

occ_br$p_adj <- p.adjust(occ_br$p, method = "BH")
occ_br$p     <- NULL

# ---------- Interior
samples <- rownames(sample_info_no_R[sample_info_no_R$Site == "Interior", ])
samples <- samples[!is.na(sample_info_no_R[samples, "Age"])]
info <- sample_info_no_R[samples, ]
info$Habitat <- factor(info$Habitat, levels = c("forest", "treeline"))
info$Age_s   <- as.numeric(scale(info$Age))

mat    <- taxon_counts[, samples]
mat    <- mat[rowSums(mat > 0) >= 3, ]                          # taxa in at least 3 trees
mat_ab <- mat[rowSums(mat > 0) >= ceiling(0.30 * ncol(mat)), ]  # taxa in at least 30% of trees

dds <- DESeqDataSetFromMatrix(mat_ab, colData = info, design = ~ Age_s + Habitat)
dds <- estimateSizeFactors(dds, type = "poscounts")
res_hab <- results(DESeq(dds, test = "LRT", reduced = ~ Age_s,   fitType = "local", quiet = TRUE),
                   pAdjustMethod = "BH", independentFiltering = FALSE, cooksCutoff = FALSE)
res_age <- results(DESeq(dds, test = "LRT", reduced = ~ Habitat, fitType = "local", quiet = TRUE),
                   pAdjustMethod = "BH", independentFiltering = FALSE, cooksCutoff = FALSE)
da_in <- data.frame(Site = "Interior", taxon = rownames(res_hab),
                  log2FC = res_hab$log2FoldChange, p_adj = res_hab$padj,
                  log2FC_age = res_age$log2FoldChange, p_adj_age = res_age$padj)

occ_in <- data.frame(Site = "Interior", taxon = rownames(mat),
                   n_forest = NA, n_treeline = NA, p = NA)

# one Fisher test per taxon
for (i in 1:nrow(mat)) {   
  tb <- table(factor(mat[i, ] > 0, levels = c(FALSE, TRUE)), info$Habitat)
  occ_in$n_forest[i]   <- tb["TRUE", "forest"]
  occ_in$n_treeline[i] <- tb["TRUE", "treeline"]
  occ_in$p[i]          <- fisher.test(tb)$p.value}

occ_in$p_adj <- p.adjust(occ_in$p, method = "BH")
occ_in$p     <- NULL

da_all  <- rbind(da_ar, da_br, da_in)
occ_all <- rbind(occ_ar, occ_br, occ_in)

da_all[!is.na(da_all$p_adj)     & da_all$p_adj     < 0.05, ]

##             Site            taxon    log2FC        p_adj log2FC_age  p_adj_age
## 2   Alaska_Range        Beauveria  4.342565 4.150449e-07   4.342565 0.81476793
## 11  Alaska_Range   Cryptocoryneum -3.939678 1.865102e-03  -3.939678 0.81476793
## 14  Alaska_Range  Entomortierella -3.410494 1.957246e-03  -3.410494 0.81476793
## 23  Alaska_Range      Kockovaella -3.216385 8.858492e-03  -3.216385 0.81476793
## 42  Alaska_Range      Polyozellus  2.783391 2.110120e-04   2.783391 0.81476793
## 48  Alaska_Range       Sistotrema -7.888183 8.858492e-03  -7.888183 0.81476793
## 50  Alaska_Range         Tausonia  2.450374 1.705057e-02   2.450374 0.81476793
## 54  Alaska_Range      Trechispora  3.088146 1.865102e-03   3.088146 0.81476793
## 62  Brooks_Range         Capronia -2.787358 1.535374e-04  -2.787358 0.34258592
## 64  Brooks_Range  Chromolaenicola -4.518377 1.349909e-02  -4.518377 0.81676531
## 65  Brooks_Range Cladophialophora -1.222143 6.299402e-03  -1.222143 0.19516946
## 69  Brooks_Range   Cryptocoryneum -4.729571 1.512520e-02  -4.729571 0.58435677
## 71  Brooks_Range        Epicoccum  2.532034 7.039523e-03   2.532034 0.90050251
## 72  Brooks_Range        Exophiala  1.511187 1.538781e-02   1.511187 0.51429585
## 76  Brooks_Range      Hygrophorus -5.093006 1.559046e-02  -5.093006 0.12208945
## 77  Brooks_Range          Inocybe  2.065431 4.694431e-02   2.065431 0.64326468
## 84  Brooks_Range          Morinia  3.788824 2.913600e-03   3.788824 0.12208945
## 88  Brooks_Range        Mycenella  3.216261 4.739507e-02   3.216261 0.14633306
## 89  Brooks_Range         Niesslia -2.020244 1.559046e-02  -2.020244 0.51429585
## 100 Brooks_Range Ramoconidiophora  2.634735 2.798298e-03   2.634735 0.01220521
## 109 Brooks_Range       Tomentella  2.318341 2.798298e-03   2.318341 0.51429585
## 150     Interior      Penicillium  1.456904 3.402592e-02   1.456904 0.90990698

occ_all[occ_all$p_adj < 0.05, ]

##             Site          taxon n_forest n_treeline        p_adj
## 16  Alaska_Range Cryptocoryneum       17          3 0.0008353244
## 20  Alaska_Range      Epicoccum        4         14 0.0375532922
## 124 Brooks_Range       Lyomyces        0         11 0.0045718417

da_all[!is.na(da_all$p_adj_age) & da_all$p_adj_age < 0.05, ]   # taxa driven by age

##             Site            taxon    log2FC       p_adj log2FC_age   p_adj_age
## 100 Brooks_Range Ramoconidiophora 2.6347351 0.002798298  2.6347351 0.012205213
## 117     Interior       Calonarius 0.9460503 0.483862591  0.9460503 0.007543203

res_all <- merge(da_all, occ_all, by = c("Site", "taxon"),
                 all = TRUE, suffixes = c("_abund", "_occ"))
res_all$log2FC      <- signif(res_all$log2FC, 3)      # signif, not round: keeps small p-values
res_all$p_adj_abund <- signif(res_all$p_adj_abund, 3)
res_all$log2FC_age  <- signif(res_all$log2FC_age, 3)
res_all$p_adj_age   <- signif(res_all$p_adj_age, 3)
res_all$p_adj_occ   <- signif(res_all$p_adj_occ, 3)
write.csv(res_all, "Habitat_tests.csv", row.names = FALSE)

OTU abundancy per Taxa per Site - treeline vs forest (age-adjusted): Species

# genus-level count matrix
rank <- "Species...14"
tax_map <- as.data.frame(tax_tab[, c("OTU", rank)])
names(tax_map) <- c("OTU", "taxon")
tax_map <- tax_map[tax_map$taxon != "-" & !is.na(tax_map$taxon), ]
tax_map <- tax_map[tax_map$OTU %in% rownames(filt_count_tab_no_r), ]
taxon_counts <- rowsum(as.matrix(filt_count_tab_no_r[tax_map$OTU, ]),
                       group = tax_map$taxon)

# ---------- Alaska_Range
samples <- rownames(sample_info_no_R[sample_info_no_R$Site == "Alaska_Range", ])
samples <- samples[!is.na(sample_info_no_R[samples, "Age"])]
info <- sample_info_no_R[samples, ]
info$Habitat <- factor(info$Habitat, levels = c("forest", "treeline"))
info$Age_s   <- as.numeric(scale(info$Age))

mat    <- taxon_counts[, samples]
mat    <- mat[rowSums(mat > 0) >= 3, ]                          # taxa in at least 3 trees
mat_ab <- mat[rowSums(mat > 0) >= ceiling(0.30 * ncol(mat)), ]  # taxa in at least 30% of trees

dds <- DESeqDataSetFromMatrix(mat_ab, colData = info, design = ~ Age_s + Habitat)
dds <- estimateSizeFactors(dds, type = "poscounts")
res_hab <- results(DESeq(dds, test = "LRT", reduced = ~ Age_s,   fitType = "local", quiet = TRUE),
                   pAdjustMethod = "BH", independentFiltering = FALSE, cooksCutoff = FALSE)
res_age <- results(DESeq(dds, test = "LRT", reduced = ~ Habitat, fitType = "local", quiet = TRUE),
                   pAdjustMethod = "BH", independentFiltering = FALSE, cooksCutoff = FALSE)
da_ar <- data.frame(Site = "Alaska_Range", taxon = rownames(res_hab),
                  log2FC = res_hab$log2FoldChange, p_adj = res_hab$padj,
                  log2FC_age = res_age$log2FoldChange, p_adj_age = res_age$padj)

occ_ar <- data.frame(Site = "Alaska_Range", taxon = rownames(mat),
                   n_forest = NA, n_treeline = NA, p = NA)

# one Fisher test per taxon
for (i in 1:nrow(mat)) {   
  tb <- table(factor(mat[i, ] > 0, levels = c(FALSE, TRUE)), info$Habitat)
  occ_ar$n_forest[i]   <- tb["TRUE", "forest"]
  occ_ar$n_treeline[i] <- tb["TRUE", "treeline"]
  occ_ar$p[i]          <- fisher.test(tb)$p.value}

occ_ar$p_adj <- p.adjust(occ_ar$p, method = "BH")
occ_ar$p     <- NULL

# ---------- Brooks_Range
samples <- rownames(sample_info_no_R[sample_info_no_R$Site == "Brooks_Range", ])
samples <- samples[!is.na(sample_info_no_R[samples, "Age"])]
info <- sample_info_no_R[samples, ]
info$Habitat <- factor(info$Habitat, levels = c("forest", "treeline"))
info$Age_s   <- as.numeric(scale(info$Age))

mat    <- taxon_counts[, samples]
mat    <- mat[rowSums(mat > 0) >= 3, ]                          # taxa in at least 3 trees
mat_ab <- mat[rowSums(mat > 0) >= ceiling(0.30 * ncol(mat)), ]  # taxa in at least 30% of trees

dds <- DESeqDataSetFromMatrix(mat_ab, colData = info, design = ~ Age_s + Habitat)
dds <- estimateSizeFactors(dds, type = "poscounts")
res_hab <- results(DESeq(dds, test = "LRT", reduced = ~ Age_s,   fitType = "local", quiet = TRUE),
                   pAdjustMethod = "BH", independentFiltering = FALSE, cooksCutoff = FALSE)
res_age <- results(DESeq(dds, test = "LRT", reduced = ~ Habitat, fitType = "local", quiet = TRUE),
                   pAdjustMethod = "BH", independentFiltering = FALSE, cooksCutoff = FALSE)
da_br <- data.frame(Site = "Brooks_Range", taxon = rownames(res_hab),
                  log2FC = res_hab$log2FoldChange, p_adj = res_hab$padj,
                  log2FC_age = res_age$log2FoldChange, p_adj_age = res_age$padj)

occ_br <- data.frame(Site = "Brooks_Range", taxon = rownames(mat),
                   n_forest = NA, n_treeline = NA, p = NA)

# one Fisher test per taxon
for (i in 1:nrow(mat)) {  
  tb <- table(factor(mat[i, ] > 0, levels = c(FALSE, TRUE)), info$Habitat)
  occ_br$n_forest[i]   <- tb["TRUE", "forest"]
  occ_br$n_treeline[i] <- tb["TRUE", "treeline"]
  occ_br$p[i]          <- fisher.test(tb)$p.value}

occ_br$p_adj <- p.adjust(occ_br$p, method = "BH")
occ_br$p     <- NULL

# ---------- Interior
samples <- rownames(sample_info_no_R[sample_info_no_R$Site == "Interior", ])
samples <- samples[!is.na(sample_info_no_R[samples, "Age"])]
info <- sample_info_no_R[samples, ]
info$Habitat <- factor(info$Habitat, levels = c("forest", "treeline"))
info$Age_s   <- as.numeric(scale(info$Age))

mat    <- taxon_counts[, samples]
mat    <- mat[rowSums(mat > 0) >= 3, ]                          # taxa in at least 3 trees
mat_ab <- mat[rowSums(mat > 0) >= ceiling(0.30 * ncol(mat)), ]  # taxa in at least 30% of trees

dds <- DESeqDataSetFromMatrix(mat_ab, colData = info, design = ~ Age_s + Habitat)
dds <- estimateSizeFactors(dds, type = "poscounts")
res_hab <- results(DESeq(dds, test = "LRT", reduced = ~ Age_s,   fitType = "local", quiet = TRUE),
                   pAdjustMethod = "BH", independentFiltering = FALSE, cooksCutoff = FALSE)
res_age <- results(DESeq(dds, test = "LRT", reduced = ~ Habitat, fitType = "local", quiet = TRUE),
                   pAdjustMethod = "BH", independentFiltering = FALSE, cooksCutoff = FALSE)
da_in <- data.frame(Site = "Interior", taxon = rownames(res_hab),
                  log2FC = res_hab$log2FoldChange, p_adj = res_hab$padj,
                  log2FC_age = res_age$log2FoldChange, p_adj_age = res_age$padj)

occ_in <- data.frame(Site = "Interior", taxon = rownames(mat),
                   n_forest = NA, n_treeline = NA, p = NA)

# one Fisher test per taxon
for (i in 1:nrow(mat)) {   
  tb <- table(factor(mat[i, ] > 0, levels = c(FALSE, TRUE)), info$Habitat)
  occ_in$n_forest[i]   <- tb["TRUE", "forest"]
  occ_in$n_treeline[i] <- tb["TRUE", "treeline"]
  occ_in$p[i]          <- fisher.test(tb)$p.value}

occ_in$p_adj <- p.adjust(occ_in$p, method = "BH")
occ_in$p     <- NULL

da_all  <- rbind(da_ar, da_br, da_in)
occ_all <- rbind(occ_ar, occ_br, occ_in)

da_all[!is.na(da_all$p_adj)     & da_all$p_adj     < 0.05, ]

##             Site                          taxon    log2FC        p_adj
## 1   Alaska_Range       Beauveria_pseudobassiana  3.649517 0.0002903636
## 16  Alaska_Range           Cryptocoryneum_rosae -4.087998 0.0011466793
## 18  Alaska_Range     Entomortierella_parvispora -3.798724 0.0005579081
## 25  Alaska_Range        Kockovaella_prillingeri -3.426342 0.0056070119
## 33  Alaska_Range     Mortierella_basiparvispora  1.880578 0.0110938764
## 47  Alaska_Range     Penicillium_penicillioides  2.373798 0.0112795229
## 50  Alaska_Range           Penicillium_zhuangii  2.952082 0.0003625886
## 53  Alaska_Range        Piloderma_sphaerosporum  7.451743 0.0028908894
## 55  Alaska_Range           Polyozellus_umbrinus  2.294560 0.0019527050
## 61  Alaska_Range        Trechispora_stevensonii  3.220696 0.0259713585
## 69  Brooks_Range            Capronia_rubiginosa -3.127476 0.0002653942
## 71  Brooks_Range Chromolaenicola_chiangraiensis -4.670736 0.0072783114
## 72  Brooks_Range      Cladophialophora_humicola -2.749632 0.0253653950
## 80  Brooks_Range            Epicoccum_pimprinum  2.537788 0.0072783114
## 87  Brooks_Range                Inocybe_laurina  5.513002 0.0448568380
## 94  Brooks_Range            Morinia_trachycarpi  3.828502 0.0027966485
## 98  Brooks_Range                  Mycena_metata -7.631310 0.0014818585
## 114 Brooks_Range    Ramoconidiophora_euphorbiae  2.584934 0.0027966485
## 116 Brooks_Range               Sebacina_epigaea  5.458064 0.0485076829
## 123 Brooks_Range          Umbelopsis_ramanniana -3.494340 0.0052055625
## 152     Interior             Inocybe_flocculosa 14.618793 0.0119182030
## 169     Interior         Penicillium_aquadulcis  2.021580 0.0012831434
## 173     Interior    Penicillium_monsserratidens  4.179226 0.0011387603
## 174     Interior        Penicillium_raistrickii  2.450242 0.0150578861
##     log2FC_age    p_adj_age
## 1     3.649517 0.7246464746
## 16   -4.087998 0.6182380084
## 18   -3.798724 0.7145823689
## 25   -3.426342 0.7246464746
## 33    1.880578 0.5187994750
## 47    2.373798 0.5187994750
## 50    2.952082 0.5187994750
## 53    7.451743 0.5187994750
## 55    2.294560 0.7246464746
## 61    3.220696 0.5197813882
## 69   -3.127476 0.2000028228
## 71   -4.670736 1.0000000000
## 72   -2.749632 0.9995696864
## 80    2.537788 0.9995696864
## 87    5.513002 0.8365682379
## 94    3.828502 0.0757154193
## 98   -7.631310 0.5848784580
## 114   2.584934 0.0007831878
## 116   5.458064 0.8365682379
## 123  -3.494340 0.9995696864
## 152  14.618793 0.0925341429
## 169   2.021580 0.4796850033
## 173   4.179226 0.7883266662
## 174   2.450242 0.8738475504

occ_all[occ_all$p_adj < 0.05, ]

##             Site                   taxon n_forest n_treeline       p_adj
## 49  Alaska_Range    Cryptocoryneum_rosae       17          3 0.001388853
## 52  Alaska_Range     Epicoccum_pimprinum        4         14 0.041625336
## 108 Alaska_Range Piloderma_sphaerosporum        3         14 0.016130321
## 228 Brooks_Range        Lyomyces_sambuci        0         11 0.008665235

da_all[!is.na(da_all$p_adj_age) & da_all$p_adj_age < 0.05, ]   # taxa driven by age

##             Site                       taxon   log2FC       p_adj log2FC_age
## 114 Brooks_Range Ramoconidiophora_euphorbiae 2.584934 0.002796649   2.584934
## 117 Brooks_Range         Sebacina_incrustans 4.772042 0.138180141   4.772042
##        p_adj_age
## 114 0.0007831878
## 117 0.0240692922

res_all <- merge(da_all, occ_all, by = c("Site", "taxon"),
                 all = TRUE, suffixes = c("_abund", "_occ"))
res_all$log2FC      <- signif(res_all$log2FC, 3)      # signif, not round: keeps small p-values
res_all$p_adj_abund <- signif(res_all$p_adj_abund, 3)
res_all$log2FC_age  <- signif(res_all$log2FC_age, 3)
res_all$p_adj_age   <- signif(res_all$p_adj_age, 3)
res_all$p_adj_occ   <- signif(res_all$p_adj_occ, 3)
write.csv(res_all, "Habitat_tests_species.csv", row.names = FALSE)

OTU abundancy per Taxa per Site

tax_data <-tax_tab
tax_data <- tax_data[tax_data$OTU %in% rownames(filt_count_tab_no_r), ] # I used row counts
counts_tree <- as.data.frame(filt_count_tab_no_r)

counts_tree$OTU<-rownames(counts_tree)
rownames(counts_tree) <- NULL

# Combine tables sorting by OTU
Tax_otu_data <- left_join(tax_data, counts_tree,  by = "OTU" ) 
Sites <- unique(sample_info_no_R$Site)

set.seed(50)

samples <- rownames(sample_info_no_R[sample_info_no_R$Site == "Alaska_Range", ])
obj_ar <- parse_tax_data(Tax_otu_data, class_cols = "taxonomy", class_sep = ";")

## New names:
## • `Species...4` -> `Species...5`
## • `Species...14` -> `Species...15`

obj_ar <- filter_taxa(obj_ar, taxon_names != "-")
obj_ar$data$tax_abund <- calc_taxon_abund(obj_ar, data = "tax_data", cols = samples)

## Summing per-taxon counts from 40 columns for 437 taxa

obj_ar$data$tax_abund$total_counts <- rowSums(obj_ar$data$tax_abund[, samples])

samples <- rownames(sample_info_no_R[sample_info_no_R$Site == "Brooks_Range", ])
obj_br <- parse_tax_data(Tax_otu_data, class_cols = "taxonomy", class_sep = ";")

## New names:
## • `Species...4` -> `Species...5`
## • `Species...14` -> `Species...15`

obj_br <- filter_taxa(obj_br, taxon_names != "-")
obj_br$data$tax_abund <- calc_taxon_abund(obj_br, data = "tax_data", cols = samples)

## Summing per-taxon counts from 41 columns for 437 taxa

obj_br$data$tax_abund$total_counts <- rowSums(obj_br$data$tax_abund[, samples])

samples <- rownames(sample_info_no_R[sample_info_no_R$Site == "Interior", ])
obj_in <- parse_tax_data(Tax_otu_data, class_cols = "taxonomy", class_sep = ";")

## New names:
## • `Species...4` -> `Species...5`
## • `Species...14` -> `Species...15`

obj_in <- filter_taxa(obj_in, taxon_names != "-")
obj_in$data$tax_abund <- calc_taxon_abund(obj_in, data = "tax_data", cols = samples)

## Summing per-taxon counts from 30 columns for 437 taxa

obj_in$data$tax_abund$total_counts <- rowSums(obj_in$data$tax_abund[, samples])

results <- list(Alaska_Range = obj_ar, Brooks_Range = obj_br, Interior = obj_in)

obj_loc <- results[["Alaska_Range"]]
obj_loc_D <- filter_taxa(obj_loc, n_obs > 1)
set.seed(50)

hD<-heat_tree(obj_loc_D, node_label = gsub(pattern = "\\[|\\]", replacement = "", taxon_names),
            node_size = log(total_counts),
            node_size_range = c(0.0005, 0.025), 
            node_color = log(total_counts), # one taxa with 0 counts
            node_color_range = c("white","#33BBEE","#EE7733"),
            node_label_size_range = c(0.01, 0.025), 
            node_color_axis_label = "Abundance log(count number)",
            layout = "davidson-harel", initial_layout = "reingold-tilford")

## Warning: Infinite values found for "node_size". These will be graphed in the same way as the largest (Inf) or smallest (-Inf) real number supplied.

## Warning: Infinite values found for "edge_size". These will be graphed in the same way as the largest (Inf) or smallest (-Inf) real number supplied.

## Warning: Infinite values found for "node_label_size". These will be graphed in the same way as the largest (Inf) or smallest (-Inf) real number supplied.

## Warning: Infinite values found for "edge_label_size". These will be graphed in the same way as the largest (Inf) or smallest (-Inf) real number supplied.

plot(hD)

obj_loc <- results[["Brooks_Range"]]
obj_loc_B <- filter_taxa(obj_loc, n_obs > 1)

set.seed(50)

hB<-heat_tree(obj_loc_B, node_label = gsub(pattern = "\\[|\\]", replacement = "", taxon_names),
            node_size = log(total_counts),
            node_size_range = c(0.0005, 0.025), 
            node_color = log(total_counts), # no zeros
            node_color_range = c("white","#33BBEE","#EE7733"),
            node_label_size_range = c(0.01, 0.025), 
            node_color_axis_label = "Abundance log(count number)",
            layout = "davidson-harel", initial_layout = "reingold-tilford")
plot(hB)

obj_loc <- results[["Interior"]]
obj_loc_I <- filter_taxa(obj_loc, n_obs > 1)

set.seed(50)

hI<-heat_tree(obj_loc_I, node_label = gsub(pattern = "\\[|\\]", replacement = "", taxon_names),
            node_size = log(total_counts),
            node_size_range = c(0.0005, 0.025), 
            node_color = log(total_counts), #three taxa with 0 counts
            node_color_range = c("white","#33BBEE","#EE7733"),
            node_label_size_range = c(0.01, 0.025), 
            node_color_axis_label = "Abundance log(count number)",
            layout = "davidson-harel", initial_layout = "reingold-tilford")

## Warning: Infinite values found for "node_size". These will be graphed in the same way as the largest (Inf) or smallest (-Inf) real number supplied.

## Warning: Infinite values found for "edge_size". These will be graphed in the same way as the largest (Inf) or smallest (-Inf) real number supplied.

## Warning: Infinite values found for "node_label_size". These will be graphed in the same way as the largest (Inf) or smallest (-Inf) real number supplied.

## Warning: Infinite values found for "edge_label_size". These will be graphed in the same way as the largest (Inf) or smallest (-Inf) real number supplied.

plot(hI)

Species diversty (All OTUs): treeline vs forest

#Hill number = effective number of species.  Hill numbers determines the measures’ sensitivity to species relative abundances.  Hill numbers include the three most widely used species diversity measures as special cases: species richness (q = 0), Shannon diversity (q = 1, as the effective number of common species in the assemblage) and Simpson diversity (q = 2, effective number of dominant species in the assemblage)

#I used row counts of reads after removing of OTUs found in NTCs and technical replicates

count_tab <-filt_count_tab_no_r
count_tab <- as.data.frame(count_tab > 0) # to make presence/absence matrix

fD_info_rows <- sample_info_no_R[sample_info_no_R$Habitat == "forest" & sample_info_no_R$Site == "Alaska_Range",]
fD_counts_col <- count_tab[colnames(count_tab) %in% rownames(fD_info_rows) ]
tD_info_rows <- sample_info_no_R[sample_info_no_R$Habitat == "treeline" & sample_info_no_R$Site == "Alaska_Range", ]
tD_counts_col <- count_tab[colnames(count_tab) %in% rownames(tD_info_rows) ]

fI_info_rows <- sample_info_no_R[sample_info_no_R$Habitat == "forest" & sample_info_no_R$Site == "Interior",]
fI_counts_col <- count_tab[colnames(count_tab) %in% rownames(fI_info_rows) ]
tI_info_rows <- sample_info_no_R[sample_info_no_R$Habitat == "treeline" & sample_info_no_R$Site == "Interior", ]
tI_counts_col <- count_tab[colnames(count_tab) %in% rownames(tI_info_rows) ]

fB_info_rows <- sample_info_no_R[sample_info_no_R$Habitat == "forest" & sample_info_no_R$Site == "Brooks_Range",]
fB_counts_col <- count_tab[colnames(count_tab) %in% rownames(fB_info_rows) ]
tB_info_rows <- sample_info_no_R[sample_info_no_R$Habitat == "treeline" & sample_info_no_R$Site == "Brooks_Range", ]
tB_counts_col <- count_tab[colnames(count_tab) %in% rownames(tB_info_rows) ]


# create a dataframe with incidence data summed by rows within plots
inc_bytype<- cbind(rowSums(fD_counts_col), rowSums(tD_counts_col), rowSums(fI_counts_col), rowSums(tI_counts_col), rowSums(fB_counts_col), rowSums(tB_counts_col))

colnames(inc_bytype) <- c("A_Forest", "A_Treeline", "I_Forest", "I_Treeline","B_Forest", "B_Treeline")
inc_bytype<-as.data.frame(inc_bytype)
n_samples <-c(ncol(fD_counts_col), ncol(tD_counts_col), ncol(fI_counts_col), ncol(tI_counts_col), ncol(fB_counts_col), ncol(tB_counts_col))
rownames(inc_bytype)<-tax_data$OTU
inc_bytype <-rbind(n_samples, inc_bytype)

# Get the diversity estimates
out <- iNEXT(inc_bytype, q=c(0,1,2), datatype ="incidence_freq", se = TRUE)

# Plot the diversity estimates
# Sample-size-based rarefaction and extrapolation (R/E) curve (type=1). This curve plots diversity estimates with 95% confidence intervals (if se=TRUE) as a function of sample size up to double the reference sample size, by default, or a user‐specified endpoint:

ggiNEXT(out, type=1, facet.var="Order.q", se = TRUE, color.var="Assemblage")

ggiNEXT(out, type=2, facet.var="None", color.var="Assemblage", se = TRUE)

Species diversty (ECM OTUs): treeline vs forest

#I used row counts of reads after removing of OTUs found in NTCs and technical replicates
count_tab <-filt_count_tab_no_r

# filter OTU from tax_data:
tax_data     <- tax_tab[tax_tab$OTU %in% rownames(count_tab), ]
tax_data_ECM <- tax_data[tax_data$Guild %in% "Ectomycorrhizal", ]

count_tab <- count_tab[rownames(count_tab) %in% tax_data_ECM$OTU, ]
count_tab <- as.data.frame(count_tab > 0)   # presence/absence matrix

fD_info_rows <- sample_info_no_R[sample_info_no_R$Habitat == "forest" & sample_info_no_R$Site == "Alaska_Range",]
fD_counts_col <- count_tab[colnames(count_tab) %in% rownames(fD_info_rows) ]
tD_info_rows <- sample_info_no_R[sample_info_no_R$Habitat == "treeline" & sample_info_no_R$Site == "Alaska_Range", ]
tD_counts_col <- count_tab[colnames(count_tab) %in% rownames(tD_info_rows) ]

fI_info_rows <- sample_info_no_R[sample_info_no_R$Habitat == "forest" & sample_info_no_R$Site == "Interior",]
fI_counts_col <- count_tab[colnames(count_tab) %in% rownames(fI_info_rows) ]
tI_info_rows <- sample_info_no_R[sample_info_no_R$Habitat == "treeline" & sample_info_no_R$Site == "Interior", ]
tI_counts_col <- count_tab[colnames(count_tab) %in% rownames(tI_info_rows) ]

fB_info_rows <- sample_info_no_R[sample_info_no_R$Habitat == "forest" & sample_info_no_R$Site == "Brooks_Range",]
fB_counts_col <- count_tab[colnames(count_tab) %in% rownames(fB_info_rows) ]
tB_info_rows <- sample_info_no_R[sample_info_no_R$Habitat == "treeline" & sample_info_no_R$Site == "Brooks_Range", ]
tB_counts_col <- count_tab[colnames(count_tab) %in% rownames(tB_info_rows) ]

# create a dataframe with incidence data summed by rows within plots
inc_bytype<- cbind(rowSums(fD_counts_col), rowSums(tD_counts_col), rowSums(fI_counts_col), rowSums(tI_counts_col), rowSums(fB_counts_col), rowSums(tB_counts_col))

colnames(inc_bytype) <- c("A_Forest", "A_Treeline", "I_Forest", "I_Treeline","B_Forest", "B_Treeline")
inc_bytype<-as.data.frame(inc_bytype)
n_samples <-c(ncol(fD_counts_col), ncol(tD_counts_col), ncol(fI_counts_col), ncol(tI_counts_col), ncol(fB_counts_col), ncol(tB_counts_col))
inc_bytype <-rbind(n_samples, inc_bytype)

# Get the diversity estimates
out <- iNEXT(inc_bytype, q=c(0,1,2), datatype ="incidence_freq", se = TRUE)

# Plot the diversity estimates
# Sample-size-based rarefaction and extrapolation (R/E) curve (type=1). This curve plots diversity estimates with 95% confidence intervals (if se=TRUE) as a function of sample size up to double the reference sample size, by default, or a user‐specified endpoint:

ggiNEXT(out, type=1, facet.var="Order.q", se = TRUE, color.var="Assemblage")

Metabarcoding

Kristina Kuprina

2026-09-30

Treatment of NTC

Normalizing for sampling depth (variance stabilizing transformation)

PCoA (Plot)

Without technical replicates

Normalizing for sampling depth (variance stabilizing transformation)

PCoA (Plot, All OTUs)

PCoA (Plot, only ECM)

Relative Abundancy

Relative abundane - Genus

Prediction of the RAF composition (all OTUs): Variance partitioning

Redundance analyses - All OTUs

Prediction of the RAF composition (only ECM OTUs): Variance partitioning

Redundance analyses - Only ECM

Mantel test (all OTUs)

Mantel test (ECM OTUs)

Prediction of Tree growth (rcsBAI) - all OTUs

Prediction of Tree growth (rcsBAI) - only ECM

OTU abundancy per Taxa per Site - treeline vs forest (age-adjusted): Genus

OTU abundancy per Taxa per Site - treeline vs forest (age-adjusted): Species

OTU abundancy per Taxa per Site

Species diversty (All OTUs): treeline vs forest

Species diversty (ECM OTUs): treeline vs forest