GO terms were identified using BLAST2GO analysis71. In short, a single representative nucleotide gene sequence from each of the 12,977 gene clusters was selected and BLASTX was performed locally against the NR database (last updated February 2016). The BLASTX output file in XML format (E-value < = 1e-05) was then used to map the GO terms using BLAST2GO under default settings71. GO term enrichment analysis was performed for all the clade-specific gene clusters, genes gained along the terminal and internal branches, genes under recombination and genes under positive selection. GO enrichment analysis was tested using the GOEAST tool72 assuming our experiment was a customized microarray platform. Because there was only one species placed in Clade 4, we could not compute enrichment for this clade for positive selection analysis. The p-value of GO ID enrichment was calculated as the hypergeometric probability of getting a sample of genes (example: number of genes gained/under recombination/selection) under the null hypothesis that they were selected randomly from the total pool of 12,977 genes. To control error rates for multiple hypothesis testing, the p-values were adjusted using the Benjamini Hochberg method implemented in the software Q-value67 where a false discovery rate (FDR)-adjusted p-value < 0.05 was considered significant.