The following code produces much the same sort of result using R

par(mfrow=c(1,2)) volef=c(47,12,3,2,0,0,0,0,0,0,0,1) volem=c(44,15,5,4,0,1,0,0,0,0,0,0) both=matrix(c(volef,volem), nr=2, byrow=TRUE) ylim=c(0,max(both[1,]+both[2,])) barplot(both, beside=FALSE, ylim=ylim, names=as.character(1:12), col=c('green','orange'), main='Bar Plot', xlab='number females/colony', ylab='Frequency') barplot(both,beside=TRUE, ylim=ylim, names=as.character(1:12), col=c('green','orange'), main='Bar Plot', xlab='number females/colony', ylab='Frequency')