####Load packages####
library(vegan)
library(readxl)
library(fossil) 
library(paleotree) 
library(divDyn)
library(chronosphere)
library(dendextend)
library(tidyverse)
library(dplyr)
library(purrr)
library(scales)
library(stats)
library(ggplot2)
library(plotrix) 
library(forcats)
library(cluster)
library(gridExtra)
library(extrafont)

# Set the working directory to where your data 
file_path <- "E:/file"
setwd(file_path)

# Read the dataset from an Excel file
data_file <- "TriaSCstg.xlsx"
sheet_name <- "Sheet1"
data <- read_excel(paste0(file_path, "/", data_file), sheet = sheet_name)

# Create a new data frame with the necessary columns and calculate new ages based on section
# Data restructuring
age <- data$`Age resolution`
stg <- data$Stage
abu <- data$Abundance
gen <- data$`Genus / subgenus or higher taxon`
spe <- data$Species
gensp <- paste(gen, spe)
sample <- data$Sample
section <- data$Section
lev <- data$`Meters in section`
mol <- data$`Mode of life`
depth <-  data$Depth
lit <- data$Lithology


sbiv <- data.frame(sample, section, lev, age, stg,gen, spe, gensp, abu, mol, depth, lit)
sbiv$n.age <- rep(NA, nrow(sbiv))

#### Setting up age models for geological sections ####
sects <- unique(section)
sects

###Liangshuichong
liang <- subset(sbiv, section=="Liangshuichong")

#add approximate time inferred from sample level 
oldest <- min(liang$lev)
youngest <- max(liang$lev)
x <- c(251.9, 248.1)
y <- c(oldest, youngest)
mod <- lm(x ~ y)
summary(mod)
l.ages <- mod$coefficients[1] + mod$coefficients[2] * liang$lev

sbiv$n.age[section == "Liangshuichong"] <- l.ages # Assign to n.age column

###Maotai
mao <- subset(sbiv, section=="Maotai")

oldest <- min(mao$lev)
youngest <- max(mao$lev)
x <- c(244.2, 242.3)
y <- c(oldest, youngest)
mod <- lm(x ~ y)
summary(mod)
l.ages <- mod$coefficients[1] + mod$coefficients[2]*mao$lev

sbiv$n.age[section == "Maotai"] <- l.ages

###Gaoqiao
gao <- subset(sbiv, section=="Gaoqiao")

#Filter data for Fassanian
gao <- subset(gao, age == "Fassanian")

oldest <- min(gao$lev)
youngest <- max(gao$lev)
x <- c(241.46, 241)
y <- c(oldest, youngest)
mod <- lm(x ~ y)
summary(mod)
l.ages_fas <- mod$coefficients[1] + mod$coefficients[2]*gao$lev

sbiv$n.age[section == "Gaoqiao"] <- l.ages_fas

###Hehuachi
hehua <- subset(sbiv, section=="Hehuachi")
oldest <- min(hehua$lev)
youngest <- max(hehua$lev)
x <- c(237, 233)
y <- c(oldest, youngest)
mod <- lm(x ~ y)
summary(mod)
l.ages <- mod$coefficients[1] + mod$coefficients[2]*hehua$lev

sbiv$n.age[section=="Hehuachi"] <- l.ages


####Statistical Analysis of South China ####

##### Alpha Diversity ####

##the stage level time scale

data(stages)
#finding which interval the 'early_interval' and 'late_interval' fields point to

data(keys)

n.sbiv <- sbiv
n.sbiv <- n.sbiv %>%
  mutate(age = gsub("Julian1|Julian2", "Julian", age))

stgMin <- categorize(n.sbiv[ ,"age"], keys$stgInt)
stgMin <- as.numeric(stgMin)

# empty container
n.sbiv$stg <- rep(NA, nrow(n.sbiv))
# Add "stg"
n.sbiv$stg <- stgMin

#load the substage data
substg <- read.csv(file = "bstages.csv", header = TRUE, sep = ",")

#add substage column to the data
n.sbiv$substg <- rep(NA, nrow(n.sbiv))
n.sbiv$substg <- rep(NA, nrow(n.sbiv))

n.sbiv <- merge(n.sbiv, substg[, c('substage', 'sstg')], 
                by.x = 'age', by.y = 'substage', all.x = TRUE) # Add substages

n.sbiv$substg <- n.sbiv$sstg



#### Samples ####
samples <- create.matrix(sbiv, tax.name = "gensp", locality="sample",
                         abund.col = "abu", abund=TRUE)


# limit analysis to some minimum value
y <- which(colSums(samples) > 30)
samples <- samples[,y]

# delete empty rows (could also be used to exclude rare species)
x <- which(rowSums(samples)> 0)
samples <- samples[x,]
r <- vegdist(samples, "bray")

# now turn the matrix around for Q-mode analyses
t.samples <- t(samples)

# compute a distance matrix

q <- vegdist(t.samples, "bray")
# perform a cluster analysis 
clu <- hclust(q, method="ward.D2")
plot(clu, hang = -1, main="Species")
p.table <- prop.table(samples, 2)

#Ecologic occurrences by sample
#This function is for the taxonomic names, may not fit ecology
e.samples <- create.matrix(sbiv, tax.name = "mol", locality="sample",
                           abund.col = "abu", abund=TRUE)

# limit analysis to some minimum value
y <- which(colSums(e.samples) > 30)
e.samples <- e.samples[,y]
# delete empty rows (could also be used to exclude rate species)
x <- which(rowSums(e.samples)>0)
e.samples <- e.samples[x,]
e.r <- vegdist(e.samples, "bray")

# now turn the matrix around for Q-mode analyses
t.e.samples <- t(e.samples)

# compute a distance matrix

e.q <- vegdist(t.e.samples, "bray")

# perform a cluster analysis 
e.clu <- hclust(e.q, method="ward.D2")
plot(e.clu, hang = -1, main="Ecology")

e.p.table <- prop.table(e.samples, 2)



##### Cluster Analysis of South China #####

# Taxonomic cluster #

###Colored Cluster of South China  ###
dend <- as.dendrogram(clu)

i = 0

colLab <<- function(n) {
  if (is.leaf(n)) {
    a = attributes(n)
    ligne = match(a$label, sbiv[, 1])
    age = sbiv[ligne, 4]
    
    col_age = switch(
      age,
      Griesbachian = "#E41A1C",
      Dienerian = "#377EB8",
      Smithian = "#4DAF4A",
      Spathian = "#00008B",
      Pelsonian = "#FFD700",
      Illyrian = "#FF7F00",
      Fassanian = "#984EA3",
      Julian1 = "#A65628",
      Julian2 = "#A65628",
      Tuvalian = "#F781BF",
      "black"  # Default color
    )
    
    #Set label color
    attr(n, "nodePar") <- c(
      a$nodePar,
      list(
        lab.col = col_age,
        lab.font = 1,
        pch = NA
        #,lab.cex = 0.8
      )
    )
    
    #set node point color
    # attr(n, "edgePar") <- list(col=col_age)
    
  }
  return(n)
}


######## Fig. 3 ########
#Taxonomy
dend<- as.dendrogram(clu)
dL <- dendrapply(dend, colLab)

pdf("Fig3A.pdf",  width = 10, height = 10, family = "ArialMT")

plot(dL, cex= 1.2, horiz = TRUE)

dev.off()

# Ecology Cluster
e.dend <- as.dendrogram(e.clu)
e.dL <- dendrapply(e.dend, colLab)

pdf("Fig3B.pdf",  width = 10, height = 10, family = "ArialMT")


