###
# ASTRAL Figures in Manuscript
###

# Produced by Sven BUERKI (svenbuerki@boisestate.edu), October 2025

### 
# Set wd
####
setwd("~/Documents/Current_projects/Turreaeae/Phylogenetic_analyses")

###
# Load package
###
library(ape)
library(phytools)

###
# Load User-Defined function
###
####
# Function to extract q1 to q3 and pp1 in ASTRAL tree in phylo class
####
nodeASTRALextractoR <- function(tree){
  # Prepare data.frame to store node data
  nodeInfoASTRAL <- as.data.frame(matrix(ncol = 5, nrow = Nnode(tree)))
  colnames(nodeInfoASTRAL) <- c("NodeID", "q1", "q2", "q3", "pp1")
  
  # Populate node ID data
  total_nodes <- Ntip(tree) + Nnode(tree)
  nodeInfoASTRAL$NodeID <- seq(from = Ntip(tree)+1, to = total_nodes, by = 1)
  
  # Populate q1
  ## Extract q2 values
  q1_values <- sapply(tree$node.label, function(label) {
    # First go at the data
    tmp <- strsplit(label, split = "q1=")[[1]][2]
    q1 <- as.numeric(strsplit(tmp, split = ";")[[1]][1]) * 100
    return(q1)
  })
  ## Update data.frame
  nodeInfoASTRAL$q1 <- as.vector(q1_values)
  
  # Populate q2
  ## Extract q2 values
  q2_values <- sapply(tree$node.label, function(label) {
    # First go at the data
    tmp <- strsplit(label, split = "q2=")[[1]][2]
    q2 <- as.numeric(strsplit(tmp, split = ";")[[1]][1]) * 100
    return(q2)
  })
  ## Update data.frame
  nodeInfoASTRAL$q2 <- as.vector(q2_values)
  
  # Populate q3
  ## Extract q3 values
  q3_values <- sapply(tree$node.label, function(label) {
    # First go at the data
    tmp <- strsplit(label, split = "q3=")[[1]][2]
    q3 <- as.numeric(strsplit(tmp, split = ";")[[1]][1]) * 100
    return(q3)
  })
  ## Update data.frame
  nodeInfoASTRAL$q3 <- as.vector(q3_values)
  
  # Populate pp1
  ## Extract pp1 values
  pp1_values <- sapply(tree$node.label, function(label) {
    # First go at the data
    tmp <- strsplit(label, split = "pp1=")[[1]][2]
    pp1 <- as.vector(strsplit(tmp, split = ";")[[1]][1])
    return(pp1)
  })
  ## Update data.frame
  nodeInfoASTRAL$pp1 <- as.vector(pp1_values)
  
  #Return data
  return(nodeInfoASTRAL)
}

###
# Load CSV with tip names
###
tip.names <- read.csv("translation_tips_Turraea_revMC.csv")
# Produce final names
tip.names$Final.Names <- gsub(" ", "_", paste(tip.names$Species.name, " - ", tip.names$Collectors, sep = ""))

###
# Load tree
###
#Please adjust path based on your computer
# - ASTRAL tree
tree <- ape::read.tree(file = "Turraea_18_AstralV2.tre")

###
# Determine MRCA of clade IV
###
clade_IV <- ape::getMRCA(tree, tip = c("Turraea_023_L001", "Turraea_015_L001"))

###
# Collapse tree based on clade_IV
###
# Get all descendants of clade_IV
descendants <- getDescendants(tree, clade_IV)
descendants.tips <- tree$tip.label[descendants]
descendants.tips <- descendants.tips[which(descendants.tips != "Turraea_023_L001")]

# Drop descendants
collapsed_tree <- drop.tip(tree, descendants.tips, subtree = FALSE)

# Rename collapsed node to clade ID
collapsed_tree$tip.label[which(collapsed_tree$tip.label == "Turraea_023_L001")] <- "Clade_IV"

###
# Include only 1 T. pulchella
###
collapsed_tree <- drop.tip(collapsed_tree, "Turraea_035_L001", subtree = FALSE)

###
# Collapse tree based on Outgroup taxa
###
ingroup <- ape::getMRCA(collapsed_tree, tip = c("Turraea_046_L001", "Gramble_15498_L001"))

