An Atlas of the Speed of Copy Number Changes in Animal Gene Families and Its Implications

Statistics of gene families

The summary statistics of gene families is shown in Table 1. Since singletons can be the result of gene loss in a gene family, we included all singleton genes and considered them as gene families of size 1. There are altogether 23,713 gene families in the 10 species. The total number of genes and gene families in most species are around 20,000 and 10,000, respectively. The number of genes per family is about 2. Comparatively, the fruitfly has the smallest number of genes and gene families, whereas the mouse has the largest number of genes and gene families. The relatively smaller number of genes and families in the fruitfly may be due to its high rate of genome-wide gene deletion [16], [17], [18], so that the fruitfly's genome may be more compact than that of the other species. Interestingly, the zebrafish has the least number of gene families but the largest average gene family size, which may be due to its frequent genome duplications [19].

The relative CNC rate patterns

In order to get a distribution of the rates of CNCs, we used humans as a reference point to calculate rates of change in the other nine species because the human genome is one of the best annotated genomes and it seems intuitively more natural for people to understand changes relative to humans. Moreover, as we do not focus on determining the direction of changes, the choice of reference species should not affect our subsequent inference.

We calculated the rates of CNCs () with respect to humans for all of the 23,713 gene families (see Materials and Methods for the definition of ), and generated rate patterns by sorting the species based on their values. An example of a rate pattern is: . This pattern means that the fruitfly has the largest (rate of copy number changes), and Rs of all the other species are the same. Here, we included humans () in the pattern to help determine the direction of gene family size changes (i.e. expanded or shrunk).

We clustered all gene families into groups based on their rate patterns, calling each group a rate pattern group (RPG). We then sorted and indexed the RPGs by descending order of the number of gene families in each group. Thus, the smaller the index of a RPG, the more gene families the RPG contains. For instance, the RPG, , is the largest RPG containing 4,543 gene families, so it is indexed as . Altogether, we obtained 2,637 RPGs. Because the RPGs containing small number of gene families may be generated by noise and lack statistical power, we mainly focused on the 21 largest RPGs () that contain no less than 100 families. These 21 RPGs contain 17,412 (73.43%) gene families, indicating that the majority of the gene families are included in these RPGs (Table 2). Based on the evolutionary scenarios inferred from these RPGs, we further classified the 21 RPGs into three categories, 1) species-specific RPGs, including P0, P1, P2, P3, P4, P5, P6, P9, P10, P11, P17, P19, and P20; 2) group-specific RPGs, including P8, P12, P13, P14, P15, and P16; and 3) conserved RPG, P7, each of which is discussed as follows.

Table 2. The first 21 major RPGs.

Pattern ID	Pattern content	# of gene families
		4543 (19.16%)
		1984 (8.37%)
		1981 (8.35%)
		1642 (6.92%)
		869 (3.66%)
		800 (3.37%)
		798 (3.37%)
		768 (3.24%)
		559 (2.36%)
		474 (2.00%)
		385 (1.62%)
		385 (1.62%)
		377 (1.59%)
		374 (1.58%)
		333 (1.40%)
		262 (1.10%)
		216 (0.91%)
		192 (0.81%)
		182 (0.77%)
		164 (0.69%)
		124 (0.52%)
Total		17412 (73.43%)

Open in a new tab

Sorted by the number of genes in each pattern in descending order. When species are connected with “ = ”, they are sorted by alphabetic order.

^*

shortened pattern IDs, which are indexed based on the number of gene families in each pattern, the smaller the ID number is the more genes the pattern has.

GO annotation of RPGs through FBUA

It is interesting to know how the RPGs differ from each other in terms of biological functions and how functional difference is related to species difference. For this purpose, we annotated the RPGs into vectors of tip GO terms using the Full Bottom-Up Annotation (FBUA) approach that we proposed (see Materials and Methods for details). The complete annotation matrices are provided in Text S1, S2, and S3. Here we only focused on the 21 major RPGs.