plot(e.dL, cex= 1.2, horiz = TRUE)

legend(x="topleft", legend=c("Griesbachian","Dienerian","Smithian", "Pelsonian","Illyrian", "Fassanian","Julian"),
       col=c("#E41A1C", "#377EB8", "#4DAF4A", "#FFD700", "#FF7F00","#984EA3", "#A65628"), 
       cex=1, lty=1,lwd=2)

dev.off()


######### Fig. S1 ########

pdf("FigS1.pdf",  width = 5, height = 10, family = "ArialMT")
plot(dend,  horiz = TRUE)
abline(v = 1, col="#E41A1C", lwd=3, lty=2)
dev.off()


##### Two-way cluster of South China #####
######### Fig. S3 ########

# Taxonomy
pdf("FigS3A.pdf",  width = 10, height = 10, family = "ArialMT")
twoWayEcologyCluster(r, q, propAbund=p.table, clustMethod="ward.D2", trimChar = 16)
dev.off()


# Ecology
pdf("FigS3B.pdf",  width = 10, height = 10, family = "ArialMT")
twoWayEcologyCluster(e.r, e.q, propAbund=e.p.table, clustMethod="ward.D2")
dev.off()


##### NMDS of south China #####

### Taxonomic NMDS ###
m <- metaMDS(t.samples)

m$stress
stressplot(m)

# Hypotesis testing
#Group by substage
#Gri:1 Die:2 Smi:3 Spa:4 Pel:5 Ill:6 Fas:7 Jul:8 Tuv:9

groups <- factor(c(rep(7, 3), 8, rep(8, 7), rep(1, 4), rep(2,6),
                   3, 1, rep(5, 6),  6))

#### Table S1 ####
#ANOSIM
asim <- anosim(t.samples,groups)
asim
plot(asim)

### Eco NMDS ###
e.m <- metaMDS(t.e.samples)

e.m$stress
stressplot(e.m)

#ANOSIM
e.asim <- anosim(t.e.samples,groups)
e.asim
plot(e.asim)

#Group By Stage#
gri <- unique(sbiv$sample[sbiv$age=="Griesbachian"])
die <- unique(sbiv$sample[sbiv$age=="Dienerian"])
smi <- unique(sbiv$sample[sbiv$age=="Smithian"])
pel <- unique(sbiv$sample[sbiv$age=="Pelsonian"])
ill <- unique(sbiv$sample[sbiv$age=="Illyrian"])
lad <- unique(sbiv$sample[sbiv$age=="Fassanian"])
jul <- unique(sbiv$sample[sbiv$age=="Julian1"|sbiv$age=="Julian2"])

age2 <- rep("", nrow(m$points))
age2[row.names(m$points) %in% gri] <- "gri"
age2[row.names(m$points) %in% die] <- "die"
age2[row.names(m$points) %in% smi] <- "smi"
age2[row.names(m$points) %in% pel] <- "pel"
age2[row.names(m$points) %in% ill] <- "ill"
age2[row.names(m$points) %in% lad] <- "lad"
age2[row.names(m$points) %in% jul] <- "jul"

age2 <- factor(age2, levels = c("gri", "die", "smi","pel", "ill", "lad", "jul"))

#Taxonomy

######## Fig. 4 #######
pdf("Fig4A.pdf",  width = 10, height = 10, family = "ArialMT")

mdscol = c(rep("#984EA3", 3),  rep("#A65628", 8), rep("#E41A1C", 4), rep("#377EB8",6),
           "#458B00", "#E41A1C", rep("#E6AB02", 6),  rep("#FF7F00", 1))

#Plot blank figure
ordiplot(m, type="none")

#add species labels
orditorp(m,display="species",col= "grey60", air=0.01,cex=1)

#add sample labels
orditorp(m,display="sites",col= mdscol, air=0.01,cex=1.25)

points(m, display = "sites", col= mdscol, pch = 19, cex = 2)


ordihull(m, groups=age2, col=c("#E41A1C","#377EB8","#458B00",
                               "#E6AB02","#FF7F00", "#984EA3",
                               "#A65628","#F781BF"))

dev.off()

#Ecology
pdf("Fig4B.pdf",  width = 10, height = 10, family = "ArialMT")

ordiplot(e.m, type="none")

#add species labels
orditorp(e.m,display="species",col= "grey60", air=0.01,cex=1)

#add sample labels
orditorp(e.m,display="sites",col= mdscol, air=0.01,cex=1.25)

points(e.m, display = "sites", col= mdscol, pch = 19,cex = 2)

ordihull(e.m, groups=age2, col=c("#E41A1C","#377EB8","#458B00",
                                 "#E6AB02","#FF7F00", "#984EA3",
                                 "#A65628","#F781BF"))

legend(x="topright", legend=c("Griesbachian","Dienerian","Smithian", "Pelsonian","Illyrian", "Fassanian", "Julian"),
col=c("#E41A1C","#377EB8","#458B00", "#E6AB02","#FF7F00", "#984EA3", "#A65628"), 
cex=1, lty=1,lwd=2)

dev.off()


##### Taxonomy and Ecology Correlation of South China #####

#Taxonomic diversity index
s.div <- diversity(samples, MARGIN=2)#Shannon
s.S <- specnumber(samples, MARGIN=2)#Number of species
s.J <- s.div/log(s.S)#Evenness

#Ecologic diversity index
e.div <- diversity(e.samples, MARGIN=2)#Shannon
e.S <- specnumber(e.samples, MARGIN=2)#Number of species
e.J <- e.div/log(e.S)#Evenness

#Correlation Between Paired Samples
cor.div <- cor.test(s.div, e.div)$estimate
cor.div
cor.S <- cor.test(s.S, e.S)$estimate
cor.S
cor.J <- cor.test(s.J, e.J)$estimate
cor.J

#X-Y Plotting and Add Straight Lines to Plots
p1 <- plot(s.S, e.S)
mode1 <- lm(e.S~s.S)
summary(mode1)
abline(lm(e.S~s.S))
text(x = max(s.S) * 0.25, y = max(e.S) * 1, labels = paste("R =", round(cor.S, 2)))

p2 <- plot(s.div, e.div)
mode2 <- lm(e.div~s.div)
abline(lm(e.div~s.div))
text(x = max(s.div) * 0.2, y = max(e.div) * 1, labels = paste("R =", round(cor.div, 2)))

p3 <- plot(s.J, e.J)
mode3 <- lm(e.J~s.J)
abline(lm(e.J~s.J))
text(x = 0.25, y = 0.95, labels = paste("R =", round(cor.J, 2)))

####### Fig. 6 #######
#create pdf
pdf("Fig6.pdf",  width = 11, height = 4, family = "ArialMT")

#Divide into three columns
par(mfrow = c(1, 3))

