Family-specific scaling laws in bacterial genomes

TL;DR

This study reveals that individual protein domain families in bacterial genomes follow specific scaling laws with genome size, indicating complex functional and evolutionary relationships beyond current annotations.

Contribution

It systematically characterizes family-specific scaling laws and introduces a heterogeneity score, deepening understanding of protein family evolution in bacterial genomes.

Findings

Protein families follow distinct scaling laws with genome size.

Functionally similar families often share similar scaling behaviors.

Heterogeneity score reveals diverse scaling patterns within categories.

Abstract

Among several quantitative invariants found in evolutionary genomics, one of the most striking is the scaling of the overall abundance of proteins, or protein domains, sharing a specific functional annotation across genomes of given size. The size of these functional categories change, on average, as power-laws in the total number of protein-coding genes. Here, we show that such regularities are not restricted to the overall behavior of high-level functional categories, but also exist systematically at the level of single evolutionary families of protein domains. Specifically, the number of proteins within each family follows family-specific scaling laws with genome size. Functionally similar sets of families tend to follow similar scaling laws, but this is not always the case. To understand this systematically, we provide a comprehensive classification of families based on their scaling properties. Additionally, we develop a quantitative score for the heterogeneity of the scaling of families belonging to a given category or predefined group. Under the common reasonable assumption that selection is driven solely or mainly by biological function, these findings point to fine-tuned and interdependent functional roles of specific protein domains, beyond our current functional annotations. This analysis provides a deeper view on the links between evolutionary expansion of protein families and the functional constraints shaping the gene repertoire of bacterial genomes.

Figures

Tables

Table 1. Table 1: Family scaling exponents can be associated to specific biological functions. Each cell in the table indicates the number of families that functional categories (rows) share with groups of families whose scaling exponents fall in pre-defined intervals (columns). The table also shows the Z-scores for a standard hypergeometric test (shown in green for over-representation and in red for under-representation, only | Z | > 𝑍 absent |Z|> 1.96 are shown).

Detailed function	${\beta }_i\le 0.6$		${\beta }_i\ge 1.4$	${\beta }_c\pm {\sigma }_{{\beta }_c}$
Translation	20($4.3$)	1($-3.7$)	0	$0.16\pm 0.03$
DNA replication/repair	11	7	0	$0.51\pm 0.07$
Transport	5	9	1	$1.1\pm 0.2$
Proteases	7	9	0	$0.9\pm 0.1$
Protein modification	8	1($-2.3$)	2	$1.06\pm 0.09$
Ion m/tr	11	3	3($-2.2$)	$1.3\pm 0.1$
Other enzymes	29	32	2	$1.04\pm 0.06$
Coenzyme m/tr	17($2.2$)	6	1	$0.85\pm 0.09$
Redox	4($-3.3$)	18($3.1$)	2	$1.2\pm 0.1$
Energy	11	7	0	$0.86\pm 0.09$
Nucleotide m/tr	16($3.1$)	3($-2.5$)	0	$0.53\pm 0.08$
Carbohydrate m/tr	4	8	0	$1.0\pm 0.2$
Transferases	5	11	1	$1.05\pm 0.07$
Amino acids m/tr	7	6	0	$0.8\pm 0.2$
DNA-binding	5	4	4($3.3$)	$1.5\pm 0.1$
Signal transduction	1($-2.7$)	5	5($5.0$)	$1.6\pm 0.2$
Unknown function	9	7	0	$0.98\pm 0.09$

Table 2. Supplementary Table S1: Symbols and codes used to identify functional categories.

symbol	category code	category name
	C	Energy
	E	Amino acids met./tr.
	F	Nucleotide met./tr.
	G	Carbohydrate met./tr.
	H	Coenzyme met./tr.
	J	Translation
	L	DNA replication
	LA	DNA binding
	O	Protein modification
	OA	Proteases
	P	Ion met./tr.
	RA	Redox
	RB	transferase
	RC	Other enzymes
	RF	Transport
	S	Unknown function
	T	Signal transductino

Table 3. Supplementary Table S2: Spearman correlations among family parameters. The table reports Spearman correlation coefficients between sets of family parameters, comparing biological/evolutionary and abundance properties. Each row describes biological parameters: for the Superfamily database we used the foldability (quantified by size-corrected contact order, SMCO  Debes2013 ) and the diversity of EC-numbers (quantifying the functional plasticity of a given family, see Section  S5 ) associated with families. For Pfam families, we considered the Hidden Markov Model sequence length (Hmm length) and the evolutionary rate (retrieved from  Ndhlovu2015 ). The parameters listed in columns are the exponent and prefactor of the family scaling law ( β i subscript 𝛽 𝑖 \beta_{i} and α i subscript 𝛼 𝑖 \alpha_{i} respectively), the mean family abundance calculated over all genomes ( ⟨ a i ⟩ delimited-⟨⟩ subscript 𝑎 𝑖 \langle a_{i}\rangle ) and the ratio between the average relative abundance (see definition of frequency in Section  S2 ) and family occurrence ( f i / o i subscript 𝑓 𝑖 subscript 𝑜 𝑖 f_{i}/o_{i} ). Relevant correlations are found for the diversity in EC numbers restricted to metabolic families and the scaling exponent β i subscript 𝛽 𝑖 \beta_{i} , as well as with the mean and relative family abundance. Family abundance and scaling exponent are also correlated (Spearman 0.72).

database	parameters	${\beta }_i$	${\alpha }_i$	$⟨a_i⟩$	$f_i/o_i$
SUPFAM	SMCO	$-0.06$	$0.07$	$0.04$	$0.04$
	EC numbers	$0.22$	$-0.14$	$0.40$	$0.35$
	(not met. families)
	EC numbers	$0.64$	$-0.50$	$0.77$	$0.74$
	(met. families)
PFAM	Hmm length	$0.13$	$-0.13$	$0.10$	$0.10$
	Ka/Ks	$0.11$	$-0.12$	$0.03$	$0.05$

Table 4. Supplementary Table S3: Scaling exponent of functional categories. The table reports the scaling exponent β c subscript 𝛽 𝑐 \beta_{c} of all functional categories examined, both for superfamilies (SUPFAM column) and clans (PFAM column). The error associated with the exponent is calculated as the root mean square deviation of the logarithm of the category abundance across all genomes from the estimated scaling law (see Methods).

cat. code	category name	${\beta }_c\pm {\sigma }_{{\beta }_c}$ (Supfam)	${\beta }_c\pm {\sigma }_{{\beta }_c}$ (Pfam)
A	RNA binding, met./tr.	$0.4\pm 0.1$	$0.21\pm 0.10$
B	Chromatin structure	$0.6\pm 0.3$	$--$
C	Energy	$0.8\pm 0.1$	$0.94\pm 0.09$
CA	E-transfer	$1.8\pm 0.3$	$1.73\pm 0.34$
CB	Photosynthesis	$0.9\pm 0.4$	$0.56\pm 0.28$
D	Cell cycle, Apoptosis	$1.2\pm 0.2$	$--$
E	Amino acids m/tr	$0.9\pm 0.2$	$1.09\pm 0.17$
EA	Nitrogen m/tr	$1.8\pm 0.2$	$--$
F	Nucleotide m/tr	$0.5\pm 0.1$	$0.70\pm 0.17$
G	Carbohydrate m/tr	$1.0\pm 0.2$	$1.38\pm 0.31$
GA	Polysaccharide m/tr	$1.1\pm 0.2$	$1.13\pm 0.20$
H	Coenzyme m/tr	$0.9\pm 0.1$	$1.27\pm 0.18$
HA	Small molecule binding	$0.9\pm 0.1$	$0.74\pm 0.06$
HD	Receptor activity	$1.3\pm 0.5$	$0.25\pm 0.20$
HE	Ligand binding	$0.3\pm 0.1$	$1.19\pm 0.31$
I	Lipid m/tr	$1.2\pm 0.2$	$1.31\pm 0.20$
IA	Phospholipid m/tr	$0.6\pm 0.3$	$0.77\pm 0.25$
J	Translation	$0.16\pm 0.03$	$0.17\pm 0.05$
K	Transcription	$0.9\pm 0.1$	$0.98\pm 0.19$
L	DNA replication/repair	$0.5\pm 0.1$	$0.63\pm 0.07$
LA	DNA-binding	$1.5\pm 0.1$	$1.49\pm 0.13$
LB	RNA processing	$0.2\pm 0.1$	$0.07\pm 0.12$
M	Cell envelope m/tr	$0.7\pm 0.4$	$--$
MA	Cell adhesion	$1.3\pm 0.3$	$1.49\pm 0.29$
N	Cell motility	$0.7\pm 0.3$	$0.48\pm 0.20$
O	Protein modification	$1.1\pm 0.1$	$0.99\pm 0.13$
OA	Proteases	$1.0\pm 0.1$	$1.01\pm 0.10$
OB	Kinases/phosphatases	$1.3\pm 0.2$	$1.02\pm 0.31$
P	Ion m/tr	$1.3\pm 0.1$	$1.04\pm 0.17$
Q	Secondary metabolism	$1.5\pm 0.2$	$2.01\pm 0.29$
R	General	$1.2\pm 0.2$	$1.61\pm 0.23$
RA	Redox	$1.2\pm 0.1$	$1.19\pm 0.13$
RB	Transferases	$1.1\pm 0.1$	$1.13\pm 0.09$
RC	Other enzymes	$1.1\pm 0.1$	$1.13\pm 0.06$
RD	Protein interaction	$1.0\pm 0.2$	$0.93\pm 0.19$
RF	Transport	$1.1\pm 0.2$	$1.10\pm 0.16$
S	Unknown function	$1.0\pm 0.1$	$0.89\pm 0.13$
SB	Toxins/defense	$1.1\pm 0.3$	$1.28\pm 0.29$
T	Signal transduction	$1.6\pm 0.2$	$1.74\pm 0.20$
TA	Other regulatory function	$1.1\pm 0.2$	$0.55\pm 0.24$

Table 5. Supplementary Table S4: Scaling exponent of superfamilies from the SUPERFAMILY database. The abundance of a (super)family scales as a power law of the genome size with family-dependent scaling exponents β i subscript 𝛽 𝑖 \beta_{i} . Each row corresponds to a domain family and shows its scaling exponent along with its error (see Methods) and the category to which the family belongs (category code). Families corresponding to the same functional category are ordered in decreasing order of abundance.

cat. code	family name	${\beta }_i\pm {\sigma }_{{\beta }_i}$
A	Alpha-L RNA-binding motif	$0.2\pm 0.1$
A	PUA domain-like	$0.7\pm 0.2$
C	6-phosphogluconate dehydrogenase C-terminal domain-like	$1.2\pm 0.2$
C	Glyceraldehyde-3-phosphate dehydrogenase-like,	$1.0\pm 0.2$
	C-terminal domain