There are altogether 6,964, 7,182, and 1,391 vector elements for biological processes, molecular function, and cellular component, respectively, for the 21 RPGs. Since we are mostly interested in identifying major differences between RPGs, we used the upper outlier of all elements in a RPG's vector to represent the RPG, that is, we selected the elements that are no less than in each vector, where is the 75% quartile and is the 25% quartile. As vector elements represent the overall probabilities that genes in a RPG are annotated to specific GO tip functions, this statistical cutoff allows us to pick out the highly enriched GO tip functions in each RPG. After removing the less frequent tips, we reassembled the remaining vectors and obtained three representative matrices that consist of 594, 355, and 165 GO tip elements for biological processes (BP), molecular functions (MF), and cellular components (CC) respectively.

To better reveal the relationships between the GO tips in the annotation matrix, we clustered the GO tips according to their similarities in the GO graphs [20]. At the same time, we also arranged the order of RPGs in the matrices based on the divergence times represented by the RPGs. Accordingly, the RPGs are ordered as

Figures 2– 4 show the full representative annotation matrices. Gray points mean that the values of matrix elements are less than the upper outlier. We expected that signals that are located closer in the matrices are basically more similar. However, overall, this expectation is very weakly reflected in BP and MF.

Vertically, non-mammals generally tend to have more representative data points (non-gray points) than mammals do especially in CC and partially in MF and BP. RPGs that represent mammals or non-mammals are more likely to be correlated within themselves, especially in CC (, spearman's correlation) and partially in MF and BP, indicating that there might be a significant change in some aspects of gene duplication between mammals and non-mammals. In MF, P7, P0, and P6 are significantly correlated with each other (, spearman's correlation), showing that conserved gene families may share a similar spectrum of functions with those fly specific gene families.

Horizontally, most of the patterns of the GO tips in BP and MF are not strictly clustered in a way that is compatible with their internal similarities in the GO graphs (showing by the trees in the left parts of Figures 2– 4). But, in some local regions this expectation holds. For example, in MF, magnesium, calcium, and zinc ion binding functions (GO:0000287; GO:0005509; zinc ion binding) matched perfectly between annotation matrices and GO DAGs. Notice that some functions are highly representative for many RPGs. More than 2/3 of the RPGs have functions that are related to transcription during RNA-mediated transposition for BP; to DNA clamp loader activity, DNA translocase activity, RNA polymerase II transcription factor activity, various binding (enhancers, magnesium, calcium, zinc, ATPs, and GTPs), DNA bending activity, negative regulation of interleukin-18 production, and actin homodimerization activity for MF; and to Nebenkern, macronucleus, micronucleus, apoplast, primary endosperm nucleus, root cap mucilage biosynthetic process, and root epithelial mucilage biosynthetic process for CC. This shows that gene families with ion and ATP or GTP binding functions are more likely to undergo heterogeneous rates of CNCs in multiple species.

Functional differences between species-specific RPGs

We summarized the highest weighted representative GO tip annotations (upper 90% quartile) in species-specific RPGs. In CC and MF, we did not find obvious species-specific differences in terms of functional categories in the top ranked GO tip annotations. In CC, almost all of the most representative GO tips are connected with nucleus, while in MF, almost all of them are related to ion binding, ATP binding or GTP binding.

In BP, the top weighted functions vary across species. About 2/3 of the differences are located in human vs. fly and human vs. zebrafish. Detailed list of functions is in Table 4. In fly vs. human RPGs (P0 and P6), we found two insect related pathways, ecdysone and chitin related pathways. But more pathways are sensory or nerve system related, such as chemical stimulus, neuropeptide signalling, corticospinal neuron axon decussation, Notch receptor processing, and sodium ion transport. We also noticed a mammalian specific pathway, alveolus development. Similar to previous studies, there are transcriptional regulatory related functions, such as the regulation of three RNA polymerase promoters, transcription anti-termination, and histone deubiquitination. RNA-mediated transposition is closely related to retroposition. There are also cancer related regulatory pathways, such as epidermal growth factor ligand processing, regulation of angiogenesis, Notch receptor processing, and Wnt signaling pathway. Maybe the most significant observation for speciation is the GO:0002077 that influences the acrosome matrix dispersal.

