R.Version()$version.string; R.Version()$platform
## [1] "R version 3.5.1 (2018-07-02)"
## [1] "x86_64-apple-darwin15.6.0"
library(picante)
packageVersion("picante")
## [1] '1.7'
library(phytools)
packageVersion("phytools")
## [1] '0.6.44'
library(PVR)
packageVersion("PVR")
## [1] '0.3'
library(adephylo)
packageVersion("adephylo")
## [1] '1.1.11'
library(phylobase)
packageVersion("phylobase")
## [1] '0.8.4'
library(geiger)
packageVersion("geiger")
## [1] '2.0.6'

Loading and processing the phylogeny

listasp<-read.table("prune_bocaina.txt")
PyronTree<-read.tree("amph_shl_new_Posterior_7238.1.tre")
filo<-prune.sample(t(listasp), PyronTree)
plot(filo, cex = 0.8);axisPhylo()

plotTree(filo)

Loading and processing trait data

trait<-read.csv2("canto.csv", h=TRUE, dec = ".")
habitat<-trait[,c(1,2)]
head(habitat)
##                           habitat microhabitat
## Bokermannohyla_ahenea           3            2
## Bokermannohyla_circumdata       3            1
## Aplastodiscus_perviridis        3            2
## Hypsiboas_bandeirantes          3            2
## Scinax_hayii                    3            2
## Scinax_duartei                  3            2
trait<-trait[,-c(1,2,4,8)]#excluding variables following the M&M of the main text
head(trait)
##                           tipos_notas    F_max     F_min    F_dom
## Bokermannohyla_ahenea               2 1717.933  536.0533 1240.313
## Bokermannohyla_circumdata           1 1597.915  629.7100 1044.375
## Aplastodiscus_perviridis            1 2371.445 1742.8327 2054.665
## Hypsiboas_bandeirantes              1 5658.853 4167.4333 4857.893
## Scinax_hayii                        1 3215.920  492.2267  861.300
## Scinax_duartei                      1 3554.983 1668.8067 3123.750
##                           duracao_canto  n_pulsos seas_repro
## Bokermannohyla_ahenea         1.2794000 66.366667          1
## Bokermannohyla_circumdata     0.1282000  3.900000          1
## Aplastodiscus_perviridis      0.1346727  1.000000          1
## Hypsiboas_bandeirantes        0.1236667 20.800000          2
## Scinax_hayii                  0.2608667 10.066667          1
## Scinax_duartei                0.3027667  8.166667          2
notes=as.matrix(trait)[,1]
names(notes)=rownames(trait)
F_max=as.matrix(trait)[,2]
names(F_max)=rownames(trait)
F_min=as.matrix(trait)[,3]
names(F_min)=rownames(trait)
F_dom=as.matrix(trait)[,4]
names(F_dom)=rownames(trait)
durcall=as.matrix(trait)[,5]
names(durcall)=rownames(trait)
pulses=as.matrix(trait)[,6]
names(pulses)=rownames(trait)
saz_rep=as.matrix(trait)[,7]
names(saz_rep)=rownames(trait)

Loading spatial and temporal distribution

spatial<-read.csv2("coespacial.csv", h=TRUE)
head(spatial)
##                           PT1 PT2 PT3 AP1 AP2 BT1 BP2 BP4 BP5 BP8 BP9 PP1
## Aplastodiscus_perviridis    0   0   0   0   1   0   9   3   2   4   0   0
## Bokermannohyla_ahenea       0   0   0   0   0   0   2   0   0   0   0   0
## Bokermannohyla_circumdata   0   0   0   0   0   0   0   0   0   0   2   0
## Dendropsophus_elegans       0   0   0   0   0   0   0   0   0   0   0   0
## Dendropsophus_microps      23  65  45  25  45  13   0   0   0   0   0   5
## Dendropsophus_minutus       7   0   0  45  12   7   0   0   0   0   0   4
##                           PP2 PP3 PP4
## Aplastodiscus_perviridis    0   0   6
## Bokermannohyla_ahenea       0   0   0
## Bokermannohyla_circumdata   0   0   2
## Dendropsophus_elegans       0   2   0
## Dendropsophus_microps       7   8   1
## Dendropsophus_minutus       0   7   5
co.temporal<-read.csv2("cotemporal.csv",h=TRUE)
head(co.temporal)
##                           Jul Aug Sep Oct Nov Dec Jan Feb Mar Apr May Jun
## Aplastodiscus_perviridis    0   0   0   5  15  17  18   9   0   0   0   0
## Bokermannohyla_ahenea       0   0   0   0   0   0   2   0   0   0   0   0
## Bokermannohyla_circumdata   0   0   0   0   0   2   2   0   0   0   0   0
## Dendropsophus_elegans       0   0   0   0   0   0   1   2   0   0   0   0
## Dendropsophus_microps       0   0  14  76  13  20 137  18   1   0   0   0
## Dendropsophus_minutus       0   9   6  42  11  32  71  34  14   0   0   0