C	Phosphoenolpyruvate/pyruvate domain	$1.0\pm 0.2$
C	SIS domain	$0.7\pm 0.2$
C	LeuD/IlvD-like	$0.9\pm 0.2$
C	Enolase C-terminal domain-like	$0.8\pm 0.2$
C	Transmembrane di-heme cytochromes	$0.5\pm 0.3$
C	Aconitase iron-sulfur domain	$0.5\pm 0.2$
C	Cytochrome c oxidase subunit I-like	$0.5\pm 0.2$
C	UDP-glucose/	$0.5\pm 0.2$
	GDP-mannose dehydrogenase C-terminal domain
C	Citrate synthase	$0.6\pm 0.2$
C	PEP carboxykinase-like	$0.2\pm 0.2$
C	Cytochrome c oxidase subunit III-like	$0.4\pm 0.2$
C	PK C-terminal domain-like	$0.2\pm 0.1$
C	Enzyme I of the PEP:sugar phosphotransferase	$0.4\pm 0.2$
	system HPr-binding (sub)domain
C	Cytochrome c oxidase subunit II-like, transmembrane region	$0.3\pm 0.2$
CA	Cytochrome c	$1.0\pm 0.5$
CA	Acyl-CoA dehydrogenase C-terminal domain-like	$2.0\pm 0.4$
CA	FMN-dependent nitroreductase-like	$0.8\pm 0.3$
CA	ISP domain	$1.2\pm 0.3$
CA	Sulfite reductase hemoprotein (SiRHP), domains 2 and 4	$0.9\pm 0.2$
CA	Succinate dehydrogenase/fumarate reductase flavoprotein,	$0.7\pm 0.2$
	catalytic domain
CB	PRC-barrel domain	$0.7\pm 0.3$
D	Rhodanese/Cell cycle control phosphatase	$1.1\pm 0.3$
E	ACT-like	$0.8\pm 0.2$
E	Tryptophan synthase beta subunit-like	$1.0\pm 0.2$
	PLP-dependent enzymes
E	Carbamate kinase-like	$0.5\pm 0.1$
E	PLP-binding barrel	$0.7\pm 0.1$
E	Glutamine synthetase/guanido kinase	$0.7\pm 0.2$
E	L-aspartase-like	$0.7\pm 0.2$
E	Diaminopimelate epimerase-like	$0.5\pm 0.2$
E	Alanine racemase C-terminal domain-like	$0.5\pm 0.2$
E	Aspartate/glutamate racemase	$0.3\pm 0.2$
E	Arginase/deacetylase	$0.8\pm 0.2$
E	Aspartate/ornithine carbamoyltransferase	$0.3\pm 0.1$
E	Serine metabolism enzymes domain	$0.3\pm 0.2$
E	Chorismate mutase II	$0.4\pm 0.2$
EA	RmlC-like cupins	$1.8\pm 0.2$
F	Ribonuclease H-like	$0.7\pm 0.3$
F	Adenine nucleotide alpha hydrolases-like	$0.9\pm 0.1$
F	Nucleotidylyl transferase	$0.17\pm 0.05$
F	PRTase-like	$0.4\pm 0.1$
F	Nucleotidyltransferase	$0.7\pm 0.2$
F	Pseudouridine synthase	$0.3\pm 0.1$
F	Ribulose-phoshate binding barrel	$0.5\pm 0.2$
F	Tetrahydrobiopterin biosynthesis enzymes-like	$0.4\pm 0.2$
F	Purine and uridine phosphorylases	$0.3\pm 0.2$
F	Nucleotidyltransferase substrate binding subunit/domain	$0.2\pm 0.2$
F	Nicotinate/Quinolinate PRTase C-terminal domain-like	$0.4\pm 0.2$
F	Nucleoside phosphorylase/	$0.3\pm 0.2$
	phosphoribosyltransferase catalytic domain
F	Nucleoside phosphorylase/	$0.4\pm 0.2$
	phosphoribosyltransferase N-terminal domain
F	NadA-like	$0.00\pm 0.03$
G	(Trans)glycosidases	$1.4\pm 0.4$
G	Aldolase	$0.9\pm 0.2$
G	Phosphoglucomutase, first 3 domains	$0.4\pm 0.1$
G	Galactose-binding domain-like	$0.8\pm 0.5$
G	Six-hairpin glycosidases	$1.2\pm 0.4$
G	Duplicated hybrid motif	$0.5\pm 0.2$
G	Xylose isomerase-like	$1.0\pm 0.3$
G	Carbohydrate phosphatase	$0.6\pm 0.2$
G	HIT-like	$0.6\pm 0.2$
G	Phosphoglucomutase, C-terminal domain	$0.4\pm 0.1$
G	PK beta-barrel domain-like	$0.9\pm 0.2$
G	HPr-like	$0.2\pm 0.2$
GA	UDP-Glycosyltransferase/glycogen phosphorylase	$1.0\pm 0.2$
GA	Pectin lyase-like	$1.3\pm 0.4$
GA	Glycosyl hydrolase domain	$0.8\pm 0.3$
GA	Barwin-like endoglucanases	$0.5\pm 0.2$
H	Glutathione synthetase ATP-binding domain-like	$0.9\pm 0.1$
H	Acyl-CoA dehydrogenase NM domain-like	$2.0\pm 0.4$
H	PreATP-grasp domain	$0.7\pm 0.1$
H	Single hybrid motif	$0.8\pm 0.2$
H	FMN-binding split barrel	$1.2\pm 0.2$
H	Riboflavin synthase domain-like	$1.0\pm 0.2$
H	Succinyl-CoA synthetase domains	$0.5\pm 0.2$
H	YrdC/RibB	$0.5\pm 0.2$
H	Molybdenum cofactor biosynthesis proteins	$0.6\pm 0.2$
H	Dihydrofolate reductase-like	$0.6\pm 0.2$
H	UROD/MetE-like	$0.5\pm 0.3$
H	Dihydropteroate synthetase-like	$0.4\pm 0.2$
H	Cobalamin (vitamin B12)-binding domain	$0.4\pm 0.3$
H	Activating enzymes of the ubiquitin-like proteins	$0.4\pm 0.2$
H	Nicotinate/Quinolinate PRTase N-terminal domain-like	$0.4\pm 0.2$
H	Glutamine synthetase, N-terminal domain	$0.6\pm 0.2$
H	Peptide deformylase	$0.2\pm 0.2$
H	RibA-like	$0.4\pm 0.2$
H	MoeA C-terminal domain-like	$0.2\pm 0.2$
H	ApbE-like	$0.3\pm 0.2$
HA	P-loop containing nucleoside triphosphate hydrolases	$0.71\pm 0.08$
HA	NAD(P)-binding Rossmann-fold domains	$1.4\pm 0.1$
HA	FAD/NAD(P)-binding domain	$1.3\pm 0.2$
HA	Thiamin diphosphate-binding fold (THDP-binding)	$0.9\pm 0.1$
HA	FAD-binding domain	$1.0\pm 0.2$
HA	Nucleotide-binding domain	$0.8\pm 0.2$
HA	Sensory domain-like	$0.7\pm 0.4$
HD	Methyl-accepting chemotaxis protein (MCP)	$1.0\pm 0.5$
	signaling domain
HD	PhoU-like	$0.3\pm 0.2$
HE	TGS-like	$0.3\pm 0.1$
I	Thioesterase/thiol ester dehydrase-isomerase	$1.5\pm 0.2$
I	Probable ACP-binding domain of	$1.0\pm 0.3$
	malonyl-CoA ACP transacylase
I	Creatinase/prolidase N-terminal domain	$0.5\pm 0.2$
I	Prokaryotic lipoproteins and	$0.4\pm 0.2$
	lipoprotein localization factors
IA	PLC-like phosphodiesterases	$0.6\pm 0.2$
J	Ribosomal protein S5 domain 2-like	$0.30\pm 0.07$
J	Translation proteins	$0.26\pm 0.06$
J	EF-G C-terminal domain-like	$0.30\pm 0.09$
J	Sm-like ribonucleoproteins	$0.8\pm 0.3$
J	Triger factor/SurA peptide-binding domain-like	$0.3\pm 0.2$
J	ValRS/IleRS/LeuRS editing domain	$0.03\pm 0.04$
J	Release factor	$0.2\pm 0.1$
J	L30e-like	$0.3\pm 0.2$
J	EF-Tu/eEF-1alpha/eIF2-gamma C-terminal domain	$0.5\pm 0.2$
J	S13-like H2TH domain	$0.3\pm 0.2$
J	NusB-like	$0.3\pm 0.1$
J	ClpS-like	$0.4\pm 0.1$
J	Ribosome binding protein Y (YfiA homologue)	$0.3\pm 0.1$
K	Tetracyclin repressor-like, C-terminal domain	$2.4\pm 0.3$
K	LexA/Signal peptidase	$0.6\pm 0.2$
K	Poly A polymerase C-terminal region-like	$0.2\pm 0.2$
K	GreA transcript cleavage protein, N-terminal domain	$0.2\pm 0.1$
K	CYTH-like phosphatases	$0.3\pm 0.1$
L	Nucleic acid-binding proteins	$0.31\pm 0.07$
L	DNA breaking-rejoining enzymes	$1.0\pm 0.3$
L	Nudix	$1.2\pm 0.2$
L	RuvA domain 2-like	$0.4\pm 0.1$
L	Restriction endonuclease-like	$0.8\pm 0.3$
L	DNA/RNA polymerases	$0.8\pm 0.2$
L	DNA-glycosylase	$0.7\pm 0.2$
L	DNase I-like	$0.8\pm 0.2$
L	DNA polymerase III clamp loader subunits,	$0.2\pm 0.1$
	C-terminal domain
L	Resolvase-like	$0.4\pm 0.4$
L	Uracil-DNA glycosylase-like	$0.6\pm 0.2$
L	GIY-YIG endonuclease	$0.4\pm 0.2$
L	DNA ligase/mRNA capping enzyme, catalytic domain	$0.7\pm 0.2$
L	HRDC-like	$0.4\pm 0.2$
L	N-terminal domain of MutM-like DNA repair proteins	$0.4\pm 0.2$
L	TRCF domain-like	$0.02\pm 0.03$
LA	Winged helix DNA-binding domain	$1.8\pm 0.2$
LA	Homeodomain-like	$2.2\pm 0.3$
LA	lambda repressor-like DNA-binding domains	$1.4\pm 0.3$
LA	C-terminal effector domain of	$1.8\pm 0.2$
	the bipartite response regulators
LA	Periplasmic binding protein-like I	$1.4\pm 0.4$
LA	Putative DNA-binding domain	$1.1\pm 0.2$
LA	Fatty acid responsive transcription factor FadR,	$1.9\pm 0.3$
	C-terminal domain