The differences between zebrafish and human in BP are mainly focused on the transcriptional regulatory related functions and signaling pathways, such as focal adhesion formation and protein kinase B signaling cascade. There are also special pathways: heart looping and ecdysone-mediated cell death, which are also found in the case of fly vs. human comparison.

In other species vs. human cases, although the overall weights of GO tips are not high, we can still find some interesting observations. In chicken, the highest ranked difference from human is the arginine biosynthetic process via ornithine, which is closely related to the urea cycle that shows great difference between land mammals and birds. In opossum, the top representative GO function is phototransduction, which may reflect the difference of retina between opossum and human. In cow, the function of ER to Golgi vesicle-mediated transport may have something to do with milking. In human, a neuropeptide signaling pathway is among the top ranked GO tips.

Functional differences between group-specific RPGs

Similar to the differences between species-specific RPGs, there is little difference between group-specific RPGs in CC and MF. All of the top ranked group-specific GO functions are related to nucleus in CC, while in MF most of them are related to zinc ion binding, ATP binding, and RNA polymerase II transcription factor activity-enhancer binding.

GO functions in BP differ greatly among species groups (Table 5). In fly vs. zebrafish (note that we showed that according to the pattern, P15 most likely reflects the difference between fly and zebrafish, rather than being compared with human), the most common top-ranked GO function differences are the regulation of transcription of all the three RNA polymerases (I, II and III). The retroposition related pathway is the top most difference between human and fly or zebrafish. We also observed functions that are explicitly related in phenotypic differences, such as blood vessels, neural crest cell migration, bone mineralization and bone morphogenetic protein pathway, and photoreceptor cell maintenance. The difference is also seen in regulation of some important pathways, such as transepithelial chloride transport, cell growth, and antisense RNA transcription.

Table 5. Top group-specific annotations (BP).

Groups	Pattern	Go id	Weight	GO Description
fly and zebrafish	P15	GO:0032199	52.20	transcription during RNA-mediated transposition
		GO:0001569	39.08	patterning of blood vessels
		GO:0007329	33.92	positive regulation of transcription from RNA polymerase II promoter by pheromones
		GO:0001755	32.47	neural crest cell migration
		GO:0007071	29.86	negative regulation of transcription from RNA polymerase III promoter, mitotic
		GO:0045494	29.42	photoreceptor cell maintenance
		GO:0046024	28.74	positive regulation of transcription from RNA polymerase III promoter, mitotic
		GO:0007069	27.89	negative regulation of transcription from RNA polymerase I promoter, mitotic
		GO:0046020	27.83	negative regulation of transcription from RNA polymerase II promoter by pheromones
		GO:0046018	26.76	positive regulation of transcription from RNA polymerase I promoter, mitotic
		GO:0030321	26.30	transepithelial chloride transport
		GO:0030502	26.04	negative regulation of bone mineralization
		GO:0001560	25.31	regulation of cell growth by extracellular stimulus
		GO:0030514	23.58	negative regulation of BMP signaling pathway
		GO:0060196	23.33	positive regulation of antisense RNA transcription
coldwarm	P8	GO:0007004	37.74	telomere maintenance via telomerase
blooded		GO:0031295	28.86	T cell costimulation
		GO:0032199	28.40	transcription during RNA-mediated transposition
		GO:0045842	27.47	positive regulation of mitotic metaphase/anaphase transition
		GO:0007257	21.80	activation of JNK activity
		GO:0007037	18.16	vacuolar phosphate transport
		GO:0050860	17.26	negative regulation of T cell receptor signaling pathway
		GO:0045086	17.07	positive regulation of interleukin-2 biosynthetic process
		GO:0002071	16.47	glandular epithelial cell maturation
		GO:0050434	15.24	positive regulation of viral transcription
