diff --git a/Supplementary_appendices/Appendix_S9_Hoepke_et_al_ITS_plotting.r b/Supplementary_appendices/Appendix_S9_Hoepke_et_al_ITS_plotting.r
new file mode 100644
index 0000000000000000000000000000000000000000..c2b764442b9a207486d1e3a461dda81dfc882fde
--- /dev/null
+++ b/Supplementary_appendices/Appendix_S9_Hoepke_et_al_ITS_plotting.r
@@ -0,0 +1,244 @@
+###### R Script for paper "Phylogenetic and morphometric analysis of Plantago section 
+###### Coronopus (Plantaginaceae)" by Jannes Hoepke, Ladislav Mucina and Dirk C. Albach
+###### Date: 11 August 2015
+###### Date of revision: 25 February 2019
+###### Script is based on ...
+# GLOR, Rich (2008). "R Tip: Labeling Trees w/ Posterior Probability and Bootstrap Support". 
+#   http://treethinkers.blogspot.de/2008/10/labeling-trees-posterior-probability.html 
+#   (access date: 4th August 2015)
+# ... but was modified by Hoepke, J.
+
+# This script was created to plot the support values (summarised by SumTrees.py) on the "best" 
+# (= most likely) tree from GARLI. For doing so, please copy supplementary Appendix files 
+# S5 to S23 in one folder (S5-S11 are for ITS, S12-S17 for trnLF, and S18-S23 for cpDNA).
+
+
+
+
+################################# Initial set-up #######################################
+
+
+### Loading packages
+library(ape)
+library(geiger)
+library(phytools)
+
+
+### Print information about used package versions (not necessarily matching with citations in the 
+#   publication but it shows the program versions during my last successful run during revision 
+#   of this file)
+sink("Appendix_S10_Hoepke_et_al_tree_plotting_RSessionInfo.txt")
+print(sessionInfo())
+sink()
+
+
+### Create (sub-)directories, e.g. for storing plots
+system("mkdir Plots")
+
+
+
+
+########################### Make preparations ##################################
+
+
+### Some plot options
+# define cex-size for tiplabel
+tcex<-0.7
+# define cex-size for support
+scex<-0.6
+# define cex-size for tree line
+lcex<-1.5
+
+### Read in trees
+read.nexus("ITS_best.tre")->bestTree #FigTree output (rooted)
+read.nexus("ITS_mb_consensus.tre")->bayesTree #SumTrees.py output
+read.nexus("ITS_garli_consensus.tre")->bootTree #SumTrees.py output
+
+### Setting plot options
+opar<-par()
+par(mar=c(0,0,0,0),xpd=T)
+
+
+
+
+########################### Root all trees ##################################
+
+
+### Reroot bestTree
+plot(bestTree,cex=tcex)
+
+
+### Adjust bestTree via "rotate"
+nodelabels()
+rrbestTree<-rotate(bestTree, 61, polytom = c(1,2))
+plot(rrbestTree,show.node.label=F,cex=tcex)
+nodelabels()
+rrbestTree<-rotate(rrbestTree, 62, polytom = c(1,2))
+plot(rrbestTree,show.node.label=F,cex=tcex)
+nodelabels()
+rrbestTree<-rotate(rrbestTree, 63, polytom = c(2,3))
+plot(rrbestTree,show.node.label=F,cex=tcex)
+nodelabels()
+rrbestTree<-rotate(rrbestTree, 72, polytom = c(1,3))
+plot(rrbestTree,show.node.label=F,cex=tcex)
+nodelabels()
+rrbestTree<-rotate(rrbestTree, 73, polytom = c(1,2))
+plot(rrbestTree,show.node.label=F,cex=tcex)
+nodelabels()
+rrbestTree<-rotate(rrbestTree, 74, polytom = c(1,3))
+plot(rrbestTree,show.node.label=F,cex=tcex)
+nodelabels()
+rrbestTree<-rotate(rrbestTree, 79, polytom = c(3,13))
+plot(rrbestTree,show.node.label=F,cex=tcex)
+nodelabels()
+rrbestTree<-rotate(rrbestTree, 89, polytom = c(3,13))
+plot(rrbestTree,show.node.label=F,cex=tcex)
+
+
+# Create new labels for bestTree
+bestLabels<-read.csv("bestTree_tips.csv", header=T, sep=";")
+str(bestLabels)
+attach(bestLabels)
+bestLabels<-data.frame(lapply(bestLabels, as.character), stringsAsFactors=FALSE)
+#converts all colums from a dataframe into characters
+detach(bestLabels)
+attach(bestLabels)
+str(bestLabels)
+bestTips <- mixedFontLabel(new1," ",new2,", ",new3," ",new4, italic = 1:3, parenthesis = 7, sep="")
+
+
+### Reroot bootTree
+plot(bootTree,cex=tcex)
+nodelabels()
+
+i<-67
+rbootTree <- reroot(bootTree, i, position = 0.5 * bootTree$edge.length[which(bootTree$edge[,2] == i)])
+plot(rbootTree,show.node.label=T,cex=tcex)
+#Adjust bootTree via "rotate"
+nodelabels()
+rrbootTree<-rotate(rbootTree, 68, polytom = c(1,2))
+plot(rrbootTree,show.node.label=T,cex=tcex)
+nodelabels()
+rrbootTree<-rotate(rrbootTree, 72, polytom = c(1,2))
+plot(rrbootTree,show.node.label=T,cex=tcex)
+nodelabels()
+rrbootTree<-rotate(rrbootTree, 74, polytom = c(1,2))
+plot(rrbootTree,show.node.label=T,cex=tcex)
+nodelabels()
+rrbootTree<-rotate(rrbootTree, 75, polytom = c(18,19))
+plot(rrbootTree,show.node.label=T,cex=tcex)
+
+
+### Delete "Root" labeling
+rrbootTree$node.label[which(rrbootTree$node.label=="Root")]<-""
+plot(rrbootTree,show.node.label=T,cex=tcex)
+
+
+### Reroot bayesTree
+plot(bayesTree, show.node.label=T,cex=tcex)
+nodelabels()
+i<-67
+rbayesTree <- reroot(bayesTree, i, position = 0.5 * bayesTree$edge.length[which(bayesTree$edge[,2] == i)])
+plot(rbayesTree,show.node.label = T,cex=tcex)
+#Adjust bootTree via "rotate"
+nodelabels()
+rrbayesTree<-rotate(rbayesTree, 68, polytom = c(1,2))
+plot(rrbayesTree,show.node.label=T,cex=tcex)
+nodelabels()
+rrbayesTree<-rotate(rrbayesTree, 75, polytom = c(1,2))
+plot(rrbayesTree,show.node.label=T,cex=tcex)
+nodelabels()
+rrbayesTree<-rotate(rrbayesTree, 80, polytom = c(1,13))
+plot(rrbayesTree,show.node.label=T,cex=tcex)
+
+
+### Delete "Root" labeling
+rrbayesTree$node.label[which(rrbayesTree$node.label=="Root")]<-""
+plot(rrbayesTree,show.node.label=T,cex=tcex)
+
+
+
+
+########################### Prepare comparisons of trees ##################################
+
+### The getAllSubTrees function below is a necessary subfunction that atomizes a tree into 
+### each individual subclade and was provided compliments of Luke Harmon.
+getAllSubtrees<-function(phy, minSize=2) {
+  res<-list()
+  count=1
+  ntip<-length(phy$tip.label)
+  for(i in 1:phy$Nnode) {
+    l<-tips(phy, ntip+i)
+    bt<-match(phy$tip.label, l)
+    if(sum(is.na(bt))==0) {
+      st<-phy
+      }  
+    else st<-drop.tip(phy, phy$tip.label[is.na(bt)])
+    if(length(st$tip.label)>=minSize) {
+      res[[count]]<-st
+      count<-count+1
+      }
+  }
+  res
+}
+
+
+### The plotBootBayes function below plots both posterior probability and bootstrap values 
+### on each node of the consensus tree obtained from your Maximum-Likelihood analysis.
+### Bootstrap values will appear above and to the left of the node they support, whereas 
+### Bayesian posterior probabilies will appear below and to the left of the node.
+
+tiff("ITS_bestTree.tiff", width = 20, height = 25, units = 'cm', res = 300)
+par(mar=c(2,0,0,9),xpd=T)
+getAllSubtrees(rrbayesTree)->bayesSub
+getAllSubtrees(rrbootTree)->bootSub
+getAllSubtrees(rrbestTree)->bestSub
+bayesList<-matrix("",Nnode(rrbestTree),1)
+bootList<-matrix("",Nnode(rrbestTree),1)
+
+
+#The commands below compare all the subclades in the Bayes tree to all the subclades
+#in the bootstrap tree, and vice versa, and identifies all those clades that are 
+#identical.
+for(i in 1:Nnode(rrbestTree)) {
+  for(j in 1:Nnode(rrbayesTree)) {
+    match(bestSub[[i]]$tip.label[order(bestSub[[i]]$tip.label)], bayesSub[[j]]$tip.label[order(bayesSub[[j]]$tip.label)])->shared
+    match(bayesSub[[j]]$tip.label[order(bayesSub[[j]]$tip.label)], bestSub[[i]]$tip.label[order(bestSub[[i]]$tip.label)])->shared2
+    if(sum(is.na(c(shared,shared2)))==0) {
+      rrbayesTree$node.label[j]->bayesList[i]
+    }}}
+for(i in 1:Nnode(rrbestTree)) {
+  for(j in 1:Nnode(rrbootTree)) {
+    match(bestSub[[i]]$tip.label[order(bestSub[[i]]$tip.label)], bootSub[[j]]$tip.label[order(bootSub[[j]]$tip.label)])->shared
+    match(bootSub[[j]]$tip.label[order(bootSub[[j]]$tip.label)], bestSub[[i]]$tip.label[order(bestSub[[i]]$tip.label)])->shared2
+    if(sum(is.na(c(shared,shared2)))==0) {
+      rrbootTree$node.label[j]->bootList[i]
+    }}}
+
+plot(rrbestTree, cex=tcex, edge.width =lcex, show.tip.label=F) #Plots your Maximum-Likelihood consensus tree
+tiplabels(bestTips, adj = -0.02, frame = "none", col = "black", bg="white",cex=tcex)
+
+bootNodes1<-bootList
+bootNodes1[c(33)]<-""
+bootNodes2<-rep("",length(bootNodes1))
+bootNodes2[33]<-"100"
+
+
+bayesNodes1<-bayesList
+bayesNodes1[c(33)]<-""
+bayesNodes2<-rep("",length(bayesList))
+bayesNodes2[33]<-"1.00"
+
+
+nodelabels(bootNodes1, adj=c(1.2, -0.3), frame="n", cex=scex, font=1) 
+nodelabels(bayesNodes1, adj=c(1.12, 1.3), frame="n", cex=scex, font=1)
+nodelabels(bootNodes2, adj=c(2.3, -0.3), frame="n", cex=scex, font=1) 
+nodelabels(bayesNodes2, adj=c(2.1, 1.3), frame="n", cex=scex, font=1)
+
+add.scale.bar(x = 0, y = -1, length = NULL, cex=tcex)
+dev.off()
+# --> graphical output!
+
+
+par(opar)
+