LA	Glucocorticoid receptor-like (DNA-binding domain)	$0.4\pm 0.2$
LA	TrpR-like	$0.4\pm 0.3$
LA	Ribbon-helix-helix	$0.8\pm 0.3$
LA	IHF-like DNA-binding proteins	$0.3\pm 0.3$
LA	ParB/Sulfiredoxin	$0.6\pm 0.3$
LA	KorB DNA-binding domain-like	$0.4\pm 0.3$
LB	EPT/RTPC-like	$0.3\pm 0.1$
M	OmpA-like	$0.9\pm 0.4$
MA	vWA-like	$1.2\pm 0.3$
MA	Pili subunits	$0.8\pm 0.4$
MA	PGBD-like	$0.7\pm 0.3$
MA	Hedgehog/DD-peptidase	$0.5\pm 0.2$
O	ATPase domain of HSP90 chaperone/	$1.5\pm 0.2$
	DNA topoisomerase II/ histidine kinase
O	GroES-like	$1.6\pm 0.3$
O	FKBP-like	$0.6\pm 0.2$
O	Chaperone J-domain	$0.5\pm 0.2$
O	Cyclophilin-like	$0.9\pm 0.2$
O	Double Clp-N motif	$0.6\pm 0.2$
O	HSP20-like chaperones	$0.6\pm 0.2$
O	GroEL equatorial domain-like	$0.2\pm 0.2$
O	GroEL apical domain-like	$0.2\pm 0.2$
O	Peptide methionine sulfoxide reductase	$0.2\pm 0.2$
O	GroEL-intermediate domain like	$0.2\pm 0.1$
OA	ClpP/crotonase	$1.0\pm 0.2$
OA	Zn-dependent exopeptidases	$1.0\pm 0.2$
OA	Metallo-dependent phosphatases	$1.0\pm 0.2$
OA	Metalloproteases (”zincins”), catalytic domain	$0.8\pm 0.3$
OA	LuxS/MPP-like metallohydrolase	$0.5\pm 0.3$
OA	Cysteine proteinases	$1.0\pm 0.3$
OA	Bacterial exopeptidase dimerisation domain	$1.1\pm 0.3$
OA	Trypsin-like serine proteases	$1.0\pm 0.3$
OA	Creatinase/aminopeptidase	$0.5\pm 0.1$
OA	HSP40/DnaJ peptide-binding domain	$0.3\pm 0.2$
OA	DPP6 N-terminal domain-like	$0.6\pm 0.4$
OA	Subtilisin-like	$0.8\pm 0.3$
OA	Rhomboid-like	$0.5\pm 0.2$
OA	Macro domain-like	$0.3\pm 0.2$
OA	Tricorn protease N-terminal domain	$0.4\pm 0.2$
OB	Protein kinase-like (PK-like)	$1.2\pm 0.3$
OB	PP2C-like	$0.8\pm 0.3$
OB	Phosphohistidine domain	$0.7\pm 0.2$
OB	Phosphotyrosine protein phosphatases I	$0.6\pm 0.2$
OB	Acylphosphatase/BLUF domain-like	$0.3\pm 0.2$
P	Periplasmic binding protein-like II	$1.6\pm 0.3$
P	MFS general substrate transporter	$1.4\pm 0.3$
P	Multidrug resistance efflux transporter EmrE	$1.3\pm 0.3$
P	HlyD-like secretion proteins	$1.5\pm 0.4$
P	Ferritin-like	$1.0\pm 0.2$
P	Cupredoxins	$1.1\pm 0.3$
P	Calcium ATPase, transduction domain A	$0.6\pm 0.2$
P	Calcium ATPase, transmembrane domain M	$0.6\pm 0.2$
P	TrkA C-terminal domain-like	$0.6\pm 0.3$
P	HMA, heavy metal-associated domain	$0.4\pm 0.2$
P	Band 7/SPFH domain	$0.5\pm 0.2$
P	Fe-S cluster assembly (FSCA) domain-like	$0.4\pm 0.2$
P	Voltage-gated potassium channels	$0.6\pm 0.2$
P	Magnesium transport protein CorA, transmembrane region	$0.5\pm 0.2$
P	CorA soluble domain-like	$0.5\pm 0.2$
P	Clc chloride channel	$0.3\pm 0.2$
Q	Dimeric alpha+beta barrel	$2.1\pm 0.3$
Q	Clavaminate synthase-like	$1.6\pm 0.3$
Q	Concanavalin A-like lectins/glucanases	$0.8\pm 0.4$
Q	Terpenoid synthases	$0.7\pm 0.2$
Q	Homo-oligomeric flavin-containing	$0.4\pm 0.2$
	Cys decarboxylases, HFCD
R	Bet v1-like	$1.5\pm 0.4$
R	Helical backbone metal receptor	$0.9\pm 0.3$
R	ADC-like	$1.0\pm 0.3$
R	ARM repeat	$0.6\pm 0.3$
R	Peripheral subunit-binding domain of 2-oxo	$0.3\pm 0.2$
	acid dehydrogenase complex
R	Pentein	$0.4\pm 0.2$
R	JAB1/MPN domain	$0.2\pm 0.2$
RA	Thioredoxin-like	$1.1\pm 0.2$
RA	4Fe-4S ferredoxins	$0.8\pm 0.4$
RA	Metallo-hydrolase/oxidoreductase	$1.1\pm 0.2$
RA	Glyoxalase/Bleomycin resistance protein/	$2.1\pm 0.3$
	Dihydroxybiphenyl dioxygenase
RA	ALDH-like	$1.5\pm 0.2$
RA	2Fe-2S ferredoxin-like	$1.0\pm 0.3$
RA	Flavoproteins	$0.9\pm 0.3$
RA	alpha-helical ferredoxin	$0.9\pm 0.3$
RA	FAD-linked reductases, C-terminal domain	$1.4\pm 0.3$
RA	Formate/glycerate dehydrogenase catalytic domain-like	$1.0\pm 0.2$
RA	NAD(P)-linked oxidoreductase	$1.3\pm 0.3$
RA	Isocitrate/Isopropylmalate dehydrogenase-like	$0.6\pm 0.2$
RA	Aminoacid dehydrogenase-like, N-terminal domain	$0.7\pm 0.1$
RA	FAD/NAD-linked reductases,	$0.7\pm 0.2$
	dimerisation (C-terminal) domain
RA	Formate dehydrogenase/DMSO reductase, domains 1-3	$1.0\pm 0.3$
RA	Ferredoxin reductase-like, C-terminal NADP-linked domain	$1.0\pm 0.3$
RA	Dehydroquinate synthase-like	$0.7\pm 0.3$
RA	Inosine monophosphate dehydrogenase (IMPDH)	$0.5\pm 0.2$
RA	Acid phosphatase/Vanadium-dependent haloperoxidase	$0.7\pm 0.2$
RA	FAD-linked oxidases, C-terminal domain	$1.0\pm 0.3$
RA	Sulfite reductase, domains 1 and 3	$0.6\pm 0.2$
RA	Succinate dehydrogenase/	$0.5\pm 0.2$
	fumarate reductase flavoprotein C-terminal domain
RA	LDH C-terminal domain-like	$0.3\pm 0.2$
RA	FAD-linked oxidoreductase	$0.6\pm 0.2$
RB	S-adenosyl-L-methionine-dependent methyltransferases	$0.9\pm 0.1$
RB	PLP-dependent transferases	$1.2\pm 0.1$
RB	Acyl-CoA N-acyltransferases (Nat)	$1.7\pm 0.2$
RB	Nucleotide-diphospho-sugar transferases	$0.9\pm 0.2$
RB	Class I glutamine amidotransferase-like	$1.0\pm 0.1$
RB	CoA-dependent acyltransferases	$1.0\pm 0.5$
RB	NagB/RpiA/CoA transferase-like	$1.1\pm 0.2$
RB	TK C-terminal domain-like	$0.7\pm 0.2$
RB	FabD/lysophospholipase-like	$1.1\pm 0.3$
RB	Tetrapyrrole methylase	$1.0\pm 0.3$
RB	Glycerol-3-phosphate (1)-acyltransferase	$0.7\pm 0.2$
RB	Formyltransferase	$0.6\pm 0.1$
RB	D-aminoacid aminotransferase-like PLP-dependent enzymes	$0.5\pm 0.2$
RB	4’-phosphopantetheinyl transferase	$0.5\pm 0.2$
RB	Methylated DNA-protein cysteine methyltransferase,	$0.8\pm 0.2$
	C-terminal domain
RB	Methylated DNA-protein cysteine methyltransferase domain	$0.6\pm 0.2$
RB	Homocysteine S-methyltransferase	$0.2\pm 0.2$
RC	alpha/beta-Hydrolases	$1.6\pm 0.3$
RC	Actin-like ATPase domain	$0.7\pm 0.1$
RC	HAD-like	$0.9\pm 0.2$
RC	Thiolase-like	$1.3\pm 0.3$
RC	Radical SAM enzymes	$0.7\pm 0.3$
RC	Acetyl-CoA synthetase-like	$1.9\pm 0.3$
RC	Metallo-dependent hydrolases	$1.3\pm 0.2$
RC	HD-domain/PDEase-like	$0.7\pm 0.3$
RC	beta-lactamase/transpeptidase-like	$0.9\pm 0.2$
RC	Trimeric LpxA-like enzymes	$0.7\pm 0.2$
RC	Lysozyme-like	$1.0\pm 0.2$
RC	Composite domain of metallo-dependent hydrolases	$1.4\pm 0.2$
RC	N-terminal nucleophile aminohydrolases (Ntn hydrolases)	$1.1\pm 0.2$
RC	Ribokinase-like	$0.9\pm 0.2$
RC	Alkaline phosphatase-like	$1.0\pm 0.3$
RC	DHS-like NAD/FAD-binding domain	$1.2\pm 0.2$
RC	Phospholipase D/nuclease	$0.8\pm 0.2$
RC	Glycoside hydrolase/deacetylase	$1.1\pm 0.2$
RC	Cytidine deaminase-like	$0.6\pm 0.1$
RC	LysM domain	$0.3\pm 0.4$
RC	SGNH hydrolase	$0.9\pm 0.3$
RC	PurM N-terminal domain-like	$0.2\pm 0.1$
RC	PurM C-terminal domain-like	$0.2\pm 0.1$
RC	Phosphoglycerate mutase-like	$0.9\pm 0.3$
RC	Galactose mutarotase-like	$0.6\pm 0.3$
RC	Carbon-nitrogen hydrolase	$0.8\pm 0.2$
RC	PHP domain-like	$0.6\pm 0.2$
RC	Enolase N-terminal domain-like	$0.9\pm 0.3$
RC	Quinoprotein alcohol dehydrogenase-like	$0.9\pm 0.3$
RC	all-alpha NTP pyrophosphatases	$0.5\pm 0.2$
RC	FAH	$1.3\pm 0.3$
RC	PFL-like glycyl radical enzymes	$0.3\pm 0.3$
RC	Amidase signature (AS) enzymes	$0.8\pm 0.3$
RC	Isochorismatase-like hydrolases	$1.2\pm 0.3$
RC	L,D-transpeptidase catalytic domain-like	$0.9\pm 0.3$
RC	Chorismate lyase-like	$0.8\pm 0.3$
RC	MoCo carrier protein-like	$0.5\pm 0.2$
RC	NAD kinase	$0.5\pm 0.2$
RC	ADC synthase	$0.4\pm 0.2$
RC	Folate-binding domain	$0.6\pm 0.2$
RC	AraD-like aldolase/epimerase	$0.7\pm 0.2$
RC	FMT C-terminal domain-like	$0.3\pm 0.2$
RC	IlvD/EDD N-terminal domain-like	$0.7\pm 0.2$
RC	Chelatase	$0.4\pm 0.2$
RC	Aminomethyltransferase beta-barrel domain	$0.6\pm 0.2$
RC	2-isopropylmalate synthase LeuA,	$10.3\pm 0.2$
	allosteric (dimerisation) domain