# Get all descendants of ingroup
ingroup.tips <- getDescendants(collapsed_tree, ingroup)
ingroup.tips <- collapsed_tree$tip.label[ingroup.tips]
outgroup.tips <- collapsed_tree$tip.label[which(!collapsed_tree$tip.label %in% ingroup.tips)]

# Drop descendants
collapsed_tree <- drop.tip(collapsed_tree, outgroup.tips[-1], subtree = FALSE)

# Rename collapsed node to clade ID
collapsed_tree$tip.label[which(collapsed_tree$tip.label == outgroup.tips[1])] <- "Outgroup_taxa"

###
# Extract ASTRAL data for plotting
###
# Extract q1 to q3 and pp1
collapsed_tree_nodeASTRAL <- nodeASTRALextractoR(collapsed_tree)
# Replace NAs by 0s
#collapsed_tree_nodeASTRAL[is.na(collapsed_tree_nodeASTRAL)] <- 0

# Updated nodelabels with pp1
collapsed_tree$node.label <- collapsed_tree_nodeASTRAL$pp1

###
# Rename tips
###
# Create vector with new names
new_tips <- sapply(collapsed_tree$tip.label, function(tip) {
  if (tip %in% tip.names$idSequencing) {
    tip.names$Final.Names[match(tip, tip.names$idSequencing)]
  } else {
    tip  # keep original name
  }
})
# Rename tips in tree 
collapsed_tree$tip.label <- as.vector(new_tips)

###
# Load UDF for Fig 1
###
# cladebox<-function(tree,node,color=NULL,...){
#   if(is.null(color)) color<-make.transparent("yellow",0.2)
#   obj<-get("last_plot.phylo",envir=.PlotPhyloEnv)
#   h<-max(nodeHeights(tree))
#   parent<-tree$edge[which(tree$edge[,2]==node),1]
#   x0<-max(c(obj$xx[node]+obj$xx[parent])/2,obj$xx[node]-0.05*h)
#   x1<-obj$x.lim[2]
#   dd<-getDescendants(tree,node)
#   y0<-min(range(obj$yy[dd]))-0.5
#   y1<-max(range(obj$yy[dd]))+0.5
#   polygon(c(x0,x1,x1,x0),c(y0,y0,y1,y1),col=color,
#           border=0)
# }

cladebox <- function(tree, node, color = NULL, cex = 1, ...) {
  if (is.null(color)) color <- make.transparent("yellow", 0.2)
  obj <- get("last_plot.phylo", envir = .PlotPhyloEnv)
  h <- max(nodeHeights(tree))
  x0 <- max(obj$xx[1:Ntip(tree)])  # Start at max tip x-position (aligned tips)
  # Estimate label width: use longest tip label length scaled by cex
  max_label_width <- max(nchar(tree$tip.label)) * 0.3 * cex  # Approximate char width
  x1 <- x0 + max_label_width  # End at approximate right edge of labels
  dd <- getDescendants(tree, node)
  y0 <- min(range(obj$yy[dd])) - 0.5
  y1 <- max(range(obj$yy[dd])) + 0.5
  polygon(c(x0, x1, x1, x0), c(y0, y0, y1, y1), col = color, border = 0)
}

###
# Figure of simplified tree
###

###
#With clades I to IV
###
# Rm Clade_IV and Outgroup taxa to facilitate figure printing
collapsed_tree$tip.label[which(collapsed_tree$tip.label == "Clade_IV")] <- ""
collapsed_tree$tip.label[which(collapsed_tree$tip.label == "Outgroup_taxa")] <- ""
pdf("Figures/Fig1_main_simplified_tree_ASTRAL_piecharts_pp1.pdf")
ape::plot.phylo(collapsed_tree, cex = 0.5, show.node.label = TRUE, 
                align.tip.label = TRUE, use.edge.length = TRUE, font = 3, no.margin = FALSE)

# Add piecharts on nodes
ape::nodelabels(pie = collapsed_tree_nodeASTRAL[,2:4], piecol = c("blue", "orange", "grey"), adj = c(0.1,0.7), cex = 0.6)

# Clade I - cladebox
node <- ape::getMRCA(collapsed_tree, tip = c("Turraea_pulchella_-_Pelger_730", "Turraea_streyi_-_Strey_6876"))
cladebox(tree = collapsed_tree, node = node, color = make.transparent("red", 0.2), cex = 0.5)