Running a phylogentic PCA to summarize the advertisement call data

obj<-phylo4d(filo, trait[,-7])
table.phylo4d(obj, box = FALSE, cex.label = 0.7)

phylo.pca<-ppca(obj, method = "patristic",scannf=FALSE, nfposi=2)#standardizing traits by default to 0 mean and unit variance
scatter(phylo.pca, useLag = TRUE)

plot(phylo.pca, useLag = TRUE)

phylo.pca$c1#loadings. PC1 => espectral eigenvector | PC2 => temporal eigenvector
##                       PA1           PA2
## tipos_notas   -0.33810637 -0.0672892960
## F_max          0.07823164 -0.5908032733
## F_min         -0.61114391 -0.3656785193
## F_dom         -0.30422851  0.7145035327
## duracao_canto  0.52403091 -0.0001389401
## n_pulsos      -0.37269778 -0.0467712079
vectors<-phylo.pca$li#eigenvectors of PCA

Naming the vectors

espectral=as.matrix(vectors)[,1]
names(pulses)=rownames(vectors)
temporal=as.matrix(vectors)[,2]
names(temporal)=rownames(vectors)
new.data2<-treedata(filo, habitat)
new.data<-treedata(filo, vectors, sort = TRUE)

Plotting trait space

phenogram(filo, espectral, fsize = 0.8, ftype = "i", ylab= "Spectral eigenvector of the call")

phenogram(filo, temporal, fsize = 0.8, ftype = "i", ylab= "Temporal eigenvector of the call")

Running PVR

Preparing data

x<-PVRdecomp(filo)
set.seed(1001)#to allow reproducilibity of random numbers

building the PSR curve (with both PCs)

psrcurve<-PSR(x, trait = vectors, null.model = TRUE, Brownian.model = TRUE, times = 999)
PSRplot(psrcurve, info = "both")#Top plot refers to Null model; bottom plot refers to neutral model (Brownian Motion)

psrcurve.esp<-PSR(x, trait = espectral, null.model = TRUE, Brownian.model = TRUE, times = 999)
PSRplot(psrcurve.esp, info = "both")

psrcurve.tem<-PSR(x, trait = temporal, null.model = TRUE, Brownian.model = TRUE, times = 999)
PSRplot(psrcurve.tem, info = "both")

Running PVR for spectral parameters

anal_pvr<-PVR(psrcurve, trait = espectral, envVar = new.data2$data)
anal_pvrphy<-PVR(psrcurve, trait = espectral)
round(anal_pvrphy@PVR$R2, 3)#phylogenetic signal
## [1] 0.21
VarPartplot(anal_pvr)#plot of variation partitioning

## [a] - Variation explained by environmental variables
## [b] - Shared variation between environmental variables and phylogeny(PVR)
## [c] - Variation explained by phylogeny(PVR)
## [d]     - Unexplained variation
str(anal_pvr@VarPart)#R2 of each component of variation partitioning
## List of 4
##  $ a: num 0.0845
##  $ b: num 0.151
##  $ c: num 0.0587
##  $ d: num 0.706

Running PVR for temporal parameters

anal_pvr2<-PVR(psrcurve, trait = temporal, envVar = new.data2$data)
anal_pvr2phy<-PVR(psrcurve, trait = temporal)
round(anal_pvr2phy@PVR$R2, 3)#phylogenetic signal
## [1] 0.124
VarPartplot(anal_pvr2)#plot of variation partitioning