RC	CNF1/YfiH-like putative cysteine hydrolases	$0.2\pm 0.2$
RC	Nqo1 middle domain-like	$0.2\pm 0.2$
RC	beta-carbonic anhydrase, cab	$0.7\pm 0.2$
RC	N-acetylmuramoyl-L-alanine amidase-like	$0.3\pm 0.2$
RC	post-HMGL domain-like	$0.4\pm 0.2$
RC	Nqo1C-terminal domain-like	$0.2\pm 0.2$
RC	DmpA/ArgJ-like	$0.4\pm 0.2$
RC	Riboflavin kinase-like	$0.0\pm 0.0$
RC	LigT-like	$0.3\pm 0.2$
RD	TPR-like	$1.2\pm 0.4$
RD	FMN-linked oxidoreductases	$1.0\pm 0.1$
RD	Nqo1 FMN-binding domain-like	$0.3\pm 0.2$
RF	Multidrug efflux transporter AcrB transmembrane domain	$1.2\pm 0.3$
RF	Multidrug efflux transporter AcrB pore domain;	$1.2\pm 0.4$
	PN1, PN2, PC1 and PC2 subdomains
RF	Multidrug efflux transporter AcrB TolC docking domain;	$1.2\pm 0.4$
	DN and DC subdomains
RF	CBS-domain	$0.9\pm 0.2$
RF	ABC transporter transmembrane region	$0.7\pm 0.3$
RF	NTF2-like	$1.6\pm 0.3$
RF	Outer membrane efflux proteins (OEP)	$1.2\pm 0.3$
RF	ABC transporter involved in vitamin B12 uptake, BtuC	$0.8\pm 0.3$
RF	Rudiment single hybrid motif	$0.8\pm 0.2$
RF	Mechanosensitive channel protein MscS (YggB),	$0.6\pm 0.3$
	C-terminal domain
RF	Mechanosensitive channel protein MscS (YggB),	$0.6\pm 0.2$
	transmembrane region
RF	Proton glutamate symport protein	$0.3\pm 0.3$
RF	Ammonium transporter	$0.3\pm 0.2$
S	Sigma2 domain of RNA polymerase sigma factors	$1.4\pm 0.3$
S	ACP-like	$1.3\pm 0.4$
S	alpha/beta knot	$0.3\pm 0.1$
S	E set domains	$1.0\pm 0.3$
S	MOP-like	$1.0\pm 0.3$
S	PIN domain-like	$0.7\pm 0.3$
S	Anti-sigma factor antagonist SpoIIaa	$1.1\pm 0.3$
S	YjgF-like	$1.3\pm 0.2$
S	HCP-like	$0.4\pm 0.4$
S	ITPase-like	$0.4\pm 0.1$
S	MoaD/ThiS	$0.5\pm 0.2$
S	YbaK/ProRS associated domain	$0.6\pm 0.2$
S	Sporulation related repeat	$0.3\pm 0.2$
S	GatB/YqeY motif	$0.3\pm 0.1$
SB	AhpD-like	$1.5\pm 0.3$
T	CheY-like	$1.7\pm 0.2$
T	PYP-like sensor domain (PAS domain)	$2.0\pm 0.5$
T	Homodimeric domain of signal transducing histidine kinase	$1.6\pm 0.3$
T	Nucleotide cyclase	$1.6\pm 0.4$
T	GAF domain-like	$1.7\pm 0.3$
T	PDZ domain-like	$0.5\pm 0.2$
T	EAL domain-like	$0.9\pm 0.4$
T	cAMP-binding domain-like	$1.4\pm 0.3$
T	Histidine-containing phosphotransfer domain, HPT domain	$1.2\pm 0.3$
T	GlnB-like	$0.6\pm 0.2$
T	Mss4-like	$0.6\pm 0.3$
TA	Sigma3 and sigma4 domains of	$1.1\pm 0.2$
	RNA polymerase sigma factors
TA	OsmC-like	$1.0\pm 0.2$
TA	CinA-like	$0.3\pm 0.1$

Table 6. Supplementary Table S5: Scaling exponent of Pfam clans. The abundance of a clan scales as a power law of the genome size with family-dependent scaling exponents β i subscript 𝛽 𝑖 \beta_{i} . Each row of the table corresponds to a clan and shows its scaling exponent along with its error (see Methods) and the corresponding functional category (category code). Clans associated to the same functional category are ordered in decreasing order of abundance.

cat. code	clan name	${\beta }_i\pm {\sigma }_{{\beta }_i}$
A	S4 domain superfamily	$0.29\pm 0.15$
C	Pyruvate kinase-like TIM barrel superfamily	$1.16\pm 0.19$
C	6-phosphogluconate dehydrogenase C-terminal-like superfamily	$1.10\pm 0.16$
C	SIS domain fold	$0.76\pm 0.21$
C	Transmembrane di-heme cytochrome superfamily	$0.82\pm 0.24$
C	Enolase like TIM barrel	$1.13\pm 0.28$
C	PFK-like superfamily	$0.50\pm 0.23$
C	LeuD/IlvD-like	$0.52\pm 0.16$
CA	Cytochrome c superfamily	$0.88\pm 0.45$
CA	Acyl-CoA dehydrogenase, C-terminal domain-like	$1.95\pm 0.42$
CA	Rieske-like iron-sulphur domain	$1.15\pm 0.30$
CA	FMN-dependent nitroreductase-like	$0.78\pm 0.24$
CB	PRC-barrel like superfamily	$0.57\pm 0.28$
E	ACT-like domain	$0.70\pm 0.20$
E	gamma-glutamylcysteine synthetase/glutamine synthetase clan	$0.91\pm 0.23$
E	DAP epimerase superfamily	$0.51\pm 0.17$
E	Arginase/deacetylase superfamily	$0.74\pm 0.24$
E	Aspartate/glutamate racemase superfamily	$0.74\pm 0.23$
F	Ribonuclease H-like superfamily	$0.78\pm 0.29$
F	Nucleotidyltransferase superfamily	$0.73\pm 0.19$
F	PRPP synthetase-associated protein 1	$0.43\pm 0.15$
F	Tetrahydrobiopterin biosynthesis-like enzyme superfamily	$0.43\pm 0.20$
F	Nucleotidyltransferase substrate binding domain	$0.37\pm 0.26$
F	Purine and uridine phosphorylase superfamily	$0.44\pm 0.21$
F	dUTPase like superfamily	$0.22\pm 0.15$
G	Tim barrel glycosyl hydrolase superfamily	$1.32\pm 0.42$
G	Six-hairpin glycosidase superfamily	$1.13\pm 0.36$
G	Galactose-binding domain-like superfamily	$0.84\pm 0.46$
G	inositol polyphosphate 1 phosphatase like superfamily	$0.65\pm 0.24$
G	HIT superfamily	$0.55\pm 0.23$
GA	Glycosyl transferase clan GT-B	$0.98\pm 0.19$
GA	Pectate lyase-like beta helix	$1.49\pm 0.52$
GA	Glycosyl hydrolase domain superfamily	$0.87\pm 0.26$
GA	Double Psi beta barrel glucanase	$0.55\pm 0.19$
H	ATP-grasp superfamily	$0.88\pm 0.15$
H	Acyl-coenzyme A oxidase/dehydrogenase N-terminal	$1.86\pm 0.42$
H	Riboflavin synthase/Ferredoxin reductase FAD binding domain	$0.90\pm 0.23$
H	FMN-binding split barrel superfamily	$1.52\pm 0.27$
H	Dihydrofolate reductase-like	$0.65\pm 0.23$
H	Release factor superfamily	$0.29\pm 0.09$
H	Succinyl-CoA synthetase flavodoxin domain superfamily	$0.34\pm 0.18$
HA	P-loop containing nucleoside	$0.70\pm 0.07$
	triphosphate hydrolase superfamily
HA	PCMH-like FAD binding	$1.37\pm 0.25$
HD	PhoU-like superfamily	$0.41\pm 0.20$
HE	Ubiquitin superfamily	$1.24\pm 0.31$
I	HotDog superfamily	$1.60\pm 0.24$
I	Creatinase/prolidase N-terminal domain superfamily	$0.54\pm 0.20$
IA	PLC-like phosphodiesterases	$0.53\pm 0.23$
J	Ribosomal protein S5 domain 2-like superfamily	$0.32\pm 0.08$
J	Transcription elongation factor G C-terminal	$0.28\pm 0.11$
J	Helix-two-turns-helix superfamily	$0.20\pm 0.17$
J	DALR superfamily	$0.20\pm 0.12$
K	Peptidase clan SF	$0.71\pm 0.20$
L	OB fold	$0.45\pm 0.07$
L	PD-(D/E)XK nuclease superfamily	$0.58\pm 0.23$
L	NUDIX superfamily	$1.19\pm 0.21$
L	DNA breaking-rejoining enzyme superfamily	$0.97\pm 0.25$
L	His-Me finger endonuclease superfamily	$0.78\pm 0.28$
L	DNase I-like	$0.75\pm 0.24$
L	GIY-YIG endonuclease superfamily	$0.36\pm 0.22$
L	DNA/RNA ligase superfamily	$0.55\pm 0.20$
L	HRDC-like superfamily	$0.36\pm 0.16$
LA	Helix-turn-helix clan	$1.64\pm 0.14$
LA	Periplasmic binding protein like	$1.49\pm 0.38$
LA	Fatty acid responsive transcription factor FadR,	$1.98\pm 0.34$
	C-terminal domain