# Clade I label - using text()
obj <- get("last_plot.phylo", envir = .PlotPhyloEnv)
x0 <- max(obj$xx[1:Ntip(collapsed_tree)])  # Aligned tip position
max_label_width <- max(nchar(collapsed_tree$tip.label)) * 0.02 * 0.5  # Same as in cladebox
x_label <- x0 + max_label_width + 6  # Place label just to the right of cladebox
y_label <- mean(range(obj$yy[getDescendants(collapsed_tree, node)]))  # Center vertically in clade
text(x = x_label, y = y_label, labels = "Clade I", cex = 1, srt = 90, adj = c(0.5, 0))


# Clade II - cladebox
node <- ape::getMRCA(collapsed_tree, tip = c("Nymania_capensis_-_Maurin_1096", "Nymania_capensis_-_Chase_270"))
cladebox(tree = collapsed_tree, node = node, color = make.transparent("blue", 0.2), cex = 0.5)

# Clade II label - using text()
obj <- get("last_plot.phylo", envir = .PlotPhyloEnv)
x0 <- max(obj$xx[1:Ntip(collapsed_tree)])  # Aligned tip position
max_label_width <- max(nchar(collapsed_tree$tip.label)) * 0.02 * 0.5  # Same as in cladebox
x_label <- x0 + max_label_width + 6  # Place label just to the right of cladebox
y_label <- mean(range(obj$yy[getDescendants(collapsed_tree, node)]))  # Center vertically in clade
text(x = x_label, y = y_label, labels = "Clade II", cex = 1, srt = 90, adj = c(0.5, 0))

# Clade III - cladebox
node <- ape::getMRCA(collapsed_tree, tip = c("Turraea_africana_-_Grossweiler_4914", "Naregamia_alata_-_Kanodia_89603"))
cladebox(tree = collapsed_tree, node = node, color = make.transparent("orange", 0.2), cex = 0.5)

# Clade III label - using text()
obj <- get("last_plot.phylo", envir = .PlotPhyloEnv)
x0 <- max(obj$xx[1:Ntip(collapsed_tree)])  # Aligned tip position
max_label_width <- max(nchar(collapsed_tree$tip.label)) * 0.02 * 0.5  # Same as in cladebox
x_label <- x0 + max_label_width + 6  # Place label just to the right of cladebox
y_label <- mean(range(obj$yy[getDescendants(collapsed_tree, node)]))  # Center vertically in clade
text(x = x_label, y = y_label, labels = "Clade III", cex = 1, srt = 90, adj = c(0.5, 0))

# Clade IV
text(x = x0+1.3, y = 10.9, labels = "Clade IV", cex = 1, srt = 0, adj = c(0.5, 0))

# Outgroup taxa
text(x = x0+2, y = 0.9, labels = "Outgroup taxa", cex = 1, srt = 0, adj = c(0.5, 0))
dev.off()

###
# Close up Clade IV
###

# Extract the subtree with all clade IV
clade_IV <- extract.clade(tree, node = clade_IV)
clade_IV <- ladderize(clade_IV, right = F)

###
# Clade IVa-d
###
###
# Collapse tree based on clade_IV
###
# Get all descendants of clade_IVfh
clade_IVeh <- ape::getMRCA(clade_IV, tip = c("1973_13967_L001", "Turraea_015_L001"))
descendants.eh <- getDescendants(clade_IV, clade_IVfh)
descendants.eh.tips <- clade_IV$tip.label[descendants.fh]
descendants.eh.tips <- descendants.eh.tips[which(descendants.eh.tips != "1973_13967_L001")]

# Drop descendants
coll_tree_IVad <- drop.tip(clade_IV, descendants.eh.tips, subtree = FALSE)

# Rename collapsed node to clade ID
coll_tree_IVad$tip.label[which(coll_tree_IVad$tip.label == "1973_13967_L001")] <- "Clade_IVeh"

###
# Extract ASTRAL data for plotting
###
# Extract q1 to q3 and pp1
collapsed_tree_nodeASTRAL_IVad <- nodeASTRALextractoR(coll_tree_IVad)

# Updated nodelabels with pp1
coll_tree_IVad$node.label <- round(as.numeric(collapsed_tree_nodeASTRAL_IVad$pp1),2)