non-mammals	P12	GO:0006120	34.89	mitochondrial electron transport, NADH to ubiquinone
mammals		GO:0032199	22.44	transcription during RNA-mediated transposition
		GO:0016578	19.51	histone deubiquitination
		GO:0031658	13.21	G1/S-specific negative regulation of cyclin-dependent protein kinase activity
		GO:0031661	13.20	G2/M-specific negative regulation of cyclin-dependent protein kinase activity
		GO:0030317	13.03	sperm motility
		GO:0042776	11.86	mitochondrial ATP synthesis coupled proton transport
		GO:0042777	11.86	plasma membrane ATP synthesis coupled proton transport
non-placental	P16	GO:0030890	22.66	positive regulation of B cell proliferation
placental		GO:0042104	18.04	positive regulation of activated T cell proliferation
		GO:0048304	17.04	positive regulation of isotype switching to IgG isotypes
		GO:0042523	13.13	positive regulation of tyrosine phosphorylation of Stat5 protein
murine	P14	GO:0006465	12.00	signal peptide processing

Open in a new tab

At the demarcation of cold and warm blooded animals, most of the top ranked differences are immune system related pathways: T cell costimulation, T cell receptor signaling pathway, interleukin-2 biosynthetic process, JNK activity, and viral transcription. The next common function differences include mitosis related functions such as telomere maintenance and mitotic metaphase/anaphase transition. Retroposition related pathway is ranked the third in this group.

Between non-mammal and mammals, the big difference in rates of CNCs involves many mitochondria related functions and energy and metabolism functions. Another major functional category in this group is cyclin-dependent protein kinase activity. Here, the retroposition related pathway (GO:0032199) is also highly ranked (second) in this group.

Between non-placental and placental animals, the big difference is related to immune systems regulating T cell proliferation, B cell proliferation, and IgG isotypes. In this group, a Stat5 protein related function is also highly ranked. And in the murine lineage, only a signal peptide processing pathway has a high rank.

Functional differences in the Conserved RPG

RPG P7 contains the gene families in which none of the species have differential rates of copy number changes. The gene families in this RPG are expected to reflect functions related to basic biological processes, so that the rates of CNCs are highly constrained by natural selection. Our results in Table 6 confirm this speculation. In BP, the highly ranked functions are related to the pathways involving histone, retroposition, RNA polymerase II, Gogi transport, splicing, snRNA, tRNA, cell cycle, docking to nuclear, and etc. Those processes are all basic and critical biological processes. In MF, the functions are similar to the highly ranked functions in other groups of RPGs.

Collecting Data

We analyzed 10 animal genomes including human (Homo sapiens), macaca (Macaca mulatta), mouse (Mus musculus), rat (Rattus norvegicus), dog (Canis familiaris), cow (Bos taurus), opossum (Monodelphis domestica), chicken (Gallus gallus), zebrafish (Danio rerio), and fruitfly (Drosophila melanogaster). The phylogeny of these species is shown in Figure 1 (adapted from Hedges[25]). We downloaded all the gene families for the ten species from ENSEMBL [26]. We required that each gene be a protein-coding gene with no premature stop codon, more than 100 amino acids, and have a known chromosomal location. We discarded genes on mitochondria. When a species does not have a gene family, we set the copy number of the family in this species to 0.

Calculating relative CNC rates and patterns

We used humans as the reference to calculate the CNC rate () of each gene family. Specifically,

(1)

where is the observed rate of family size change in family of species relative to the human, is the number of genes in family of species (when , the species is human), and is the divergence time in million years (MY) between species and human. The divergence times were obtained mainly from Hedges[25]. The opossum-eutheria divergence time (∼155 MY) was computed as the average of divergence time estimates in several studies [27], [28], [29]. The divergence times used in this study are shown in Figure 1. Note, the definition of does not require any predefined models, it is the observed CNC divided by the species divergence time.