# Create scatter plot for s.S vs. e.S
p1 <- ggplot(data = data.frame(s.S, e.S), aes(x = s.S, y = e.S)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +  # Add confidence interval
  labs(x = "Taxonomic Richness (s.S)", y = "Ecologic Richenss (e.S)") +
  theme_minimal() +
  geom_text(
    # x = max(s.S) * 0.2, y = max(e.S) * 1, 
    x = 5, y =8, 
    label = paste("R =", round(cor.S, 2)), hjust = 0, vjust = 1)+
  theme(panel.grid = element_blank(),  # Remove lattice-like background
        panel.border = element_rect(color = "black", fill = NA, size = 1))  # Add square border

# Create scatter plot for s.div vs. e.div
p2 <- ggplot(data = data.frame(s.div, e.div), aes(x = s.div, y = e.div)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +  # Add confidence interval
  labs(x = "Taxonomic Shannon Diversity (s.div)", y = "Ecologic Shannon Diversity (e.div)") +
  theme_minimal() +
  geom_text(
    #x = max(s.div) * 0.2, y = max(e.div) * 1,
    x = 0.25, y =1.65,  
    label = paste("R =", round(cor.div, 2)), hjust = 0, vjust = 1)+
  theme(panel.grid = element_blank(),  # Remove lattice-like background
        panel.border = element_rect(color = "black", fill = NA, size = 1))  # Add square border  # Remove lattice-like background

# Create scatter plot for s.J vs. e.J
p3 <- ggplot(data = data.frame(s.J, e.J), aes(x = s.J, y = e.J)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +  # Add confidence interval
  labs(x = "Taxonomic Evenness Index (s.J)", y = "Ecologic Evenness Index (e.J)") +
  theme_minimal() +
  geom_text(
    #x = max(s.J) * 0.2, y = max(e.J) * 1,
    x = 0.25, y =0.95,   
    label = paste("R =", round(cor.J, 2)), hjust = 0, vjust = 1)+
  theme(panel.grid = element_blank(),  # Remove lattice-like background
        panel.border = element_rect(color = "black", fill = NA, size = 1))  # Add square border

# Combine the plots (optional)
grid.arrange(p1, p2, p3, ncol = 3)

dev.off()

#### Table S2 ####
#Spearman rank-order correlations

src.S <- cor.test(s.S, e.S, method = "spearman", exact = FALSE)
src.S
src.div <- cor.test(s.div, e.div, method = "spearman")
src.div
src.J <- cor.test(s.J, e.J, method = "spearman")
src.J


#####Dynamic Diversity of South China#####
#the stage level time scale
data(stages)
#finding which interval the 'early_interval' and 'late_interval' fields point to
data(keys)

n.sbiv <- sbiv
n.sbiv <- n.sbiv %>%
  mutate(age = gsub("Julian1|Julian2", "Julian", age))

stgMin <- categorize(n.sbiv[ ,"age"], keys$stgInt)
stgMin <- as.numeric(stgMin)

# empty container
n.sbiv$stg <- rep(NA, nrow(n.sbiv))
# Add "stg"
n.sbiv$stg <- stgMin

#load the substage data
substg <- read.csv(file = "bstages.csv", header = TRUE, sep = ",")

#add substage column to the data
n.sbiv$substg <- rep(NA, nrow(n.sbiv))
n.sbiv$substg <- rep(NA, nrow(n.sbiv))

n.sbiv <- merge(n.sbiv, substg[, c('substage', 'sstg')], 
                by.x = 'age', by.y = 'substage', all.x = TRUE) # Add substages

n.sbiv$substg <- n.sbiv$sstg

#Delete the data over species (sp.), because the diversity based on species-level
o.sbiv <- n.sbiv[!grepl("sp\\..*", n.sbiv$spe), ]

#calculate taxonomic and ecological diversity
t.div <- divDyn(o.sbiv, tax = "gensp", bin="substg")#Species
f.div <- divDyn(o.sbiv, tax = "mol", bin="substg")#Function

# Add the substages to the diveristy
t.div$mid <- substg$mid[1:62]
f.div$mid <- substg$mid[1:62]


####Error bar####

# function

calculate_metrics <- function(samp) {
  
  # Initialize the result vector
  richness_indices <- rep(NA, 1000)
  
  n.row <- nrow(samp)
  n.col <- ncol(samp)
  
  for (i in 1:1000) {
    # random sample
    ssample <- sample(n.row, size = round(0.9 * nrow(samp)), replace = FALSE)
    bootsp <- samp[ssample, ]
    
    # Remove taxa that do not appear in the subsample
    subsamcolsum <- colSums(bootsp)
    if (any(subsamcolsum == 0)) {
      bootsp <- bootsp[, subsamcolsum != 0]
    }
    
    # Calculation of diversity indicators (e.g., species richness)

    if (ncol(bootsp) > 0) {
      richness_indices[i] <- as.numeric(specnumber(bootsp)[1])  # set the first value
    } else {
      richness_indices[i] <- NA  # if null set NA
    }
    
    
  }
  
  # Calculation of summary statistics for species richness
  richness_mean <- mean(richness_indices, na.rm = TRUE)
  richness_sd <- sd(richness_indices, na.rm = TRUE)
  richness_max <- max(richness_indices, na.rm = TRUE)
  richness_min <- min(richness_indices, na.rm = TRUE)
  
  #output
  cat("Mean: ", richness_mean, "\n")
  cat("Standard Deviation: ", richness_sd, "\n")
  cat("Max: ", richness_max, "\n")
  cat("Min: ", richness_min, "\n")
  
  return(list(
    richness_mean = richness_mean,
    richness_sd = richness_sd,
    richness_max = richness_max,
    richness_min = richness_min
  ))
}


#Substage error bar

#Griesbachian
dat1 <- subset(o.sbiv, age == "Griesbachian")

samp1 <- create.matrix(dat1, tax.name = "gensp", locality = "sample",
                       abund.col = "abu", abund = TRUE)

e.samp1 <- create.matrix(dat1, tax.name = "mol", locality = "sample",
                         abund.col = "abu", abund = TRUE)

g_res <- calculate_metrics(samp1)

e.g_res <- calculate_metrics(e.samp1)

#Dienerian
dat2 <- subset(o.sbiv, age == "Dienerian")

samp2 <- create.matrix(dat2, tax.name = "gensp", locality = "sample",
                       abund.col = "abu", abund = TRUE)

e.samp2 <- create.matrix(dat2, tax.name = "mol", locality = "sample",
                         abund.col = "abu", abund = TRUE)

d_res <- calculate_metrics(samp2)

e.d_res <- calculate_metrics(e.samp2)

#Smithian
dat3 <- subset(o.sbiv, age == "Smithian")

samp3 <- create.matrix(dat3, tax.name = "gensp", locality = "sample",
                       abund.col = "abu", abund = TRUE)

e.samp3 <- create.matrix(dat3, tax.name = "mol", locality = "sample",
                         abund.col = "abu", abund = TRUE)

smi_res <- calculate_metrics(samp3)

e.smi_res <- calculate_metrics(e.samp3)

#Pelsonian
dat4 <- subset(o.sbiv, age == "Pelsonian")

samp4 <- create.matrix(dat4, tax.name = "gensp", locality = "sample",
                       abund.col = "abu", abund = TRUE)

e.samp4 <- create.matrix(dat4, tax.name = "mol", locality = "sample",
                         abund.col = "abu", abund = TRUE)

p_res <- calculate_metrics(samp4)

e.p_res <- calculate_metrics(e.samp4)

#Illyrian
dat5 <- subset(o.sbiv, age == "Illyrian")

samp5 <- create.matrix(dat5, tax.name = "gensp", locality = "sample",
                       abund.col = "abu", abund = TRUE)

e.samp5 <- create.matrix(dat5, tax.name = "mol", locality = "sample",
                         abund.col = "abu", abund = TRUE)

i_res <- calculate_metrics(samp5)

e.i_res <- calculate_metrics(e.samp5)

#Fassanian

dat6 <- subset(o.sbiv, age == "Fassanian")

samp6 <- create.matrix(dat6, tax.name = "gensp", locality = "sample",
                       abund.col = "abu", abund = TRUE)

e.samp6 <- create.matrix(dat6, tax.name = "mol", locality = "sample",
                         abund.col = "abu", abund = TRUE)

f_res <- calculate_metrics(samp6)

e.f_res <- calculate_metrics(e.samp6)

#Julian
dat7 <- subset(o.sbiv, age == "Julian")
#dat7 <- subset(o.sbiv, age == "Julian1")

samp7 <- create.matrix(dat7, tax.name = "gensp", locality = "sample",
                       abund.col = "abu", abund = TRUE)

e.samp7 <- create.matrix(dat7, tax.name = "mol", locality = "sample",
                         abund.col = "abu", abund = TRUE)

j_res <- calculate_metrics(samp7)

e.j_res <- calculate_metrics(e.samp7)


#Plot the diversity curve
###### Subsampling of samples- Species #####
samp <- o.sbiv
samp[] <- NA

abus <- o.sbiv$abu
k <- 1

for (i in 1:nrow(o.sbiv)) {
  samp[k:(k+abus[i]-1),] <- o.sbiv[i,]
  k <- k+abus[i]
}

samp <- samp[,-9]

##substg##

#Shareholder Quorum Subsampling (SQS)
valid_keep_values <- c(52, 62)

ssqs.div <- subsample(samp, tax = "gensp", bin="substg", iter=100, q=0.7, 
                      output="dist", keep= valid_keep_values, type="sqs")

samp$f.mol <- rep(NA, nrow(samp))
samp$f.mol <- as.factor(samp$mol)

e.ssqs.div <- subsample(samp, tax = "f.mol", bin="substg", iter=100, q=0.7, 
                        output="dist", keep= valid_keep_values, type="sqs")

###### Fig. 7 A & B ######
#create a pdf
pdf("Fig7A.pdf",  width = 10, height = 10, family = "ArialMT")

par(mar = c(5, 4, 4, 5))

tsplot(substg, boxes="short", boxes.col="col", 
       shading="stage", xlim=52:63, ylim=c(0,20), ylab="Species Richness", xlab = "Age (Ma)")
lines (t.div$mid, t.div$divRT, pch = 16, type = "b", col="#E41A1C", lwd = 4, cex = 2)

#Add error bar
x= t.div$mid
richness_sd = c(rep(0,51),g_res$richness_sd,d_res$richness_sd,smi_res$richness_sd,0,
                0,0,p_res$richness_sd,i_res$richness_sd,f_res$richness_sd,0,
                j_res$richness_sd,0)
y_sd = richness_sd

arrows(x, t.div$divRT - y_sd, x, t.div$divRT + y_sd, angle = 90, code = 3, length = 0.05, col = "#E41A1C", lwd = 2)

#global
ssqs.div$mid <- substg$mid[1:62]
e.ssqs.div$mid <- substg$mid[1:62]

lines (ssqs.div$mid[52:63], ssqs.div$divRT[52:63], pch = 16, type = "b", col="#FF7F00", lwd = 4, cex = 2)
shades(ssqs.div$mid, ssqs.div$divRT, col= "#FF7F00")


par(new = TRUE)

tsplot(substg, boxes="short", ylim=c(-1,12), xlim=52:63, plot.args=list(axes=FALSE))
axis(side=4)
mtext(4,text="Functional richness",line=3)

lines (f.div$mid, f.div$divRT, pch = 17, type = "b", col="blue", lwd = 4, cex = 2)

#add error bar
f.x= f.div$mid
richness_sd = c(rep(0,51),e.g_res$richness_sd,e.d_res$richness_sd,e.smi_res$richness_sd,0,
                0,0,e.p_res$richness_sd,e.i_res$richness_sd,e.f_res$richness_sd,0,
                e.j_res$richness_sd,0)
f.y_sd = richness_sd

arrows(f.x, f.div$divRT - f.y_sd, f.x, f.div$divRT + f.y_sd, angle = 90, code = 3, length = 0.05, col = "blue", lwd = 2)

#global
lines (e.ssqs.div$mid[52:63], e.ssqs.div$divRT[52:63], pch = 17, type = "b", col="purple", lwd = 4, cex = 2)
shades(e.ssqs.div$mid, e.ssqs.div$divRT, col= "purple")


legend("topright", bg="white", legend=c("Raw Species", "Raw Function", "SQS Species", "SQS Function"), 
       col=c("#E41A1C", "blue", "#FF7F00", "purple"),pch = c(16, 17, 16, 17), lwd=2, 
       #inset=c(0.02,0.02), 
       cex=1)

dev.off()

###Turnover diversity####

t.sbiv <- sbiv


t.sbiv <- t.sbiv %>%
  mutate(age = gsub("Julian1", "Julian", age))

t.sbiv <- t.sbiv %>%
  mutate(age = gsub("Julian2", "Tuvalian", age))

stgMin <- categorize(t.sbiv[ ,"age"], keys$stgInt)
stgMin <- as.numeric(stgMin)

# empty container
t.sbiv$stg <- rep(NA, nrow(t.sbiv))
# Add "stg"
t.sbiv$stg <- stgMin

#Delete the data over species (sp.)
t.o.sbiv <- t.sbiv[!grepl("sp\\..*", t.sbiv$spe), ]

#add substage column to the data
t.o.sbiv$substg <- rep(NA, nrow(t.o.sbiv))
t.o.sbiv$substg <- rep(NA, nrow(t.o.sbiv))

t.o.sbiv <- merge(t.o.sbiv, substg[, c('substage', 'sstg')], 
                  by.x = 'age', by.y = 'substage', all.x = TRUE) # Add substages

t.o.sbiv$substg <- t.o.sbiv$sstg

t.t.div <- divDyn(t.o.sbiv, tax = "gensp", bin="substg")#Species
t.f.div <- divDyn(t.o.sbiv, tax = "mol", bin="substg")#Function

#net diversification#
#p, origination; q, extinction; r, net diversification.

o <- t.t.div$oriPC[53:62]
e <- t.t.div$extPC[53:62]

d<-o-e # diversification rate

#ecology
e.o <- t.f.div$oriPC[53:62]
e.e <- t.f.div$extPC[53:62]

e.d<-e.o-e.e # diversification rate

#create a pdf

pdf("Fig7B.pdf",  width = 10, height = 10, family = "ArialMT")

par(mar = c(5, 4, 4, 5))

tsplot(substg, boxes="short", boxes.col="col", 
       shading="stage", xlim=52:63, ylim=c(-0.5,2.5), ylab="Per capita turnover rates", xlab = "Age (Ma)")

lines(substg$mid[53:62], t.t.div$extPC[53:62], pch = 5, type = "b", col="#A65628", lty = 2, lwd = 4, cex = 1.5)

#lines(substg$mid[53:62], t.t.div$oriPC[53:62], pch = 5, type = "b", col="#F781BF", lty = 2, lwd = 2, cex = 1.5)

#Function turnover
lines(substg$mid[53:62], t.f.div$extPC[53:62], pch = 6, type = "b", col="#FF7F00", lty = 2, lwd = 4, cex = 1.5)

#lines(substg$mid[53:62], t.f.div$oriPC[53:62], pch = 6, type = "b", col="purple", lty = 2, lwd = 2, cex = 1.5)

par(new = TRUE)

tsplot(substg, boxes="short", xlim=52:63,  ylim=c(-3, 3),  xlab = "Age (Ma)", plot.args=list(axes=FALSE))

axis(side=4)
mtext(4,text="Diversification rate",line=3)

#species
lines(substg$mid[53:62], d , pch = 23, type = "b", col="#E41A1C", lty = 1, lwd = 4, cex = 1.5)


#function
lines(substg$mid[53:62], e.d, pch = 25, type = "b", col="blue", lty = 1, lwd = 4, cex = 1.5)

legend("topright", bg="white",legend=c("Species extinction","Functional extinction", "Species diversifaction","Functional diversifaction"), 
       col=c("#A65628", "#FF7F00","#E41A1C","blue"),pch = c(5, 6, 19, 17), lwd=2, cex=1)



dev.off()



###################################


####################################################################

#####Test sample basis####

###sample>20####

samples <- create.matrix(sbiv, tax.name = "gensp", locality="sample",
                         abund.col = "abu", abund=TRUE)

# limit analysis to some minimum value
y <- which(colSums(samples) > 20)
samples <- samples[,y]

# delete empty rows (could also be used to exclude rare species)
x <- which(rowSums(samples)> 0)
samples <- samples[x,]
r <- vegdist(samples, "bray")

# now turn the matrix around for Q-mode analyses
t.samples <- t(samples)

# compute a distance matrix
q <- vegdist(t.samples, "bray")
# perform a cluster analysis 
clu <- hclust(q, method="ward.D2")
plot(clu, hang = -1, main="Species")
p.table <- prop.table(samples, 2)

#Ecologic occurrences by sample
#This function is for the taxonomic names, may not fit ecology
e.samples <- create.matrix(sbiv, tax.name = "mol", locality="sample",
                           abund.col = "abu", abund=TRUE)

# limit analysis to some minimum value
y <- which(colSums(e.samples) > 20)
e.samples <- e.samples[,y]
# delete empty rows (could also be used to exclude rate species)
x <- which(rowSums(e.samples)>0)
e.samples <- e.samples[x,]
e.r <- vegdist(e.samples, "bray")

# now turn the matrix around for Q-mode analyses
t.e.samples <- t(e.samples)

# compute a distance matrix
e.q <- vegdist(t.e.samples, "bray")

# perform a cluster analysis 
e.clu <- hclust(e.q, method="ward.D2")
plot(e.clu, hang = -1, main="Ecology")

e.p.table <- prop.table(e.samples, 2)

###Colored Cluster of South China  ###
dend <- as.dendrogram(clu)

i = 0

colLab <<- function(n) {
  if (is.leaf(n)) {
    a = attributes(n)
    ligne = match(a$label, sbiv[, 1])
    age = sbiv[ligne, 4]
    
    col_age = switch(
      age,
      Griesbachian = "#E41A1C",
      Dienerian = "#377EB8",
      Smithian = "#4DAF4A",
      Spathian = "#984EA3",
      Pelsonian = "#FFD700",
      Illyrian = "#FF7F00",
      Fassanian = "#00008B",
      Julian1 = "#A65628",
      Julian2 = "#A65628",
      "black"  # Default color
    )
    
    #Set label color
    attr(n, "nodePar") <- c(
      a$nodePar,
      list(
        lab.col = col_age,
        lab.font = 1,
        pch = NA
        #,lab.cex = 0.8
      )
    )
    
    #set node point color
    # attr(n, "edgePar") <- list(col=col_age)
    
  }
  return(n)
}

#Taxonomy
dend<- as.dendrogram(clu)


dL <- dendrapply(dend, colLab)


##### Fig S4 #####
pdf("FigS4A.pdf",  width = 10, height = 10, family = "ArialMT")


plot(dL, 
     #main = "Species Dendrogram", 
     horiz = TRUE)
dev.off()


# Ecology Cluster
e.dend <- as.dendrogram(e.clu)
e.dL <- dendrapply(e.dend, colLab)

pdf("FigS4B.pdf",  width = 10, height = 10, family = "ArialMT")
plot(e.dL,
     # main = "Function Dendrogram", 
     horiz = TRUE)
legend(x="topleft", legend=c("Griesbachian","Dienerian","Smithian",
                             "Spathian","Pelsonian","Illyrian", "Fassanian","Julian"),
       col=c("#E41A1C", "#377EB8", "#4DAF4A", 
             "#984EA3", "#FFD700", "#FF7F00","#00008B",  "#A65628"), 
       cex=1, lty=1,lwd=2)

dev.off()

# Taxonomic NMDS

m <- metaMDS(t.samples)
ordiplot(m)
text(m, "species", col="blue", cex=0.4)
text(m, labels = rownames(t.samples), cex = 0.6, pos = 3)

m$stress
stressplot(m)


### Eco NMDS ###
#m <- metaMDS(e.q)
e.m <- metaMDS(t.e.samples)
ordiplot(e.m)

text(e.m, labels = rownames(t.e.samples), cex = 0.6, pos = 3)
text(e.m, "species", col="blue", cex=0.6)

e.m$stress
stressplot(e.m)


#Group By Stage#
gri <- unique(sbiv$sample[sbiv$age=="Griesbachian"])
die <- unique(sbiv$sample[sbiv$age=="Dienerian"])
smi <- unique(sbiv$sample[sbiv$age=="Smithian"])
pel <- unique(sbiv$sample[sbiv$age=="Pelsonian"])
ill <- unique(sbiv$sample[sbiv$age=="Illyrian"])
lad <- unique(sbiv$sample[sbiv$age=="Fassanian"])
jul <- unique(sbiv$sample[sbiv$age=="Julian1"|sbiv$age=="Julian2"])

age2 <- rep("", nrow(m$points))
age2[row.names(m$points) %in% gri] <- "gri"
age2[row.names(m$points) %in% die] <- "die"
age2[row.names(m$points) %in% smi] <- "smi"
age2[row.names(m$points) %in% pel] <- "pel"
age2[row.names(m$points) %in% ill] <- "ill"
age2[row.names(m$points) %in% lad] <- "lad"
age2[row.names(m$points) %in% jul] <- "jul"

age2 <- factor(age2, levels = c("gri", "die", "smi",
                                "pel","ill","lad", "jul"))

##### Fig. S5 #####
#Taxonomy
pdf("FigS5A.pdf",  width = 10, height = 10, family = "ArialMT")

par(mar = c(4, 4, 4, 4))  

mdscol = c(rep("#984EA3", 3),  rep("#A65628", 9), rep("#E41A1C", 5), rep("#377EB8",7),
           rep( "#458B00", 3), "#E41A1C","#FF7F00", rep("#E6AB02", 6),  "#FF7F00")

#Plot blank figure
ordiplot(m, type="none")

#add species labels
orditorp(m,display="species",col= "grey60",
         air=0.01,cex=1)

#add sample labels
orditorp(m,display="sites",col= mdscol,
         air=0.01,cex=1.25)

points(m, display = "sites", col= mdscol, pch = 19, cex = 2)

ordihull(m, groups=age2, col=c("#E41A1C","#377EB8","#458B00",
                               "#E6AB02","#FF7F00", "#984EA3",
                               "#A65628","#F781BF"))

dev.off()


#Ecology
pdf("FigS5B.pdf",  width = 10, height = 10, family = "ArialMT")

par(mar = c(4, 4, 4, 4)) 

ordiplot(e.m, type="none")

#add species labels
orditorp(e.m,display="species",col= "grey60",
         air=0.01,cex=1)

#add sample labels
orditorp(e.m,display="sites",col= mdscol,
         air=0.01,cex=1.25)

points(e.m, display = "sites", col= mdscol, pch = 19,cex = 2)


ordihull(e.m, groups=age2, col=c("#E41A1C","#377EB8","#458B00",
                                 "#E6AB02","#FF7F00", "#984EA3",
                                 "#A65628","#F781BF"))

legend(x="topright", legend=c("Griesbachian","Dienerian","Smithian",
                              "Pelsonian","Illyrian", "Fassanian", "Julian"),
       col=c("#E41A1C","#377EB8","#458B00",
             "#E6AB02","#FF7F00", "#984EA3", "#A65628"), 
       cex=1, lty=1,lwd=2)

dev.off()


###sample>40####

samples <- create.matrix(sbiv, tax.name = "gensp", locality="sample",
                         abund.col = "abu", abund=TRUE)

# limit analysis to some minimum value
y <- which(colSums(samples) > 40)
samples <- samples[,y]

# delete empty rows (could also be used to exclude rare species)
x <- which(rowSums(samples)> 0)
samples <- samples[x,]
r <- vegdist(samples, "bray")

# now turn the matrix around for Q-mode analyses
t.samples <- t(samples)

# compute a distance matrix

q <- vegdist(t.samples, "bray")
# perform a cluster analysis 
clu <- hclust(q, method="ward.D2")
plot(clu, hang = -1, main="Species")
p.table <- prop.table(samples, 2)

#Ecologic occurrences by sample
#This function is for the taxonomic names, may not fit ecology
e.samples <- create.matrix(sbiv, tax.name = "mol", locality="sample",
                           abund.col = "abu", abund=TRUE)

# limit analysis to some minimum value
y <- which(colSums(e.samples) > 40)
e.samples <- e.samples[,y]
# delete empty rows (could also be used to exclude rate species)
x <- which(rowSums(e.samples)>0)
e.samples <- e.samples[x,]
e.r <- vegdist(e.samples, "bray")

# now turn the matrix around for Q-mode analyses
t.e.samples <- t(e.samples)

# compute a distance matrix

e.q <- vegdist(t.e.samples, "bray")

# perform a cluster analysis 
e.clu <- hclust(e.q, method="ward.D2")
plot(e.clu, hang = -1, main="Ecology")

e.p.table <- prop.table(e.samples, 2)

###Colored Cluster of South China  ###
dend <- as.dendrogram(clu)

i = 0

colLab <<- function(n) {
  if (is.leaf(n)) {
    a = attributes(n)
    ligne = match(a$label, sbiv[, 1])
    age = sbiv[ligne, 4]
    
    col_age = switch(
      age,
      Griesbachian = "#E41A1C",
      Dienerian = "#377EB8",
      Smithian = "#4DAF4A",
      Spathian = "#984EA3",
      Pelsonian = "#FFD700",
      Illyrian = "#FF7F00",
      Fassanian = "#00008B",
      Julian1 = "#A65628",
      Julian2 = "#A65628",
      "black"  # Default color
    )
    
    #Set label color
    attr(n, "nodePar") <- c(
      a$nodePar,
      list(
        lab.col = col_age,
        lab.font = 1,
        pch = NA
        #,lab.cex = 0.8
      )
    )
    
    #set node point color
    # attr(n, "edgePar") <- list(col=col_age)
    
  }
  return(n)
}

#Taxonomy
dend<- as.dendrogram(clu)

dL <- dendrapply(dend, colLab)

##### Fig. S6 #####
pdf("FigS6A.pdf",  width = 10, height = 10, family = "ArialMT")


plot(dL, 
     #main = "Species Dendrogram", 
     horiz = TRUE)
dev.off()


# Ecology Cluster
e.dend <- as.dendrogram(e.clu)
e.dL <- dendrapply(e.dend, colLab)

pdf("FigS6B.pdf",  width = 10, height = 10, family = "ArialMT")
plot(e.dL,
     # main = "Function Dendrogram", 
     horiz = TRUE)
legend(x="topleft", legend=c("Griesbachian","Dienerian","Smithian",
                             "Spathian","Pelsonian","Illyrian", "Fassanian", "Julian"),
       
       col=c("#E41A1C", "#377EB8", "#4DAF4A", 
             "#984EA3", "#FFD700", "#FF7F00","#00008B", "#A65628"), 
       cex=1, lty=1,lwd=2)

dev.off()


##### Fig. S7 #####
# Taxonomic NMDS

m <- metaMDS(t.samples)
ordiplot(m)
text(m, "species", col="blue", cex=0.4)
text(m, labels = rownames(t.samples), cex = 0.6, pos = 3)

m$stress
stressplot(m)


### Eco NMDS ###
#m <- metaMDS(e.q)
e.m <- metaMDS(t.e.samples)
ordiplot(e.m)

text(e.m, labels = rownames(t.e.samples), cex = 0.6, pos = 3)
text(e.m, "species", col="blue", cex=0.6)

e.m$stress
stressplot(e.m)


#Group By Stage#
gri <- unique(sbiv$sample[sbiv$age=="Griesbachian"])
die <- unique(sbiv$sample[sbiv$age=="Dienerian"])
smi <- unique(sbiv$sample[sbiv$age=="Smithian"])
pel <- unique(sbiv$sample[sbiv$age=="Pelsonian"])
ill <- unique(sbiv$sample[sbiv$age=="Illyrian"])
lad <- unique(sbiv$sample[sbiv$age=="Fassanian"])
jul <- unique(sbiv$sample[sbiv$age=="Julian1"|sbiv$age=="Julian2"])


age2 <- rep("", nrow(m$points))
age2[row.names(m$points) %in% gri] <- "gri"
age2[row.names(m$points) %in% die] <- "die"
age2[row.names(m$points) %in% smi] <- "smi"
age2[row.names(m$points) %in% pel] <- "pel"
age2[row.names(m$points) %in% ill] <- "ill"
age2[row.names(m$points) %in% lad] <- "lad"
age2[row.names(m$points) %in% jul] <- "jul"
age2 <- factor(age2, levels = c("gri", "die", "smi",
                                "pel", "ill", "lad", "jul"))

pdf("FigS7A.pdf",  width = 10, height = 10, family = "ArialMT")

par(mar = c(4, 4, 4, 4))  

mdscol = c(rep("#984EA3", 2),  rep("#A65628", 8), rep("#E41A1C", 4), rep("#377EB8",6),
           rep( "#458B00", 1), "#E41A1C", rep("#E6AB02", 6),  "#FF7F00")

#Plot blank figure
ordiplot(m, type="none")

#add species labels
orditorp(m,display="species",col= "grey60",
         air=0.01,cex=1)

#add sample labels
orditorp(m,display="sites",col= mdscol,
         air=0.01,cex=1.25)

points(m, display = "sites", col= mdscol, pch = 19, cex = 2)

ordihull(m, groups=age2, col=c("#E41A1C","#377EB8","#458B00",
                               "#E6AB02","#FF7F00", "#984EA3",
                               "#A65628","#F781BF"))

dev.off()


#Ecology

pdf("FigS7B.pdf",  width = 10, height = 10, family = "ArialMT")

par(mar = c(4, 4, 4, 4)) 

ordiplot(e.m, type="none")

#add species labels
orditorp(e.m,display="species",col= "grey60",
         air=0.01,cex=1)

#add sample labels
orditorp(e.m,display="sites",col= mdscol,
         air=0.01,cex=1.25)

points(e.m, display = "sites", col= mdscol, pch = 19,cex = 2)

ordihull(e.m, groups=age2, col=c("#E41A1C","#377EB8","#458B00",
                                 "#E6AB02","#FF7F00", "#984EA3",
                                 "#A65628","#F781BF"))

legend(x="topright", legend=c("Griesbachian","Dienerian","Smithian",
                              "Pelsonian","Illyrian", "Fassanian", "Julian"),
       col=c("#E41A1C","#377EB8","#458B00",
             "#E6AB02","#FF7F00", "#984EA3", "#A65628"), 
       cex=1, lty=1,lwd=2)

dev.off()


#####Test lithology basis####
lbiv <- subset(sbiv, lit=="siliciclastic")

samples <- create.matrix(lbiv, tax.name = "gensp", locality="sample",
                         abund.col = "abu", abund=TRUE)

# limit analysis to some minimum value
y <- which(colSums(samples) > 30)
samples <- samples[,y]

# delete empty rows (could also be used to exclude rare species)
x <- which(rowSums(samples)> 0)
samples <- samples[x,]
r <- vegdist(samples, "bray")

# now turn the matrix around for Q-mode analyses
t.samples <- t(samples)

# compute a distance matrix

q <- vegdist(t.samples, "bray")
# perform a cluster analysis 
clu <- hclust(q, method="ward.D2")
plot(clu, hang = -1, main="Species")
p.table <- prop.table(samples, 2)

#Ecologic occurrences by sample
#This function is for the taxonomic names, may not fit ecology
e.samples <- create.matrix(lbiv, tax.name = "mol", locality="sample",
                           abund.col = "abu", abund=TRUE)

# limit analysis to some minimum value
y <- which(colSums(e.samples) > 30)
e.samples <- e.samples[,y]
# delete empty rows (could also be used to exclude rate species)
x <- which(rowSums(e.samples)>0)
e.samples <- e.samples[x,]
e.r <- vegdist(e.samples, "bray")

# now turn the matrix around for Q-mode analyses
t.e.samples <- t(e.samples)

# compute a distance matrix

e.q <- vegdist(t.e.samples, "bray")

# perform a cluster analysis 
e.clu <- hclust(e.q, method="ward.D2")
plot(e.clu, hang = -1, main="Ecology")

e.p.table <- prop.table(e.samples, 2)

###Colored Cluster of South China  ###
dend <- as.dendrogram(clu)

i = 0

colLab <<- function(n) {
  if (is.leaf(n)) {
    a = attributes(n)
    ligne = match(a$label, sbiv[, 1])
    age = sbiv[ligne, 4]
    
    col_age = switch(
      age,
      Griesbachian = "#E41A1C",
      Dienerian = "#377EB8",
      Smithian = "#4DAF4A",
      Spathian = "#984EA3",
      Pelsonian = "#FFD700",
      Illyrian = "#FF7F00",
      Fassanian = "#00008B",
      Julian1 = "#A65628",
      Julian2 = "#A65628",
      "black"  # Default color
    )
    
    #Set label color
    attr(n, "nodePar") <- c(
      a$nodePar,
      list(
        lab.col = col_age,
        lab.font = 1,
        pch = NA
        #,lab.cex = 0.8
      )
    )
    
    #set node point color
    # attr(n, "edgePar") <- list(col=col_age)
    
  }
  return(n)
}

#Taxonomy
dend<- as.dendrogram(clu)


dL <- dendrapply(dend, colLab)


#####Fig. S8#####
pdf("FigS78.pdf",  width = 10, height = 10, family = "ArialMT")


plot(dL, horiz = TRUE)
dev.off()


# Ecology Cluster
e.dend <- as.dendrogram(e.clu)
e.dL <- dendrapply(e.dend, colLab)

pdf("FigS8B.pdf",  width = 10, height = 10, family = "ArialMT")
plot(e.dL, horiz = TRUE)
legend(x="topleft", legend=c("Griesbachian","Dienerian","Smithian","Spathian", 
                             "Pelsonian","Illyrian", "Fassanian", "Julian"),
       
       col=c("#E41A1C", "#377EB8", "#4DAF4A", "#FFD700", 
             "#FF7F00","#00008B", "#A65628"), 
       
       cex=1, lty=1,lwd=2)

dev.off()


# Taxonomic NMDS

m <- metaMDS(t.samples)
ordiplot(m)

m$stress
stressplot(m)

### Eco NMDS ###

e.m <- metaMDS(t.e.samples)

e.m$stress
stressplot(e.m)

#Group By Stage#
gri <- unique(sbiv$sample[sbiv$age=="Griesbachian"])
die <- unique(sbiv$sample[sbiv$age=="Dienerian"])
smi <- unique(sbiv$sample[sbiv$age=="Smithian"])
pel <- unique(sbiv$sample[sbiv$age=="Pelsonian"])
ill <- unique(sbiv$sample[sbiv$age=="Illyrian"])
lad <- unique(sbiv$sample[sbiv$age=="Fassanian"])
jul <- unique(sbiv$sample[sbiv$age=="Julian1"|sbiv$age=="Julian2"])

age2 <- rep("", nrow(m$points))
age2[row.names(m$points) %in% gri] <- "gri"
age2[row.names(m$points) %in% die] <- "die"
age2[row.names(m$points) %in% smi] <- "smi"
age2[row.names(m$points) %in% pel] <- "pel"
age2[row.names(m$points) %in% ill] <- "ill"
age2[row.names(m$points) %in% lad] <- "lad"
age2[row.names(m$points) %in% jul] <- "jul"
age2 <- factor(age2, levels = c("gri", "die", "smi",
                                "pel","ill","lad","jul"))


######## Fig. S9 #######
#Taxonomy
pdf("FigS9A.pdf",  width = 10, height = 10, family = "ArialMT")


mdscol = c(rep("#A65628", 8),  rep("#E41A1C", 4), rep("#377EB8", 6), 
           "#458B00","#E41A1C")


#Plot blank figure
ordiplot(m, type="none")

#add species labels
orditorp(m,display="species",col= "grey60", air=0.01,cex=1)

#add sample labels
orditorp(m,display="sites",col= mdscol, air=0.01,cex=1.25)

points(m, display = "sites", col= mdscol, pch = 19, cex = 2)

ordihull(m, groups=age2, col=c("#E41A1C","#377EB8","#458B00","#E6AB02"))
dev.off()

#Ecology

pdf("FigS9B.pdf",  width = 10, height = 10, family = "ArialMT")

par(mar = c(4, 4, 4, 4)) 

ordiplot(e.m, type="none")

#add species labels
orditorp(e.m,display="species",col= "grey60",
         air=0.01,cex=1)

#add sample labels
orditorp(e.m,display="sites",col= mdscol,
         air=0.01,cex=1.25)

points(e.m, display = "sites", col= mdscol, pch = 19,cex = 2)

ordihull(e.m, groups=age2, col=c("#E41A1C","#377EB8","#458B00",
                                 "#E6AB02","#FF7F00", "#984EA3", "#A65628",
                                 "#F781BF"))

legend(x="topright", legend=c("Griesbachian","Dienerian","Smithian",
                              "Pelsonian","Illyrian", "Fassanian", "Julian"),
       col=c("#E41A1C","#377EB8","#458B00",
             "#E6AB02","#FF7F00", "#984EA3", "#A65628"), 
       cex=1, lty=1,lwd=2)

dev.off()


#### PBDB Diversity ####

#input the data
# Set the working directory to where your data is stored

data_file <- "pbdb_data_TriaBiv.xlsx"
sheet_name <- "Data"
g.dat <- read_excel(paste0(file_path, "/", data_file), sheet = sheet_name)

#View(subset(g.dat, is.na(mol)))

##### Diversity Curve of World ####

####stage####

#curve by stage

g.div <- divDyn(g.dat, tax = "genus", bin="stg")#Species
e.div <- divDyn(g.dat, tax = "mol", bin="stg")#Function

# Add the substages
g.div$mid <- stages$mid[1:62]
e.div$mid <- stages$mid[1:62]

### Subsampling Curve of PBDB ###

#Shareholder Quorum Subsampling (SQS)
valid_keep_values <- c(50, 60)

p.sqs.div <- subsample(g.dat,  iter=100, q=0.7,tax = "genus", bin="stg", 
                       output= "dist", keep=valid_keep_values, type="sqs")

g.dat$f.mol <- rep(NA, nrow(g.dat))
g.dat$f.mol <- paste(g.dat$mol, "MOL", sep="")

ep.sqs.div <- subsample(g.dat, iter=100, q=0.7, tax = "f.mol", bin="stg", 
                        output= "dist", keep=valid_keep_values,type="sqs")

#p.rr.div$mid <- stages$mid[1:62]
p.sqs.div$mid <- stages$mid[1:62]
ep.sqs.div$mid <- stages$mid[1:62]

##### Fig. 7C #####
#create a pdf
pdf("Fig7C.pdf",  width = 10, height = 10, family = "ArialMT")

par(mar = c(5, 4, 4, 5))

#Plot the curve
tsplot(stages, boxes="short", boxes.col="col", 
       shading="stage", xlim=51:59, ylim=c(0,300), ylab="Genera Richness", xlab = "Age (Ma)")

#raw
lines (g.div$mid, g.div$divRT, pch = 16, type = "b", col="#E41A1C", cex= 1.5, lwd = 4)

#SQS
lines (p.sqs.div$mid[50:62], p.sqs.div$divRT[50:62], pch = 16, type = "b", col="#FF7F00", cex= 1.5 , lwd = 4)

shades(p.sqs.div$mid, p.sqs.div$divRT, col="#FF7F00")

par(new = TRUE)

tsplot(stages, boxes="short", ylim=c(10,22), xlim=51:59, plot.args=list(axes=FALSE))
axis(side=4)
mtext(4,text="Functional richness",line=3)

lines (e.div$mid, e.div$divRT, pch = 17, type = "b", col="blue", cex= 1.5, lwd = 4) 

lines (ep.sqs.div$mid[50:62], ep.sqs.div$divRT[50:62], pch = 17, type = "b", col="#984EA3", cex= 1.5, lwd = 4)

shades(ep.sqs.div$mid, ep.sqs.div$divRT, col="#984EA3")

legend("topleft", bg="white", legend=c("Raw Species", "Raw Function", "SQS Species", "SQS Function"), 
       col=c("#E41A1C", "blue","#FF7F00", "#984EA3"), pch = c(16,17,16,17), lwd=2, 
       #inset=c(0.02,0.02), 
       cex=1)

dev.off()

####Turnover diversity####

#diversification rate#
#p, origination; q, extinction; r, diversification rate.
p.o <- g.div$oriPC[50:62]
p.e <- g.div$extPC[50:62]
p.d<-p.o-p.e # diversification rate
p.d[is.na(p.d)] <- 0

p.e.o <- e.div$oriPC[50:62]
p.e.e <- e.div$extPC[50:62]
p.e.d<- p.e.o-p.e.e # diversification rate
p.e.d[is.na(p.e.d)] <- 0

##### Fig. 7D #####

pdf("Fig7D.pdf",  width = 10, height = 10, family = "ArialMT")

par(mar = c(5, 4, 4, 5))

tsplot(stages, boxes="short", boxes.col="col", 
       shading="stage", xlim=51:59, ylim=c(-0.2,1), ylab="Per capita turnover rates", xlab = "Age (Ma)")

lines(stages$mid[50:62], g.div$extPC[50:62], pch = 5, type = "b", col="#A65628", lty = 2, lwd = 4, cex = 1.5)
#lines(stages$mid[50:62], g.div$oriPC[50:62], pch = 5, type = "b", col="#F781BF", lty = 2, lwd = 2, cex = 1.5)
lines(stages$mid[50:62], e.div$oriPC[50:62], pch = 6, type = "b", col="#FF7F00", lty = 2, lwd = 4, cex = 1.5)
#lines(stages$mid[50:62], e.div$extPC[50:62], pch = 6, type = "b", col="purple", lty = 2, lwd = 2, cex = 1.5)

par(new = TRUE)

tsplot(stages, ylim=c(-1,0.6), xlim=51:59, plot.args=list(axes=FALSE))
axis(side=4, ylab="Diversification rate")

lines(stages$mid[50:62], p.d , pch = 23, type = "b", col="#E41A1C", lty = 1, lwd = 4, cex = 1.5)
lines(stages$mid[50:62], p.e.d, pch = 25, type = "b", col="blue", lty = 1, lwd = 4, cex = 1.5)


legend("topright", bg="white",
       legend=c("Species extinction", 
                #"Species origination",
                "Functional extinction", 
                #"Functional origination", 
                "Species net diversifaction","Functional net diversifaction"), 
       col=c("#A65628", 
             #"#F781BF",
             "#FF7F00",
             #"purple",
             "#E41A1C","blue"
             ),
       pch = c(5, 6, 
              # 5, 6, 
               23, 25), lwd=2, cex=1)

dev.off()

##### PBDB Ecology Change #####

# Define the desired order for mol

# Filter the data to include only stages (stg) from 51 to 59
filtered_data <- g.dat %>%
  filter(stg >= 51 & stg <= 59) %>%
  group_by(stg, mol) %>%
  summarise(genus_count = n_distinct(genus)) %>%  # Count distinct genera
  ungroup()

# Calculate the total number of genera for each stage (stg)
filtered_data <- filtered_data %>%
  group_by(stg) %>%
  mutate(total_genus = sum(genus_count)) %>%
  ungroup()

# Calculate the relative abundance
filtered_data <- filtered_data %>%
  mutate(relative_abundance = genus_count / total_genus)

mol_order <- c(1, 4, 7, 9, 12, 19, 21, 16, 20, 22, 2, 3, 5, 6, 15, 17, 18, 10, 11, 14)

# Convert mol to a factor with the specified levels
filtered_data$mol <- factor(filtered_data$mol, levels = mol_order)

# Arrange the data by mol and stg to calculate the changes in relative abundance
filtered_data <- filtered_data %>%
  arrange(mol, stg) %>%
  group_by(mol) %>%
  mutate(change = relative_abundance - lag(relative_abundance)) %>%  # Calculate the difference
  ungroup()

# Assign colors based on the change in relative abundance
# Red for increases (> 0.001), blue for decreases (< -0.001), grey for the start (NA change), and yellow for no significant change
filtered_data <- filtered_data %>%
  mutate(color = case_when(
    is.na(change) ~ "grey",       # Starting value, no previous data
    change > 0.001 ~ "red",        # Increase
    change < -0.001 ~ "blue",      # Decrease
    TRUE ~ "yellow"                # No significant change
  ))

#### Fig. 8 ####
# Create the PDF plot
pdf("Fig8.pdf",  width = 10, height = 10, family = "ArialMT")

# Plot with circles representing relative abundance and color for changes
ggplot(filtered_data, aes(x = mol, y = factor(stg), size = relative_abundance)) +
  geom_point(aes(fill = color), shape = 21, color = "black") +  # Use circles to represent relative abundance
  scale_fill_identity() +  # Use the colors directly
  scale_size_continuous(range = c(1, 10), name = "Relative Abundance") +  # Adjust circle size range
  labs(x = "Mode of Life", y = "Stage", size = "Relative Abundance", fill = "Change") +
  theme_minimal() +
  theme(
    text = element_text(size = 12),
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "right"
  )

dev.off()