###
# Rename tips
###
# Create vector with new names
new_tips.IVad <- sapply(coll_tree_IVad$tip.label, function(tip) {
  if (tip %in% tip.names$idSequencing) {
    tip.names$Final.Names[match(tip, tip.names$idSequencing)]
  } else {
    tip  # keep original name
  }
})
# Rename tips in tree 
coll_tree_IVad$tip.label <- as.vector(new_tips.IVad)

###
# Figure Clade IVad
###

cladebox <- function(tree, node, color = NULL, cex = 1, ...) {
  if (is.null(color)) color <- make.transparent("yellow", 0.2)
  obj <- get("last_plot.phylo", envir = .PlotPhyloEnv)
  h <- max(nodeHeights(tree))
  x0 <- max(obj$xx[1:Ntip(tree)])  # Start at max tip x-position (aligned tips)
  # Estimate label width: use longest tip label length scaled by cex
  max_label_width <- max(nchar(tree$tip.label)) * 1.5 * cex  # Approximate char width
  x1 <- x0 + max_label_width  # End at approximate right edge of labels
  dd <- getDescendants(tree, node)
  y0 <- min(range(obj$yy[dd])) - 0.5
  y1 <- max(range(obj$yy[dd])) + 0.5
  polygon(c(x0, x1, x1, x0), c(y0, y0, y1, y1), col = color, border = 0)
}

###
#With clades IVad
###
# Rm Clade_IV and Outgroup taxa to facilitate figure printing
coll_tree_IVad$tip.label[which(coll_tree_IVad$tip.label == "Clade_IVeh")] <- ""

pdf("Figures/Fig2_clade_IVad_tree_ASTRAL_piecharts_pp1.pdf")
ape::plot.phylo(coll_tree_IVad, cex = 0.5, show.node.label = TRUE, 
                align.tip.label = TRUE, use.edge.length = TRUE, font = 3, no.margin = FALSE)

# Add piecharts on nodes
ape::nodelabels(pie = collapsed_tree_nodeASTRAL_IVad[,2:4], piecol = c("blue", "orange", "grey"), adj = c(-0.3,1), cex = 0.3)

# Clade IVa - cladebox
node <- ape::getMRCA(coll_tree_IVad, tip = c("Turraea_sp._nov._1_-_Leroy_16", "Turraea_sp._nov._3_-_Nusbaumer_1283"))
cladebox(tree = coll_tree_IVad, node = node, color = make.transparent("gold", 0.2), cex = 0.5)

# Clade IVa label - using text()
obj <- get("last_plot.phylo", envir = .PlotPhyloEnv)
x0 <- max(obj$xx[1:Ntip(coll_tree_IVad)])  # Aligned tip position
max_label_width <- max(nchar(coll_tree_IVad$tip.label)) * 0.02 * 0.5  # Same as in cladebox
x_label <- x0 + max_label_width + 35  # Place label just to the right of cladebox
y_label <- mean(range(obj$yy[getDescendants(coll_tree_IVad, node)]))  # Center vertically in clade
text(x = x_label, y = y_label, labels = "IVa", cex = 1, srt = 90, adj = c(0.5, 0))

# Clade IVb - cladebox
node <- ape::getMRCA(coll_tree_IVad, tip = c("Turraea_humberti_-_Léandri_2042", "Turraea_silvestris_-_Antilahimena_7148"))
cladebox(tree = coll_tree_IVad, node = node, color = make.transparent("green", 0.2), cex = 0.5)

# Clade IVb label - using text()
obj <- get("last_plot.phylo", envir = .PlotPhyloEnv)
x0 <- max(obj$xx[1:Ntip(coll_tree_IVad)])  # Aligned tip position
max_label_width <- max(nchar(coll_tree_IVad$tip.label)) * 0.02 * 0.5  # Same as in cladebox
x_label <- x0 + max_label_width + 35  # Place label just to the right of cladebox
y_label <- mean(range(obj$yy[getDescendants(coll_tree_IVad, node)]))  # Center vertically in clade
text(x = x_label, y = y_label, labels = "IVb", cex = 1, srt = 90, adj = c(0.5, 0))

# Clade IVc - cladebox
node <- ape::getMRCA(coll_tree_IVad, tip = c("Turraea_cadetii_-_Friedman_2607", "Turraea_ovata_-_Strasberg_752"))
cladebox(tree = coll_tree_IVad, node = node, color = make.transparent("violet", 0.2), cex = 0.5)

