#library and dataset
library(ape)
data(bird.orders)

#slightly better than default plot
phylo = plot.phylo(bird.orders, edge.width = 2, label.offset = 1)

###
###OPTION 1
###
#first find out where are the groups, by adding y-axis and looking at coordinates
phylo = plot.phylo(bird.orders, edge.width = 3, label.offset = 1)
axis(2,at =seq(from =0, to = phylo$y.lim[2]),cex.axis = 0.5)
axis(1,at =seq(from =0, to = phylo$x.lim[2]),cex.axis = 0.5)
xmin = phylo$x.lim[1]
xmax = phylo$x.lim[2] - 10 #so the names are not boxed

#add boxes based on the y axis coordinates, the x-axis coord. are based on the phylogeny, and your groups of interest
rect(xmin,13.5,xmax,23.5, col = "red",density = 30) #group 1
rect(xmin,0.5,xmax,5.5,col = "blue",density = 30) #group 2
rect(xmin,6.5,xmax,12.5,col = "purple",density = 30) #group 3