## [a] - Variation explained by environmental variables
## [b] - Shared variation between environmental variables and phylogeny(PVR)
## [c] - Variation explained by phylogeny(PVR)
## [d]     - Unexplained variation
str(anal_pvr2@VarPart)#R2 of each component of variation partitioning
## List of 4
##  $ a: num 0.171
##  $ b: num -0.00512
##  $ c: num 0.129
##  $ d: num 0.705

Analysis of trait diversity and community assembly

Correlation between species co-occurrence and trait distance

Building the species co-occurrence matrix based on Schoener’s index and C-Score. The higher the Schoener’s index, the higher the co-occurrence.

cooc.s<-species.dist(t(spatial), "cij")#Spatial co-occurrence
cooc.t<-species.dist(t(co.temporal), "cij")#Temporal co-occurrence
C_score_spa<-designdist(spatial, "(A-J)*(B-J)", "binary")
C_score_temp<-designdist(spatial, "(A-J)*(B-J)", "binary")

Calculating the trait distance

tr.dist<-vegdist(trait, "gower")#All call parameters
tr.dist.esp<-vegdist(trait[,1:4], "gower")#spectral parameters
tr.dist.temp<-vegdist(trait[,5:6], "gower")#temporal parameters

Correlating the trait distance and co-occurrence matrix

mantel(cooc.s, tr.dist)#correlation between spatial co-occurrence and all call attributes
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = cooc.s, ydis = tr.dist) 
## 
## Mantel statistic r: 0.04322 
##       Significance: 0.339 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.106 0.127 0.147 0.171 
## Permutation: free
## Number of permutations: 999
plot(cooc.s, tr.dist)
lines(lowess(cooc.s, tr.dist), col="black", lty=2)

mantel(cooc.t, tr.dist)#correlation between temporal co-occurrence and all call attributes
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = cooc.t, ydis = tr.dist) 
## 
## Mantel statistic r: 0.1777 
##       Significance: 0.073 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.158 0.197 0.227 0.264 
## Permutation: free
## Number of permutations: 999
plot(cooc.t, tr.dist)
lines(lowess(cooc.t, tr.dist), col="black", lty=2)

mantel(cooc.s, tr.dist.esp)#correlation between spatial co-occurrence and spectral traits
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = cooc.s, ydis = tr.dist.esp) 
## 
## Mantel statistic r: -0.01684 
##       Significance: 0.573 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.104 0.128 0.149 0.184 
## Permutation: free
## Number of permutations: 999
plot(cooc.s, tr.dist.esp)
lines(lowess(cooc.s, tr.dist.esp), col="black", lty=2)

mantel(cooc.t, tr.dist.esp)#correlation between temporal co-occurrence and spectral attributes
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = cooc.t, ydis = tr.dist.esp) 
## 
## Mantel statistic r: 0.0559 
##       Significance: 0.35 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.166 0.202 0.224 0.265 
## Permutation: free
## Number of permutations: 999
plot(cooc.t, tr.dist.esp)
lines(lowess(cooc.t, tr.dist.esp), col="black", lty=2)

mantel(cooc.s, tr.dist.temp)#correlation between spatial co-occurrence and temporal attributes
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = cooc.s, ydis = tr.dist.temp) 
## 
## Mantel statistic r: 0.1129 
##       Significance: 0.12 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.120 0.148 0.164 0.176 
## Permutation: free
## Number of permutations: 999
plot(cooc.s, tr.dist.temp)
lines(lowess(cooc.s, tr.dist.temp), col="black", lty=2)

mantel(cooc.t, tr.dist.temp)#correlation between temporal co-occurrence and temporal attributes
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = cooc.t, ydis = tr.dist.temp) 
## 
## Mantel statistic r: 0.2836 
##       Significance: 0.033 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.229 0.263 0.292 0.318 
## Permutation: free
## Number of permutations: 999
plot(cooc.t, tr.dist.temp)
lines(lowess(cooc.t, tr.dist.temp), col="black", lty=2)

As you can see, no correlation was significant, but the relationship between temporal co-occurrence and temporal parameters of the call.