LA	lambda integrase N-terminal domain	$0.49\pm 0.22$
LA	MetJ/Arc repressor superfamily	$0.56\pm 0.34$
LA	ParB-like superfamily	$0.70\pm 0.26$
LA	IHF-like DNA-binding protein supewrfamily	$0.39\pm 0.28$
LB	EPT/RTPC-like superfamily	$0.33\pm 0.13$
MA	Ig-like fold superfamily (E-set)	$1.49\pm 0.43$
MA	von Willebrand factor type A	$1.17\pm 0.30$
MA	Pilus subunit	$0.83\pm 0.41$
MA	PGBD superfamily	$0.75\pm 0.33$
MA	Peptidase MD	$0.48\pm 0.24$
N	Flagellar motor switch family	$0.34\pm 0.30$
O	GroES-like superfamily	$1.56\pm 0.28$
O	FKBP-like superfamily	$0.62\pm 0.28$
O	Chaperone J-domain superfamily	$0.45\pm 0.24$
O	Cyclophilin-like superfamily	$0.74\pm 0.19$
O	HSP20-like chaperone superfamily	$0.65\pm 0.24$
OA	Peptidase clan MA	$0.97\pm 0.15$
OA	ClpP/Crotonase superfamily	$1.06\pm 0.21$
OA	Peptidase clan MH/MC/MF	$1.03\pm 0.19$
OA	Calcineurin-like phosphoesterase superfamily	$1.00\pm 0.18$
OA	Peptidase clan CA	$1.09\pm 0.24$
OA	LuxS/MPP-like metallohydrolase	$0.48\pm 0.26$
OA	Peptidase clan PA	$0.79\pm 0.26$
OA	MACRO domain superfamily	$0.33\pm 0.18$
OB	PP2C-like superfamily	$0.97\pm 0.34$
P	Ferritin-like Superfamily	$1.05\pm 0.20$
P	Multicopper oxidase-like domain	$0.93\pm 0.30$
P	SPFH superfamily	$0.34\pm 0.23$
P	SufE/NifU superfamily	$0.23\pm 0.16$
Q	Dimeric alpha/beta barrel superfamily	$2.01\pm 0.29$
R	Bet V 1 like	$1.50\pm 0.37$
R	Acetyl-decarboxylase like superfamily	$1.06\pm 0.29$
R	Helical backbone metal receptor superfamily	$0.81\pm 0.34$
R	GME superfamily	$0.44\pm 0.21$
RA	4Fe-4S ferredoxins	$0.94\pm 0.35$
RA	Thioredoxin-like	$1.16\pm 0.20$
RA	VOC superfamily	$2.11\pm 0.32$
RA	Metallo-hydrolase/oxidoreductase superfamily	$1.13\pm 0.17$
RA	ALDH-like superfamily	$1.64\pm 0.25$
RA	2Fe-2S iron-sulfur cluster binding domain	$1.05\pm 0.26$
RA	Transthyretin superfamily	$0.94\pm 0.53$
RA	Flavoprotein	$0.89\pm 0.27$
RA	Isocitrate/Isopropylmalate dehydrogenase-like superfamily	$0.64\pm 0.17$
RA	Formate/glycerate dehydrogenase	$0.99\pm 0.19$
	catalytic domain-like superfamily
RA	Ferredoxin / Ferric reductase-like NAD binding	$1.02\pm 0.26$
RA	Dehydroquinate synthase-like superfamily	$0.68\pm 0.29$
RA	FAD-linked oxidase C-terminal domain superfamily	$1.21\pm 0.25$
RA	Acid phosphatase/Vanadium-dependent haloperoxidase	$0.77\pm 0.24$
RA	LDH C-terminal domain-like superfamily	$0.23\pm 0.25$
RA	FAD-linked oxidoreductase	$0.52\pm 0.16$
RB	PLP dependent aminotransferase superfamily	$1.22\pm 0.13$
RB	N-acetyltransferase like	$1.70\pm 0.23$
RB	Glycosyl transferase clan GT-A	$0.91\pm 0.20$
RB	Class-I Glutamine amidotransferase superfamily	$0.99\pm 0.14$
RB	Isomerase,CoA transferase &	$1.04\pm 0.22$
	Translation initiation factor Superfamily
RB	CoA-dependent acyltransferase superfamily	$0.96\pm 0.39$
RB	Patatin/FabD/lysophospholipase-like superfamily	$1.12\pm 0.26$
RB	Acyltransferase clan	$0.67\pm 0.28$
RC	FAD/NAD(P)-binding Rossmann fold Superfamily	$1.12\pm 0.08$
RC	Alpha/Beta hydrolase fold	$1.54\pm 0.29$
RC	Actin-like ATPase Superfamily	$0.75\pm 0.14$
RC	Thiolase-like Superfamily	$1.22\pm 0.25$
RC	HAD superfamily	$0.79\pm 0.19$
RC	Amidohydrolase superfamily	$1.07\pm 0.17$
RC	Hexapeptide repeat superfamily	$0.58\pm 0.19$
RC	ANL superfamily	$1.91\pm 0.30$
RC	Serine beta-lactamase-like superfamily	$0.90\pm 0.18$
RC	HD/PDEase superfamily	$0.60\pm 0.27$
RC	NTN hydrolase superfamily	$1.06\pm 0.17$
RC	Ribokinase-like superfamily	$0.87\pm 0.20$
RC	Alkaline phosphatase-like	$1.01\pm 0.30$
RC	Lysozyme-like superfamily	$0.89\pm 0.26$
RC	LysM-like domain	$0.53\pm 0.34$
RC	DHS-like NAD/FAD-binding domain	$1.20\pm 0.18$
RC	Cytidine deaminase-like (CDA) superfamily	$0.60\pm 0.14$
RC	Phospholipase D superfamily	$0.55\pm 0.28$
RC	Histidine phosphatase superfamily	$0.77\pm 0.26$
RC	Glycoside hydrolase/deacetylase superfamily	$1.01\pm 0.25$
RC	Galactose Mutarotase-like superfamily	$0.86\pm 0.31$
RC	SGNH hydrolase superfamily	$0.91\pm 0.32$
RC	PFL-like glycyl radical enzyme superfamily	$0.26\pm 0.31$
RC	Enolase N-terminal domain-like superfamily	$1.06\pm 0.27$
RC	Fumarylacetoacetate hydrolase,	$1.34\pm 0.29$
	C-terminal domain, superfamily
RC	L,D-transpeptidase catalytic domain	$0.94\pm 0.30$
RC	Chorismate lyase/UTRA superfamily	$0.82\pm 0.31$
RC	MoCo carrier protein-like superfamily	$0.44\pm 0.16$
RC	Chelatase Superfamily	$0.45\pm 0.24$
RC	Fumarate reductase respiratory	$0.31\pm 0.19$
	complex transmembrane subunits
RD	Tetratrico peptide repeat superfamily	$1.07\pm 0.44$
RD	Common phosphate binding-site TIM barrel superfamily	$0.91\pm 0.11$
RF	Membrane and transport protein	$1.01\pm 0.30$
RF	ABC transporter membrane domain clan	$0.84\pm 0.26$
RF	NTF2-like superfamily	$1.07\pm 0.36$
S	Zinc beta-ribbon	$0.63\pm 0.20$
S	ACP-like superfamily	$1.68\pm 0.41$
S	SPOUT Methyltransferase Superfamily	$0.31\pm 0.11$
S	PIN domain superfamily	$0.87\pm 0.27$
S	STAS domain superfamily	$1.17\pm 0.30$
S	YjgF-like superfamily	$1.32\pm 0.24$
S	Phenylalanine- and lysidine-tRNA	$0.26\pm 0.17$
	synthetase domain superfamily
S	YqeY-like superfamily	$0.27\pm 0.14$
S	Maf/Ham1 superfamily	$0.30\pm 0.17$
SB	AhpD-like superfamily	$1.31\pm 0.28$
ST	Type III antifreeze and spore coat polysaccharide	$0.62\pm 0.23$
T	His Kinase A (phospho-acceptor) domain	$1.72\pm 0.19$
T	CheY-like superfamily	$1.72\pm 0.23$
T	PAS domain clan	$1.92\pm 0.44$
T	Nucleotide cyclase superfamily	$1.63\pm 0.44$
T	GAF domain-like	$2.15\pm 0.29$
T	PDZ domain-like peptide-binding superfamily	$0.59\pm 0.18$
T	GlnB-like superfamily	$0.71\pm 0.25$
T	Src homology-3 domain	$0.25\pm 0.36$

Equations

$$s_i=\frac{1}{1+\sqrt{L{S}_i}}$$

$$h_i=|{\beta }_c-{\beta }_i|,$$

$$H_c=\frac{1}{F_c}\sum _i{h}_i,$$