# Clade IVc label - using text()
obj <- get("last_plot.phylo", envir = .PlotPhyloEnv)
x0 <- max(obj$xx[1:Ntip(coll_tree_IVad)])  # Aligned tip position
max_label_width <- max(nchar(coll_tree_IVad$tip.label)) * 0.02 * 0.5  # Same as in cladebox
x_label <- x0 + max_label_width + 35  # Place label just to the right of cladebox
y_label <- mean(range(obj$yy[getDescendants(coll_tree_IVad, node)]))  # Center vertically in clade
text(x = x_label, y = y_label, labels = "IVc", cex = 1, srt = 90, adj = c(0.5, 0))

# Clade IVd - cladebox
node <- ape::getMRCA(coll_tree_IVad, tip = c("Turraea_sp._nov._5_-_Zjhra_290", "Turraea_virens_-_Ralimanana_260"))
cladebox(tree = coll_tree_IVad, node = node, color = make.transparent("brown", 0.2), cex = 0.5)

# Clade IVd label - using text()
obj <- get("last_plot.phylo", envir = .PlotPhyloEnv)
x0 <- max(obj$xx[1:Ntip(coll_tree_IVad)])  # Aligned tip position
max_label_width <- max(nchar(coll_tree_IVad$tip.label)) * 0.02 * 0.5  # Same as in cladebox
x_label <- x0 + max_label_width + 35  # Place label just to the right of cladebox
y_label <- mean(range(obj$yy[getDescendants(coll_tree_IVad, node)]))  # Center vertically in clade
text(x = x_label, y = y_label, labels = "IVd", cex = 1, srt = 90, adj = c(0.5, 0))

# Clade IVeh
text(x = x0+8, y = 52.9, labels = "Clade IVeh", cex = 1, srt = 0, adj = c(0.5, 0))
dev.off()


###
# Clade IVeh
###
# Get all descendants of clade_IVfh
clade_IVeh <- ape::getMRCA(clade_IV, tip = c("1973_13967_L001", "Turraea_015_L001"))

# Extract the subtree with all clade IV
coll_tree_IVeh <- extract.clade(clade_IV, node = clade_IVeh)
coll_tree_IVeh <- ladderize(coll_tree_IVeh, right = F)

###
# Extract ASTRAL data for plotting
###
# Extract q1 to q3 and pp1
collapsed_tree_nodeASTRAL_IVeh <- nodeASTRALextractoR(coll_tree_IVeh)

# Updated nodelabels with pp1
coll_tree_IVeh$node.label <- round(as.numeric(collapsed_tree_nodeASTRAL_IVeh$pp1),2)

###
# Rename tips
###
# Create vector with new names
new_tips.IVeh <- sapply(coll_tree_IVeh$tip.label, function(tip) {
  if (tip %in% tip.names$idSequencing) {
    tip.names$Final.Names[match(tip, tip.names$idSequencing)]
  } else {
    tip  # keep original name
  }
})
# Rename tips in tree 
coll_tree_IVeh$tip.label <- as.vector(new_tips.IVeh)

###
# Figure Clade IVeh
###

cladebox <- function(tree, node, color = NULL, cex = 1, ...) {
  if (is.null(color)) color <- make.transparent("yellow", 0.2)
  obj <- get("last_plot.phylo", envir = .PlotPhyloEnv)
  h <- max(nodeHeights(tree))
  x0 <- max(obj$xx[1:Ntip(tree)])  # Start at max tip x-position (aligned tips)
  # Estimate label width: use longest tip label length scaled by cex
  max_label_width <- max(nchar(tree$tip.label)) * 1.7 * cex  # Approximate char width
  x1 <- x0 + max_label_width  # End at approximate right edge of labels
  dd <- getDescendants(tree, node)
  y0 <- min(range(obj$yy[dd])) - 0.5
  y1 <- max(range(obj$yy[dd])) + 0.5
  polygon(c(x0, x1, x1, x0), c(y0, y0, y1, y1), col = color, border = 0)
}

pdf("Figures/Fig2_clade_IVeh_tree_ASTRAL_piecharts_pp1.pdf")
ape::plot.phylo(coll_tree_IVeh, cex = 0.5, show.node.label = TRUE, 
                align.tip.label = TRUE, use.edge.length = TRUE, font = 3, no.margin = FALSE)