For each gene family, we generated a rate pattern by sorting the species based on their values. Here, a pattern is an ordered sequence of species according to their relative CNC rates to human. After that, we clustered the gene families of the same pattern into groups, called rate pattern groups (RPG), and thus a rate pattern group contains all the gene families that have the same rate pattern. Table 2 shows the 21 largest RPGs (i.e. containing the highest number of gene families).

Annotating RPGs with gene ontology using FBUA

The GO terms of all the genes were extracted from Ensembl for the three categories–Biological Process (BP), Molecular Function (MF), and Cellular Component (CC). Since not all the genes have GO terms in each category and genes of the same family should perform more or less similar functions, we combined the GO terms of the genes belonging to the same gene family and removed duplicated GO terms. Then we pooled the family GO terms within each RPG. Note that, at this level, we preserved the frequencies of GO terms, which will be used as one of the parameters to determine the importance of a certain function.

The original collection of GO terms in each RPG is not suitable for representing the function spectrum of RPGs. This is because the GO terms are organized as directed acyclic graphs (DAGs) and GO terms at the same level (here level means the number of steps descending from the root) are not intrinsically comparable. To deal with this problem, for example, Karuppasamy et al.[15] annotated and compared GO terms at multiple levels. Lopez-Bigas et al.[10] simply used the GO terms that have more than 100 human genes. Both solutions are not satisfactory. In fact, Lopez-Bigas et al.[10]'s method may be significantly influenced by the systematic bias generated by the fact that certain genes tend to be better experimentally annotated than others. Also the method loses much resolution as the nodes that satisfy the criterion tend to be low level nodes (nodes that are very close to the root).

Here, we designed an annotation pipeline called Full Bottom-Up Annotation (FBUA) to deal with the aforementioned problems. We noticed that the tip GO terms (i.e. the leaves of GO graphs) contain all the information of the paths to the root. Assuming that the entire GO graphs are fully resolved, if we require that a gene with a certain GO term be fully annotated, it should be labeled as one of the tips in the current GO graphs. Therefore, we can use the tip GO terms plus their coalescent information to represent the whole DAGs. In this way, we can preserve all the information and do not have to be bothered by the bias from the selection of GO terms. However, it should be pointed out that this method can be only as good as the resolution of the entire GO graphs and carry the same bias as the GO graphs.

Formally, we assume that an internal node (GO term ) is equally likely to be annotated to any of its children, and all the paths from to the tip GO terms are independent of one another. The probability of the GO term being annotated to any of the tip GO terms , , can be computed as

(2)

is the path number, taking values from 1 to , assuming that there are different paths going from node to the tip , is the number of children that the parent node has in path , and is the collection of parental nodes involved in the path going from to , including .

Figure 5 shows an example of how we compute the probability. The likelihood for GO:0005496 to be annotated to the tip GO:0005497 is simply as the parent GO term has four children and is equally likely to be annotated to each one of them. The likelihood for GO:0005488 to be annotated to the tip GO:0005497 is computed as follows: there are altogether two paths from GO:0005488 to GO:0005497. The probability of GO:0005488 being annotated to GO:0005497 in path 1 is , and the probability of GO:0005488 being annotated to GO:0005497 in path 2 is . Assuming that two paths are independent, we have the final probability of GO:0005488 being annotated to GO:0005497 equal to .

Let be the set of tips of the entire GO graph and the collection of GO terms in . Then, can be annotated as a vector of GO tips , where is the frequency of in and each element in the vector represents the probability of the collection of G being annotated to a specific GO tip. Therefore, for each RPG, we can obtain an annotation vector .

Clustering and plotting annotation matrix

To better reveal the relationships between the GO tips in the annotation matrix, we calculated the Jiang and Conrath's pairwise similarity distance [20] between any two GO tips using the GOSim package [30]. We then constructed Neighbor-Joining trees [31] based on the distances using the Phylip program [32]. We plotted the annotation matrix using matrix2png program [33] and combined it with the GO tip trees.

Other data analyses

All the text parsing and processing procedures were done using a series of programs coded in OCAML language. All the statistical analyses were performed in R [34].