$$p_{i,c}\propto \{\begin{array}{ll}\frac{1}{N_g^{(i)}}\sum _{g^{\prime }=1}^{N_G^{(i)}}\frac{n_i^{g^{\prime }}}{{\sum }_{k\in c}{n}_k^{g^{\prime }}}& ,\text{mbox}if{n}_i^g\ne 0\\ 0& ,\text{mbox}if{n}_i^g=0.\end{array}$$

$$o_i=\frac{N_G^{(i)}}{N_G},$$

$$L{S}_i=\frac{1}{N_G^{(i)}}\sum _{g=1}^{N_G^{(i)}}\frac{{\left[y_{fit,i}-y_i\right]}^2}{y_{fit}},$$

$$y_{fit,i}={\alpha }_i+{\beta }_i\log \left(\sum _{i=1}^{N_F}{n}_i^g\right).$$

$$s_i=\frac{1}{1+\sqrt{L{S}_i}}$$

$$f_c(i)=\frac{1}{N_G^{(i)}}\sum _g\frac{n_i^g}{n_c^g},$$

$$p_c(i)=\{\begin{array}{ll}\frac{f_c(i)}{{\sum }_{i\in c}{f}_c(i)},& \text{mbox}if{n}_i^g\ne 0\\ 0,& \text{mbox}if{n}_i^g=0.\end{array}$$

$${\mu }_{bin,c}=F_{bin}\cdot \frac{F_c}{F_{tot}},$$

$${\sigma }_{bin,c}^2=F_{bin}\cdot \frac{F_c}{F_{tot}}\cdot \left(1-\frac{F_c}{F_{tot}}\right)\cdot \left(1-\frac{F_{bin}-1}{F_{tot}-1}\right),$$

$$Z_{bin,c}=\frac{x_{bin,c}-{\mu }_{bin,c}}{{\sigma }_{bin,c}},$$

Full text

I Introduction

As demonstrated by van Nimwegen vanNimwegen03 and confirmed by a series of follow-up studies Molina08 ; MOLINA09 ; cordero2009regulome ; Grilli2012 ; charoensawan2010genomic , striking quantitative laws exist for high-level functional categories of genes. Specifically, the number of genes within individual functional categories such as e.g. that of transcriptional regulators Stover2000 ; vanNimwegen03 ; Maslov2009 exhibit clear power-laws, when plotted as a function of genome size measured in terms of its number of protein-coding genes or, at a finer level of resolution, of their constitutive domains. In prokaryotes, such scaling laws appear well conserved across clades and lifestyles molina2009scaling , supporting the simple hypothesis that these scaling laws are universally shared by this group.

From the evolutionary genomics viewpoint Koonin2011 , these laws have been explained as a byproduct of specific “evolutionary potentials”, i.e., per-category-member rates of additions/deletions fixed in the population over evolution. As predicted by quantitative arguments, estimates of such rates correlate well with the category scaling exponents vanNimwegen03 ; Molina08 . A complementary point of view Maslov2009 ; Pang2011 ; Grilli2012 focuses on the existence of universal “recipes” determining ratios of proteins between different functions. Such recipes should mirror the “dependency structure” or network operating within genomes as well as other complex systems Pang2013 . According to this point of view the usefulness, and thus the occurence, of a given functional component depends on the presence of a set of other components, which are necessary for it to be operational.

Beyond functional categories, protein coding genes can be classified in “evolutionary families” defined by the homology of their sequences. Functional categories routinely contain genes from tens or more of distinct evolutionary families. The statistics of gene families also exhibits quantitative laws and regularities starting from a universal distribution of their per-genome abundance Huynen98 , explained by evolutionary models accounting for birth, death, and expansion of individual families Qian2001a ; Karev2004 ; Cosentino09 .

While some earlier work connects per-genome abundance statistics of families with functional scaling laws Grilli2012 , the link between functional category scaling and evolutionary expansion of gene families that build them remains relatively unexplored. Clearly, selective pressure is driven by functional constraints, and thus selection cannot in principle recognize families with identical functional roles. On the other hand, slight differences in the functional spectrum of different protein domains, and interdependency of different functions can make the scenario more complex. Thus, one central question is how the abundance of genes performing a specific function emerges from the evolutionary dynamics at the family level.

Two alternative extreme scenarios can be put forward: (i) the high-level scaling laws could emerge only at the level of functions, and be “combinatorially neutral” at the level of the evolutionary families building up a particular function, or, vice versa, (ii) they could be the result of the sum family-specific scaling laws. In the first scenario all or most of the families performing a particular function would be mutually interchangeable. In the second scenario, the evolutionary potentials would be family-specific and coincide with family evolutionary expansion rates, possibly emerging from the complex dependency structure cited above, and from fine-tuned functional specificity of distinct families. An intermediate possibility is that an interplay of constraints acts on both functional and evolutionary families. The first test for the feasibility of the second scenario is the existence of scaling laws for individual families. Here, focusing on bacteria, and using protein domains to define families, we present a clear evidence for family-specific scaling laws with genome size. We show that the abundance of the families follows power laws with genome size. Comparing functional categories with a suitable null model, we show that family-specific exponents may deviate significantly from the exponent of the associated functional category. We provide a comprehensive classification of families based on common scaling exponents, which recovers the known functional associations as well as revealing new ones, and may be used to detect possible misannotations. Finally, we develop quantitative tools to measure the heterogeneity of the scaling of families belonging to a given category or predefined group of families.

II MATERIALS AND METHODS

Data sources

We considered bacterial proteomes retrieved from the SUPERFAMILY (release 1.75 downloaded in October 2014, Gough01 ) and PFAM (release 27.0 downloaded in October 2014, bateman2004pfam ; Finn2014 ) database. Evolutionary families were defined from the domain assignments of 1535 superfamilies (SUPERFAMILY database) and 446 clans (PFAM database) on all protein sequences in completed genomes. We focused the analysis on the 1112 bacterial proteomes used as species reference in the SUPERFAMILY database. For the functional annotations of the SUPERFAMILY data, we considered annotation of SCOP domains as a scheme of 50 more detailed functional categories, mapped to 7 more general function categories, developed by C. Vogel vogel2006protein . PFAM clans were annotated on the same scheme of 50 functional categories, using the mapping of clans into superfamilies available from the PFAM website http://pfam.xfam.org/clan/browse#numbers finn2006pfam .

Data analysis

For each evolutionary domain family (or a functional category consisting of multiple evolutionary families), genome sizes (measured in the overall number of domains) were logarithmically binned. For each bin we calculated mean and standard deviation of the given family abundance (number of domains) within the bin. The estimated scaling exponent $\beta_{i}$ for family $i$ is the result of the non-linear least squares fitting of the binned data weighted by the standard error of family abundance. Genome size bins containing less than $10$ genomes were not taken into account. To filter out the data that, due to low-abundance or rare families, were affected by sampling problems, we considered three independent parameters, (i) the “occurrence”, i.e. the fraction of genomes where family $i$ is present, $o_{i}=N_{G}^{(i)}/N_{G}$, where $N_{G}$ is the total number of genomes in the sample, and $N_{G}^{(i)}$ is the number of genomes where the family has non-zero abundance, (ii) the goodness of fit index

[TABLE]

where $LS_{i}$ is the error associated with the exponent $\beta_{i}$, measured as the average squared deviation between the fit and the logarithm of the empirical abundance (see SI sec. S1), and (iii) the Pearson correlation coefficient $\rho_{i}$ between the logarithm of the family abundance and the logarithm of the genome size. The index $s_{i}$ puts on the same ground families with different exponents, but generally decreases as the scaling exponent increases, in accordance with the growth of fluctuations in families with higher exponents observed in ref. Grilli2014 . Hence, we decided to use it only for low exponents, where the Pearson correlation is a bad proxy of scaling. We considered families with $s_{i}>0.9$ and $o_{i}>0.6$ for exponents lower than 0.2, otherwise families with $\rho_{i}>0.4$ and $o_{i}>0.6$ reducing the dataset to 357 superfamilies and 178 clans that satisfy both requirements. As shown in Fig. S1A, $s_{i}$ and $o_{i}$ are not mutually correlated across the genomes, implying that the two requirements are in fact independent, the same is valid for $\rho_{i}$ and $o_{i}$, see Fig. S1B. We verified that the removed families with the procedure described above do not influence the scaling of the category. Supplementary Fig. S2 reports the exponent of the category scaling before the thresholding (where all the families are considered) and after (where the domains belonging to the removed families are not considered in the category scaling), showing that the values are consistent for all the categories studied.

For each family within a given functional category, we defined a “heterogeneity score” $h_{i}$ as follows:

[TABLE]

where $\beta_{i}$ and $\beta_{c}$ are,respectively, the scaling exponents of family $i$ and functional category $c$. The heterogeneity measure for each functional category was defined as the average of the per-family heterogeneity scores $h_{i}$:

[TABLE]

where $F_{c}$ is the number of families in category $c$.

The significance of the values found with this formula was assessed against a null model assuming that the total abundance of a category is distributed randomly across the associated families. The average abundance (i.e. the fraction of domains belonging to a family averaged over genomes) and occurrence (fraction of genomes where the family is present) of each family are both conserved (note that these two properties are uncorrelated in the data, hence we chose to conserve both in the null model, assuming that they are independent, see Supplementary Fig. S3).

Given a genome $g$ with $n^{g}_{c}$ elements in the functional category $c$, divided into $F^{g}_{c}$ associated families, we redistributed the $n^{g}_{c}$ members among the $F^{g}_{c}$ sets conserving the average relative abundance of each family (see SI sec. S2). A member of family $i$ belonging to category $c$ was therefore added with probability

[TABLE]

The resulting set of $F_{c}$ artificially built evolutionary categories constrains the occurrence pattern and the average abundance of the original ones. Scaling exponents for families in the null model are extracted with the procedure described above. Only functional categories containing domains from more than $10$ distinct families were compared to the null model. All procedures were implemented as custom Python 2.7 scripts.

III RESULTS

Families have individual scaling exponents, reflected by

family-specific scaling laws

We started by addressing the question of whether individual families show scaling laws, and thus can be associated to specific scaling exponents. In order to do so, we isolated domains belonging to the same family across the sample of 1112 species-representative bacterial proteomes and plotted their abundance against the total number of domains in the corresponding proteome.

When the abundance is sufficiently high to overcome sampling problems, most families show a clearly identifiable individual scaling when plotted as a function of genome size. As an example, Fig. 1 shows the scaling of a set of chosen families in four selected functional categories. Additionally, some low-abundance families that occur in all genomes with a very consistent number of copies show definite scaling with exponents close to zero Grilli2014 , being clearly constant with size, with little or no fluctuations.

Given that functional categories follow specific scaling laws, likely related to function-specific evolutionary trends vanNimwegen03 ; Molina08 , there remain different open possibilities for the behavior of the evolutionary families composing the functional categories. One simple scenario is that family scalings are family-specific, thus validating the existence of family evolutionary expansion rates that are quantitatively different to the one of their functional category. In the opposite extreme scenario the scaling is only function-specific, and individual families performing similar functions are interchangeable. If this were the case, family diversity in scaling exponent would be only due to sampling effects, and the null model would fully reproduce the diversity in family scaling observed in empirical data. To address this question, we randomized the families within a category conserving their occurrence patterns and the category average abundance. The randomized families always show very similar scaling as the one of the corresponding category (see Supplementary Fig. S4). Hence, this analysis strongly supports the existence of family-specific scaling exponents that do not simply descend from the category scaling.

Fig. 1 shows that the presence of “outlier families” is common among functional categories. In most categories, we found families where the deviations from the category exponents is clear, beyond the uncertainty due to the errors from the fits. Fig. 1 shows some examples where in each of the shown categories $\beta_{i}$ may be higher, lower or comparable to $\beta_{c}$. A table containing all the family and category exponents is available as supplementary information (Supplementary Table S1 and Table S2).

Finally, we considered the correlation of family scaling exponents with relevant biological and evolutionary parameters such as foldability (quantified by size-corrected contact order, SMCO Debes2013 ), the diversity of EC-numbers associated with families (quantifying the functional plasticity of a given family), selective pressure (quantified by the ratio of nonsynonimous to synonymous $K_{a}/K_{s}$ substitution rates Ndhlovu2015 ) and overall family abundance. The results are summarized in Supplementary Table S2. Foldability and $K_{a}/K_{s}$ appear to have little correlation with scaling exponents. Instead, we found a significant positive correlation of exponents with family abundance, and both quantities are correlated with diversity of EC-numbers in metabolic families. This suggests that, at least for metabolism, functional properties of a fold play a role in family scaling, and that beyond metabolism, abundance and scaling are, on average, not unrelated.

The heterogenetiy in scaling exponents is

function-specific.

The analyses presented above support the hypothesis that functional categories contain families with specific scaling exponents. Indeed, the scaling exponents $\beta_{i}$ of the families can be significantly different from the category exponent $\beta_{c}$ with deviations that are much larger than predicted by randomizing the categories according to the null model (see Supplementary Fig. S4).

In order to quantify this “scaling heterogeneity” of functional categories, we computed for each family $i$ the distance between its scaling exponent $\beta_{i}$ and the category exponent $\beta_{c}$ (see Methods). We defined an index $H_{c}$ quantifying the heterogeneity of the scaling of the families within a category by averaging this distance over the families associated to a given category $c$.

Figure 2A shows the relation between the heterogeneity $H_{c}$ and the category exponent $\beta_{c}$. Interestingly, these two quantities are correlated, with categories with larger values of $\beta_{c}$ being more heterogeneous. This result can be intuitively rationalized in terms of the degrees of freedom imposed by the category exponent to the scaling of single families. Categories with small exponents are incompatible with extremely large fluctuations of family exponents, while categories with larger exponents can contain families with small $\beta_{i}$. Indeed, this trend of heterogeneity with exponents is also observed in the null model, where the hetereogenity of null categories is much smaller than empirical ones, since all families tend to take the exponent category (Supplementary Fig. S4).

Figure 2B allows a direct comparison of the heterogeneity of different categories by subtracting the mean trend. It is noteworthy that the Signal Transduction functional category, which also has clear superlinear scaling, has much lower heterogeneity than DNA-binding/transcription factors. Among the categories with linear scaling, Transferases is one of the least heterogeneous ones, while the categories Protein Modification and Ion metabolism and Transport show a large variability in the exponents of the associated families. For Protein Modification, this signal is essentially due to the Gro-ES superfamily and to the HFSP90 ATP-ase domain, which have a clear superlinear scaling, while other chaperone families, such as FKBP, HSP20-like and J-domain are clearly sublinear with exponents close to zero. Interestingly, the Gro-EL domains, functionally associated to the Gro-EL, are part of this second class (exponent close to 0.2), showing very different abundance scaling to the Gro-EL partner domains. Conversely, the category Ion Metabolism and Transport is divided equally into linearly scaling (e.g., Ferritin-like Iron homeostasis domains) and markedly sublinear families, such as SUF (sulphur assimilation) / NIF (nitrogen fixation) domains. On the other hand, categories with small values of heterogeneity are made of families with exponents close to the one of the category, as shown in Table 1 in the case of, e.g., Transferases.

Determinants of the scaling exponent of a functional category

We have shown that scaling exponents of individual families may correspond to a vsriable extent to the exponent of the corresponding functional category. However, since categories are groups of families, the scaling of the former cannot be independent of the scaling of the latter. This section explores systematically the connection between the two. As detailed below, we find that in some cases the scaling exponent of functional categories is determined by few outlier families, while in other cases most of the families within a category contribute to the category scaling exponent.

While many families have a clear power-law scaling, functional categories may contain many low-abundance families with unclear scaling properties. When considered individually, these families do not contribute much to the total number of domains of a category, but their joined effect on the scaling of the category could be potentially important. Supplementary Fig. S5 shows that the sum of these low-abundance families does not suffer from sampling problems and shows a clear scaling. Interestingly, the scaling exponents for these sums once again does not necessarily coincide with the category exponents.

Figure 3A illustrates the systematic procedure that we used in order to understand how the scaling of categories emerges from the scaling of the associated families. Families were ranked by total abundance across all genomes (from the most to the least abundant) and removed one by one from the category. At each removal step in this procedure, both the scaling exponent of the removed family and the exponent of the remainder of the category are considered. In other words, the $i$-th step evaluates the exponent of the $i$-th ranking family (in order of overall abundance) and of the set of families obtained by removing the $i$ top-ranking families (with highest abundance) from the category. The resulting exponents quantify the contribution of each family to the global category scaling, as well as the collective contribution of all the families with increasingly lower overall abundance.

The results (Fig. 3B and Supplementary Fig. S6 and S7), show how the heterogeneity features described above are related to family abundance. The collective behavior of low-abundance families may deviate sensibly from that of the functional category and families follow scaling laws that sensibly deviate from the one of the corresponding functional categories. One notable example of this are Transcription-Factor DNA-binding domains. If the abundance of the outliers families is large enough in terms of the fraction of domains in the functional category, they might be responsible for determining the scaling of the entire category, as it happens in the case of DNA-binding (which is more extensively discussed in the following section).

Overall, one can distinguish between two main behaviors, either a category scaling is driven by a low number of highly populated “outlier” families (e.g. DNA binding and Protein Modification in Fig. 3B), or the category scaling is coherent, and robust to family subtraction (e.g. Transferases and Translation in Fig. 3B). While the first behavior appears to be more common for functional categories with higher scaling exponent, there are some exceptions. Notably, the scaling of strongly super-linear categories is not always driven by a few families. For example, the functional category Signal Transduction has an exponent $\beta_{c}=1.7$, which remains stable after the removal of the largest families (Supplementary Fig. S5 and S7). Both behaviors are clearly visible for intermediate exponents (in order to appreciate this, compare the Transferases and Protein Modification categories in Fig. 3B).

Super-linear scaling of transcription factors is

determined by the behavior of a few specific highly populated families.

We considered, in particular, the case of DNA-binding / transcription factors charoensawan2010genomic , which are known to exhibit peculiar scaling in bacteria Ranea2004 ; Maslov2009 . The abundance of domains in this functional category increases superlinearly (almost quadratically) with the total domain counts vanNimwegen03 ; Pang2011 ; Grilli2014 . As shown in the first row of Fig. 1B, not all the families in this functional category display a superlinear scaling charoensawan2010genomic , and the collective scaling of the low-abundance families with genome size is much slower (see Fig. 3 and Supplementary Fig. S5). Fig. 3B shows that only the most 5-6 abundant families display a super-linear scaling ($\beta_{i}>1$). These are Winged helix DNA-binding domains (34.8% of abundance), Homeodomain-like (23.3 %) lambda repressor-like DNA-binding domains (9.5%) bipartite Response regulators (7.7%) Periplasmic binding protein-like (6.2%), and FadR-like (2.4%). The remaining 16.1% of the DNA-binding regulatory domains follows a clear sublinear scaling with genome size (exponent 0.7, see Supplementary Fig. S5).

Grouping families with similar scaling exponents shows

known associations with biological function and reveals new ones.

The above analyses show that the range of scaling exponents of families within the same functional categories is generally wide and that the scaling behavior of some families sensibly deviates from their category. At the same time, functional categories show clear characteristic scaling laws, with well-defined exponents $\beta_{c}$ molina2009scaling . We, therefore, asked to what extent a range of family scaling exponents $\beta_{i}$ is peculiar to a functional category and how this compares to the category exponent $\beta_{c}$. To this end, we grouped families based on their scaling exponents. We then used those groups to test how much specific range of exponents define specific functions by an enrichment test of functional annotations.

Table 1 shows that in most cases functional categories are over-represented in the exponent range where their scaling exponents $\beta_{c}$ is found. This confirms and puts in a wider perspective the previously reported strong association between abundance scaling with size and functional annotation. As can be expected from previous results, the functional category Protein Modification is an exception: this category is under-represented in the linear region even though its category exponent is $\sim 1.06$, since it contains two strongly superlinear families and a bulk of families with sublinear scaling. This strong heterogeneity in scaling exponents is also visible in Fig. 3B.

The results of this analysis are not sensitive to the chosen intervals for the scaling exponents. In order to show this, we performed a more systematic enrichment analysis, using sliding windows of exponents of width 0.4, and step 0.1, and plotting the Z-score for the enrichment as a function of the representative family exponent for each window (Supplementary Fig. S9). The maxima of this plot define a representative exponent for each functional category, and can be compared to the exponent $\beta_{c}$ measured directly from the plot of category abundance vs genome size (see Supplementary Figure S10). Interestingly, this analysis also shows that in many cases a single functional category is enriched for multiple groups of families with well-defined exponents, as in the case of the Protein Modification category. The cases of Ion Metabolism and Transport (already discussed), Coenzyme Metabolism and Transport, Redox also shows clear indications of enrichment for two or more exponent groups. For the category Coenzyme Metabolism and Transport this is due to the presence of a single abundant family with scaling exponent close to 2, the acyl-CoA dehydrogenase NM domain-like, whose functional annotation is still not well defined. In the case of Redox, the most abundant families (Thioredoxin-like, 4Fe-4S ferredoxins, Metallo-hydrolase/Oxydoreductase) scale linearly, but there is a wide range of families with exponents between 0.5 and 1, and once again two fairly abundant outlier families with superlinear scaling (Glyoxalase/Bleomycin resistance protein/Dioxygenase, and ALDH-like), both with a fairly wide range of functional annotations.

IV DISCUSSION AND CONCLUSIONS

Our results gather a critical mass of evidence in the direction of family-specific expansion rules for the families of protein domains found in a genome. Although previous work had focused on individual transcription factor families charoensawan2010genomic , finding in some cases some definite scaling, no attempts were made to address this question systematically. The scaling laws for domain families appear to be very robust, despite of the limited sampling of families compared to functional annotations (which are super-aggregates of families and hence have by definition higher abundance). In particular, the results are consistent between the different classifications of families we tested (SUPERFAMILY and PFAM, see SI, sec. S4).

Overall, our results indicate that scaling laws are measurable at the family level, and, given the heterogeneous scaling of families with the same functional annotations, families are likely a more reliable description level for the scaling laws than functional annotations. The interpretation of these scaling laws is related to the evolutionary dynamics of family expansion by horizontal transfer or gene duplication, and gene loss Koonin2011 ; vanNimwegen03 ; Grassi2012 . Scaling exponents are seen as “evolutionary potentials” Molina08 , is based on a model of function-specific (multiplicative) family expansion rates. Assuming this interpretation, then our result that these rates may be different for different domain families having the same functional annotation may seem puzzling. Clearly, selective pressure can only act at the functional level, and if two folds were functionally identical, there should be reasonably no advantage selecting one with respect to the other, and doing so at different specific rates. For example, a transcription factor using one fold to bind DNA rather than another one should be indistinguishable from one using a different fold, provided binding specificity and regulatory action are the same.

In view of these considerations, we believe that our findings support a more complex scenario for the interplay between domain families and their functions. Specifically, we put forward two complementary rationalizations, both of which are probably in part verified in the data. The first is that functional annotations group together different domains whose abundance is linked in different ways to genome size because of their different biochemical and biological functional roles. Such differences may range from slight biochemical specificities of different folds to plain misannotations. This is possible, e.g., with enzymes, where the biochemical range of two different folds is generally different. This obervation might be related to the positive correlation we found between the number of EC numbers corresponding to a metabolic domain and its scaling exponent. However, such interpretation might be less likely applicable to, e.g., transcription factor DNA-binding domains, where functional annotation is fairly straightforward Babu03 , but different scaling behaviors with genome size are nevertheless found.

The second potential explanation assumes the point of view where scaling laws are the result of functional interdependency between different domain families Maslov2009 ; Grassi2012a , then correlated fluctuations around the mean of family pairs should carry memory of such dependency structures Pang2013 . More in detail, there may be specific dependencies connecting the relative proportions of domains with both different and equal functional annotations that are present in the same genome, which might determine the family-specific behavior Grilli2012 . While further analysis is required to elucidate these trends, we believe that gaining knowledge on functional dependencies would be an important step to understand the functional design principles of genomes.

Of notable importance is the case of the superlinear scaling of transcription factors, which has created notable debate in the past Ranea2005 ; Maslov2009 . For the first time, we look into how this trend is subdivided between the different DBDs Babu03 . Our analysis indicates that the superlinear scaling is driven by the few most abundant superfamilies (mostly winged-helix, homeodomain, lambda repressor). However, the remaining 10-20% of the functional category gives a clear sublinear scaling with genome size, which emerges beyond any sampling problems. We speculate that these other regulatory DNA-binding domains may be functionally different or behave differently over evolutionary time scales. Hence, the scaling of transcription factors with size in bacteria is driven by a small set of domain families with scaling exponent close to two, which take up most of the abundance, but does not appear to be peculiar of all transcription factors. A “toolbox” model considering the role of transcription factors as regulator of metabolic pathways and the finite universe of metabolic reactions Maslov2009 ; Pang2011 predicts scaling exponents close to two for transcription factor families. According to our results, such model should be applicable to the leading TF families. Interestingly, the heterogeneity in the behavior in transcription factor DNA-binding domains is much higher than that of the other notable superlinear functional category, signal transduction, where removal of the leading families does not significantly affect the observed scaling of abundance with genome size. Given the clarity and uniformity of the scaling exponent, we speculate that possibly a toolbox-like model may be applicable to understand the overall scaling of this category.

Other categories clearly contain multiple sets of families with coherent exponents or single outlier families. In some cases, two main groups of families with different scaling behavior clearly emerge, and higher observed scaling exponents may be related to a wider range of functional annotations. We propose that such easily detectable trends can be used to revise and refine functional annotations of protein domains. Such functional annotations are currently largely curated by humans, and based on subjective and/or biased criteria. The analysis of family scaling gives an additional objective test to define the coherence of the families that are annotated under the same function. While yet-to-be-developed automated inference methods based on our observations could serve this purpose, the quantitative scores defined here already provide useful information. The heterogeneity of a functional category is an indication of how likely that group of domain families follows a coherent expansion rate over evolution. The enrichment scores for sets of families with a given range of scaling exponent helps to pinpoint the sets of families within the functional category that expand coherently with genome size.

V ACKNOWLEDGEMENTS

We thank Erik van Nimwegen, Madan Babu, Otto Cordero and Purushottam Dixit for helpful discussions.

Supplementary Notes

S1 Filters for reliable family scaling

Many domain families are found only in a few genomes and/or in very few copies. For this reason, they might not show clear scaling properties. We excluded such families from the analysis with some filtering criteria. In order to filter out these families, we used three independent parameters. The first one is the occurrence of a family (i.e., the fraction of genomes where it is present, and the number of points available for the fit), the second one is the Pearson correlation of its abundance with genome size and the third one is the goodness of fit to a power law.

The occurence of the family $i$ is defined as

[TABLE]

while $N_{G}^{(i)}$ is the number of genomes in which family $i$ is present and $N_{G}$ is the total number of genomes in the sample.

The Pearson correlation coefficient $\rho_{i}$ between the logarithm of the family abundance and the logarithm of the genome size quantifies the existence of a relation between family abundance and genome size. It should be noted that this quantity differs from $o_{i}$: these two properties are not correlated in the data (Fig. S1B), which indicates that filters on their value should be applied independently.

Considering families with clear but shallow scaling or constant abundance across genomes, the Pearson correlation index gives values close to zero, or slighlty negative, therefore another parameter is required to assess the accuracy of the fit results. For each family whose estimated scaling exponent is lower than 0.2, we defined the quantity $LS_{i}$ as

[TABLE]

where $y_{i}$ is the logarithm of the empirical abundance of family $i$ and $y_{fit,i}$ is the abundance calculated with the fit parameters, i.e.,

[TABLE]

The goodness of fit index $s_{i}$ was defined as

[TABLE]

so that the values of $s_{i}$ close to 1 correspond to the minimum value of the average squared deviation between the fit and the empirical values. The goodness of fit index $s_{i}$ is independent from the occurrence $o_{i}$ as shown in Fig. S1A.

We considered only the families with $o_{i}>0.6$. If the fitted scaling exponent is higher than $0.2$ then we excluded the families with $\rho_{i}<0.4$, while, for exponents lower than $0.2$, only families with $s_{i}>0.9$ were taken into account. After this thresholding, we removed 1179 families out of 1536. While the fraction of such small and sparse families is large, we verified that they do not contribute significantly to the category scaling (see Figure S2).

S2 Formulation of the null model for family-specific scaling

The null model assumes that the total abundance of a category is distributed randomly across the families belonging to it. Both the average relative abundance (i.e. fraction of domains belonging to a family averaged over genomes) and occurrence (fraction of genomes where the family is present) of each family are conserved (they are in fact two independent properties, see Fig. S3).

The null model is based on the following ingredients.

• •

The number of domains belonging to a category $c$ in genome $g$, $n_{i}^{g}$, is conserved.

• •

For each genome, domains are not assigned to families that are not present in that genome.

• •

The average frequency for each family with respect to the category is conserved.

The average frequency $f_{c}(i)$ of the family $i$ with respect to the category $c$ is defined as

[TABLE]

where the family index $i$ belongs to the set in category $c$ and the sum over $g$ is carried over all the genomes, while $N_{G}^{(i)}$ is the number of genomes in which family $i$ is present, $n_{i}^{g}$ is the abundance of family $i$ in genome $g$ and $n^{g}_{c}=\sum_{i\in c}n_{i}^{g}$ is the abundance of category $c$ in genome $g$.

Given a genome $g$, each realization of the null model redistributes randomly the $n^{g}_{c}$ domains of the functional category $c$ arranged in the $F^{g}_{c}$ families belonging to category $c$ in genome $g$. Each one of the $n^{g}_{c}$ domains is assigned to family $i$ with probability

[TABLE]

S3 Enrichment analysis for functional categories and families

with similar scaling exponent

All families passing the filters described in section S1 were divided into three groups based on the values of their exponent $\beta_{i}$:

• •

sub-linearly scaling families, $\beta_{i}\leq 0.6$

• •

linearly scaling families, $0.6<\beta_{i}<1.4$

• •

super-linearly scaling families, $\beta_{c}\geq 1.4$

We used hypergeometric tests to asses over- or under-representation of functional categories in these family groups. Given that $F_{c}$ is the number of families that belong to the category $c$, $F_{bin}$ is the number of families in either of the three groups defined above and $F_{tot}$ is the total number of families involved in this analysis, the mean and the variance of the hypergeometric distribution are:

[TABLE]

[TABLE]

For each combination of functional category and family group, the quantity $x_{bin,c}$ is the number of families that lie in the intersection of category $c$ with family group $bin$. The functional category $c$ is under-represented in the group $bin$ if $Z_{bin,c}<-1.96$, over-represented if $Z_{bin,c}>-1.96$, where $Z_{bin,c}$ is the Z-score:

[TABLE]

The resulting intersection values $x_{bin,c}$ and the significant Z-scores are reported in Table 1 of the main text.

In order to prove that the results are independent from the chosen interval of the exponents, we substituted the three groups with sliding intervals of amplitude 0.4 and step 0.1 and repeated the same process. Only intervals with more than 10 families are considered. Fig. S9 shows that the results are consistent with the previous analysis. Fig. S10 shows how the exponent corresponding to the maximum Z-score differs in some cases from the category exponent.

S4 Main results hold also for PFAM clans

We chose PFAM clans as an alternative database to test our results. PFAM clans were annotated on the same scheme of 50 functional categories used for superfamilies, using the mapping of clans into superfamilies available from the PFAM website http://pfam.xfam.org/clan/browse#numbers finn2006pfam . The scaling laws for functional categories are recovered also for clans (Fig. S11 and Table LABEL:tab:TableS1cat) and are consistent with previous results vanNimwegen03 ; Molina08 ; MOLINA09 ; cordero2009regulome ; Grilli2012 ; charoensawan2010genomic .

The following main results were recovered for Pfam clans.

• •

For each clan, the abundance across genomes scaleS as a power law of the genome size. Equally to SCOP superfamilies, Pfam clans have individual scaling exponents that may or may not follow the one of the associated functional category (Table. LABEL:tab:TableS2clan). The fitting method and threshold values are the same used for supefamilies (sec. S1). 178 clans out of 446 passed the filters and were employed for further analysis.

• •

The heterogeneity (average of the distance between the category exponent and the clan exponent), positively correlates with the category exponent (Fig. S12). Functional categories with superlinear scaling tend to be more heterogeneous and, as found for superfamilies, the functional category Signal Transduction is less heterogeneous than DNA-binding, although having the largest exponent. Unlike the case of superfamilies, Protein Modification does not have high heterogeneity score, but the difference in scaling between the (strongly superlinear) outlier family Gro-ES and the remaining ones is observed. For clans, the scaling of Protein Modification is once again strongly biased by the clan “GroES-like superfamily” (20% of the total domains).

• •

Either few or most of the clans determine the scaling exponent of the functional category they belong to. Figure S8 is coherent with what observed for superfamilies, in particular the functional category of DNA-binding is dominated by one clan (the “Helix-turn-helix” clan) that accounts for 83% of the total domains. As for superfamilies, Signal Trasduction is robust to the progressive removal of families confirming that the presence of dominant clans is not related to the superlinear scaling of the category.

• •

Grouping clans with similar scaling exponents recovers known associations between the category exponent and the biological function (Fig. S13).

S5 Correlation between family scaling and EC numbers

The Enzyme Commission (EC) number is a classification scheme for enzyme-catalyzed chemical reactions. It is built as a four-levels tree where the top nodes are six main groups of reactions, namely Oxidoreductases, Transferases, Hydrolases, Lyases, Isomerases and Ligases. We used the mapping between Superfamilies and EC terms Gough01 , to investigate the correlation between the Superfamily scaling and the number of different reactions in which the family is involved. This quantity is the count of distinct EC numbers corresponding to the finest level of the EC classification. This number shows a positive correlation (Spearman 0.74) with the scaling exponent of metabolic family (see Table S2). The diversity of EC numbers in metabolic families is also correlated with the mean total abundance of a family, since family abundance and scaling exponent are also correlated (Spearman 0.72).