# Add piecharts on nodes
ape::nodelabels(pie = collapsed_tree_nodeASTRAL_IVeh[,2:4], piecol = c("blue", "orange", "grey"), adj = c(-0.3,1), cex = 0.3)

# Clade IVe - cladebox
node <- ape::getMRCA(coll_tree_IVeh, tip = c("Turraea_pervillei_-_Jongkind_3615", "Turraea_fockei_-_Razafitsalama_1239"))
cladebox(tree = coll_tree_IVeh, node = node, color = make.transparent("plum1", 0.2), cex = 0.5)

# Clade IVe label - using text()
obj <- get("last_plot.phylo", envir = .PlotPhyloEnv)
x0 <- max(obj$xx[1:Ntip(coll_tree_IVeh)])  # Aligned tip position
max_label_width <- max(nchar(coll_tree_IVeh$tip.label)) * 0.02 * 0.5  # Same as in cladebox
x_label <- x0 + max_label_width + 47  # Place label just to the right of cladebox
y_label <- mean(range(obj$yy[getDescendants(coll_tree_IVeh, node)]))  # Center vertically in clade
text(x = x_label, y = y_label, labels = "IVe", cex = 1, srt = 90, adj = c(0.5, 0))

# Clade IVf - cladebox
node <- ape::getMRCA(coll_tree_IVeh, tip = c("Turraea_sambavensis_-_Service_Forestier_27118", "Turraea_sericea_-_Nusbaumer_1924"))
cladebox(tree = coll_tree_IVeh, node = node, color = make.transparent("grey", 0.2), cex = 0.5)

# Clade IVf label - using text()
obj <- get("last_plot.phylo", envir = .PlotPhyloEnv)
x0 <- max(obj$xx[1:Ntip(coll_tree_IVeh)])  # Aligned tip position
max_label_width <- max(nchar(coll_tree_IVeh$tip.label)) * 0.02 * 0.5  # Same as in cladebox
x_label <- x0 + max_label_width + 47  # Place label just to the right of cladebox
y_label <- mean(range(obj$yy[getDescendants(coll_tree_IVeh, node)]))  # Center vertically in clade
text(x = x_label, y = y_label, labels = "IVf", cex = 1, srt = 90, adj = c(0.5, 0))

# Clade IVg - cladebox
node <- ape::getMRCA(coll_tree_IVeh, tip = c("Humbertioturraea_labatii_-_Labat_2788", "Humbertioturraea_grandidieri_-_Jongkind_et_al._3605"))
cladebox(tree = coll_tree_IVeh, node = node, color = make.transparent("deeppink", 0.2), cex = 0.5)

# Clade IVg label - using text()
obj <- get("last_plot.phylo", envir = .PlotPhyloEnv)
x0 <- max(obj$xx[1:Ntip(coll_tree_IVeh)])  # Aligned tip position
max_label_width <- max(nchar(coll_tree_IVeh$tip.label)) * 0.02 * 0.5  # Same as in cladebox
x_label <- x0 + max_label_width + 47  # Place label just to the right of cladebox
y_label <- mean(range(obj$yy[getDescendants(coll_tree_IVeh, node)]))  # Center vertically in clade
text(x = x_label, y = y_label, labels = "IVg", cex = 1, srt = 90, adj = c(0.5, 0))

# Clade IVh - cladebox
node <- ape::getMRCA(coll_tree_IVeh, tip = c("Turraea_wakefieldii_-_Burrows_8596", "Turraea_pubescens_-_Frodin_843"))
cladebox(tree = coll_tree_IVeh, node = node, color = make.transparent("cadetblue1", 0.2), cex = 0.5)

# Clade IVh label - using text()
obj <- get("last_plot.phylo", envir = .PlotPhyloEnv)
x0 <- max(obj$xx[1:Ntip(coll_tree_IVeh)])  # Aligned tip position
max_label_width <- max(nchar(coll_tree_IVeh$tip.label)) * 0.02 * 0.5  # Same as in cladebox
x_label <- x0 + max_label_width + 47  # Place label just to the right of cladebox
y_label <- mean(range(obj$yy[getDescendants(coll_tree_IVeh, node)]))  # Center vertically in clade
text(x = x_label, y = y_label, labels = "IVh", cex = 1, srt = 90, adj = c(0.5, 0))
dev.off()
