BioSystems 82 (2005) 226–234

Detecting correlation between sequence and expression divergences in a comparative analysis of human serpin genes Zuofeng Li a , Qi Liu b , Mangen Song c , Ying Zheng c , Peng Nan a,d , Ying Cao e , Guoqiang Chen c , Yixue Li b,d , Yang Zhong a,d,∗ a School of Life Sciences, Fudan University, Shanghai 200433, China Bioinformation Center, Shanghai Institute for Biological Sciences, Chinese Academy of Sciences, Shanghai 200031, China Department of Pathophysiology, Shanghai Institute of Hematology, Rui-Jin Hospital, Shanghai Second Medical University, and Health Science Center, Shanghai Institute for Biological Sciences, Chinese Academy of Sciences, Shanghai 200025, China d Shanghai Center for Bioinformation Technology, Shanghai 201203, China e The Institute of Statistical Mathematics, 4-6-7 Minami-Azabu, Minato-ku, Tokyo 106 and Department of Biosystems Science, Graduate University for Advanced Studies, Hayama 240, Japan b

c

Received 18 June 2004; received in revised form 14 June 2005; accepted 9 July 2005

Abstract Physiological functions and characteristic structures of the serpin gene superfamily have been studied extensively, yet the evolution of the serpin genes remains unclear. Gene duplication in this superfamily may shed light on this issue. Two models are used to predict the preservation of duplicated genes: the classical model and the duplication–degeneration–complementation (DDC) model. In this study, we analyzed the phylogenetic relationships of 33 human serpin genes and the expression data of some members of the serpin superfamily from a DNA microarray of human leukemia U937 cells with stably inducible expression of the leukemia-related AML1-ETO gene. We then determined the utility of the DDC model by mapping serpin superfamily expression data to the phylogenetic tree. The correlation between sequence and expression divergences as measured by the Pearson correlation coefficient indicated that human serpin genes evolved under the DDC model. Our study provides a new strategy for comparative analysis of gene sequences and microarray data. © 2005 Elsevier Ireland Ltd. All rights reserved. Keywords: Serpin; Gene duplication; Microarray; Sequence divergence; Expression divergence

1. Introduction

∗ Corresponding author. Tel.: +86 21 55664436; fax: +86 21 65642468. E-mail address: [email protected] (Y. Zhong).

The serpin superfamily of serine protease inhibitors is widely distributed from viruses to vertebrates. They are involved in many biological processes such as blood coagulation, fibrinolysis, synaptic plasticity, apoptosis,

0303-2647/$ – see front matter © 2005 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.biosystems.2005.07.004

Z. Li et al. / BioSystems 82 (2005) 226–234

development, inflammation, and tumor development (van Gent et al., 2003). The structural characteristics common to the serpin superfamily are the reactive centerloop (RCL), with P1 and P1 residues acting as “bait” for cognate proteases, where the P1 residues control binding specificity. After RCL binding by the protease, the bond between the two residues (i.e., P1 and P1 ) is cleaved. Then the RCL moves to another part of the serpin and the protease is trapped in an irreversible conformational rearrangement (Hagglof et al., 2004). However, some members of the serpin superfamily seem to lose the inhibition function, and some even play a converse role. For example, although many serpin genes have been identified as having tumor suppressive activity by the inhibition of cancer cell invasion and metastasis, some serpins have been shown to be oncogenic (Ito et al., 2004). What drives the functional diversity of the serpin superfamily is still of great interest (Maxwell and Enrico, 2003). In particular, investigation of duplicated genes within the superfamily may shed some light on this issue. Two models are commonly used to predict the preservation of duplicated genes: the classical model and the duplication–degeneration–complementation (DDC) model (Force et al., 1999). The classical model suggests that the fate of most duplicated genes is to be fixed as pseudogenes. However, empirical data run contrary to this hypothesis. For example, about half of all duplicated genes of Xenopus laevis have been preserved for 30 million years (Hughes and Hughes, 1993). This percentage is larger than the percentage predicted by the classical model. The DDC model predicts that degenerate mutations in regulatory regions cause the loss and/or preservation of duplicated genes. Preserved genes partition ancestral functions rather than assist with the evolution of novel function. Spatial and temporal divergence of duplicated gene expression has long been a focus of evolutionists (Gu et al., 2002) and may be the first step in functional divergence. According to the DDC model, spatiotemporal gene expression patterns may be as similar as the predicted gradual sequence divergence after gene duplication by the model of neutral mutation. Therefore, it is possible to detect whether a gene evolved under the DDC model by comparing the sequence and expression divergences. For

227

example, the coupled evolution of coding region and gene expression of human genes has been detected through a correlation analysis of the sequence divergence and the expression pattern (Makova and Li, 2003). In this paper, the correlation between sequence and expression divergences in human serpins is investigated using the temporal microarray data of U937 cells containing the transcription factor AML-ETO. We reconstructed the phylogeny of the human serpins using protein-coding regions and then mapped gene expression data to this phylogeny to determine any correlations. Correlations between the rate of synonymous substitution and temporal divergence of the expression data were determined to test the utility of gene duplication models.

2. Materials and methods 2.1. Sequences and sequence divergence The nucleotide sequences used in this study were retrieved from GenBank using Batch Entrez, and were stored in one file. Formatdb (Altschul et al., 1997) was used to format the file to be a searchable database for BLAST (Altschul et al., 1997). The coding nucleotide sequences of 34 members of the serpin superfamily (including 33 human serpin genes and the serpin gene from Hordeum vulgare as an outgroup for further phylogenetic analysis) were retrieved from GenBank (Table 1). Clustal X was used for multiple sequence alignment (Thompson et al., 1997) with pairwise sequence identifies 36%. The 33 members of the human serpin superfamily segregate evolutionarily into nine (A-I) distinct clades (Silverman et al., 2004) as shown in Fig. 1. The coding sequences of the four members in Clade B and the five members in Clade A were analyzed, respectively. The protein sequences were also aligned with Clustal X and the gaps were deleted from the nucleotide alignment by imposing the amino acid residues for the resulting protein. The codons of the amino acid residues were aligned and the PAML program was used to estimate synonymous (dS ) and nonsynonymous (dN ) substitution rates with standard error (SE) (Yang, 1997). dS and dN were used as proxies for sequence divergence.

228

Z. Li et al. / BioSystems 82 (2005) 226–234

Fig. 1. A maximum parsimonious phylogenetic tree of 33 human serpin genes. Bootstrap values were calculated from 100 replicates. A-I represents the nine clades of vertebrate serpin superfamily classified by Silverman et al. (2001). This tree is rooted using the Z serprin of Hordeum vulgare as an outgroup.

Z. Li et al. / BioSystems 82 (2005) 226–234

229

Table 1 Gene names, functions, and GenBank accession numbers of the serpin superfamily members used in this study Name

Gene name

Accession No.

Length (nt)

SERPINA1 SERPINA3 SERPINA4 SERPINA5 SERPINA6 SERPINA7 SERPINA8 SERPINA10 SERPINA11 SERPINB1 SERPINB2 SERPINB3 SERPINB4 SERPINB5 SERPINB6 SERPINB7 SERPINB8 SERPINB9 SERPINB10 SERPINB11 SERPINB12 SERPINB13 SERPINC1 SERPIND1 SERPINE1 SERPINE2 SERPINF1 SERPINF2 SERPING1 SERPINH1 SERPINH2 SERPINI1 SERPINI2 H. vulgare

Alpha-1 antiproteinase Alpha-1-antichymotrypsin Kallistatin Protein C inhibitor Corticosteroid-binding globulin Thyroxine-binding globulin Angiotensinogen Protein Z-dependent protease inhibitor Serpin A9a Monocyte/neutrophil elastase inhibitor Placental plasminogen activator inhibitor Squamous cell carcinoma antigen 1 Squamous cell carcinoma antigen 2 Maspin Cytoplasmic antiprotease Megsin Cytoplasmic antiproteinase 2 Cytoplasmic antiproteinase 3 Bomapin Epipin Yukopin Headpin Antithrombin III Heparin cofactor II Plasminogen activator inhibitor type 1 Glia-derived neurite-promoting factor Pigment epithelium-differentiation factor Alpha-2-plasmin inhibitor C1 inhibitor Collagen-binding protein Collagen-binding protein 2 Neuroserpin Pancpin

NM 000295 NM 001085 L19684 NM 000624 J02943 NM 000354 NM 000029 AF181467 NM 175739 NM 030666 J02685 U19556 U19557 U04313 NM 004568 AF027866 L40377 L40378 U35459 NM 080475 NM 080474 AF169949 X68793 NM 000185 NM 000602 A03911 M76979 D00174 NM 000062 X61598 D83174 NM 005025 NM 006217 X97636

1257 1302 1284 1221 1218 1248 1458 1335 1113 1140 1248 1173 1173 1128 1131 1143 1125 1131 1194 1179 1218 1176 1395 1500 1209 1194 1257 1476 1503 1254 1257 1233 1218 1203

a

A cross-annotation of Genbank.

2.2. Phylogenetic analysis A phylogenetic tree of the human serpin genes was generated with the maximum parsimony (MP) method, implemented in the PHYLIP 3.75 package (Felsenstein, 1997). The nucleotide sequence of Z serpin of Hordeum vulgare was selected as the outgroup. 2.3. Microarray data Acute monocytic leukemia cell line U937 cells with conditional AML-ETO expression, called U937AE cells, which had been generated by the ecdysone inducible system, were maintained in RPMI-1640

medium (Sigma, St. Louis, USA) supplemented with 10% heated-inactivated fetal calf serum (Gibco BRL, Gaithersburg, USA) and 400–600 ␮g/ml G418 and 100–150 ␮g/ml Zeocin in a 5% CO2 –95% air humidified atmosphere at 37 ◦ C. The medium was changed every second day. For induction of AML1/ETO expression, 5 ␮M ponasterone A (1–5 mM stock solution in ethanol) was added directly to the cell culture medium to rapidly induce AML1-ETO expression. At 2, 6, 9, 12, 24, 48, and 72 h after ponasterone A treatment, total RNAs were extracted from cells by TRIzol reagent (Gibco BRL, Gaithersburg, USA), and reverse transcriptions were performed by TaKaRa RNA PCR kit (Takara, Dalian, China) following the

230

Z. Li et al. / BioSystems 82 (2005) 226–234

manufacturer’s instructions. For every experiment, a Western blot was performed in parallel, to confirm high-level AML1/ETO protein expression. Microarrays with 12,600 unique cDNA clones, including 9436 UniGene hits and over 600 ESTs, were fabricated in-house using a Generation III spotter (Molecular Dynamics, Sunnyvale, USA). These clones were verified by sequencing. Total RNA of U937-AE cells was prepared with Trizol Reagent (Invitrogen, San Diego, USA), and quantified using an RNA LabChip (Agilent, Palo Alto, USA) before being subjected to reverse transcription labeling. Approximately 20 ␮g of each RNA sample was reversely transcribed into cDNA primed with oligo (dT) and labeled with either Cy3dCTP or Cy5-dCTP through Superscript II (Invitrogen, San Diego, USA) reactions. Hybridization was performed in the presence of 2 ␮g poly d (A), 3XSSC, 0.1% SDS, 10 ␮g human Cot-1 and equal amounts of labeled target and reference cDNA at 65 ◦ C overnight. Data acquisition was performed through a laser scanner (Axon, Foster City, USA). The resulting images were analyzed using GenePix Pro Software (Axon, Foster City, USA). The ratios of Cy3–5 were computed using three methods: the medium of ratio, the ratio of median and the regression ratio. The point where the variances among the three ratios were more than 20% was defined as an unreliable point and omitted from this study. Log2 -transformed ratios (each gene’s mRNA levels to its mRNA levels before ponasterone A was added) of gene expression were used for further analysis. 2.4. Expression divergence and mapping The Pearson correlation coefficient (r) was used to define the expression divergence between two genes. For two duplicated genes, r is initially +1 (Wagner, 2000), however, functional divergence is indicated as r approaches −1. Between +1 and −1 lower r-value indicate greater expression divergence. The SPSS 10.0 package (SPSS, Inc, Chicago, USA) was used to compute the r-value of gene pairs. A simplified phylogenetic tree was constructed by deleting the clades without corresponding expression data from the MP tree. The expression data were then mapped onto the simplified tree. In order to show the noise ratios of the experiment, the expression data of three housekeeping genes (HK genes): beta

actin (ACTB), glyceraldehyde-3-phosphate dehydrogenase (GAPD) and lactate dehydrogenase A (LDHA), were also labeled on the mapping figure (Eisenberg and Levanon, 2003). Compared to the expression pattern of the three HK genes, serpin genes determined in this study exhibit an inducible expression pattern. 2.5. Statistical analysis The Origin 6.0 (OriginLab Corporation, Northampton, USA) algorithm was used to construct the scatter plot graph, and to process linear or polynomial regression analyses. If p < 0.05, the correlation was considered statistically significant.

3. Results 3.1. Phylogenetic relationships among human serpin genes The human genome encodes for at least 35 members of the serpin superfamily. There are 11 and 13 members for Clades A and B, respectively (Silverman et al., 2004). The nomenclature of the serpin superfamily suggested by Silverman classified the vertebrate serpin into nine (A-I) distinct clades based on the number of exons and protein construction (Silverman et al., 2001). The topology of the maximum parsimony tree developed in this study (Fig. 1) differs from the topology of the tree constructed by Silverman et al. (2004) and van Gent et al. (2003). The Clades B, H, E, and I were monophyletic group. However, Clade A was divided into two groups and was shown to be a polyphyletic group. 3.2. Mapping the expression data of the serpin superfamily onto the phylogenetic tree The expression data were mapped onto the simplified tree (Fig. 2). No clear association between sequence and expression divergences can be seen from the mapping pattern. For Clade A, there is a positive correlation between evolutionary distance and expression divergences. Clade B serpins/ovalbumin family (so named after chicken albumin protein) reside primarily within the

Z. Li et al. / BioSystems 82 (2005) 226–234

231

Fig. 2. Mapping expression data onto the phylogenetic tree of the serpin superfamily. The positive histogram columns indicate up-regulated expression whereas negative histogram columns indicate gene expression. A histogram of three housekeeping genes (HK genes): ACTB, GAPD and LDHA, is labeled on the mapping figure.

232

Z. Li et al. / BioSystems 82 (2005) 226–234

cytoplasm because of a lack of signal peptides. The difference between Clades A and B is that the latter clade conserves the function of enzyme inhibition. Some members of Clade B serpins play a key role in tumor development or suppression. Clade B serpins are tandemly distributed on human chromosome 18q and 6p and their tissue expression is diverse. For example, although Serpins B1 and B6 distribute broadly in many tissues, the distribution of Serpins B10 and B13 is tissue-specific. Therefore, different tissues might have different gene expression patterns. This indicates that the evolution of the Clade B serpin after gene duplication may follow the DDC model, implying that the mutation of the cis-regulation region will cause the divergence of tissue expression and preserve the duplicated genes.

The synonymous (dS ) and non-synonymous (dN ) substitution rates of Clades B and A were measured. The values of the Pearson correlation coefficient (r) of gene expression were calculated for each gene pair. In the Clade B serpins, a negative correlation which was statistically significant was observed between dS and the expression divergence of the gene pair (r2 = 0.73, p < 0.05) (Fig. 3a). There was also a weak, albeit significant, association between dS and the expression divergence within the Clade A serpin (r2 = 0.46, p < 0.05) (Fig. 3c). dS is therefore negatively correlated with expression divergence. In other words, spatiotemporal expression similarity of the Clade A or B serpins declines gradually after duplication. Mutations in the cis-regulatory sequence results in the spatiotemporal expression divergence. Interestingly, there is a

Fig. 3. Mutation rates plotted over expression correlation. (a) Clade B serpins synonymous substitution rate (dS ), n = 6; (b) Clade B serpins non-synonymous substitution rate (dN ) n = 6; (c) Clade A serpins dS , n = 10; (d) Clade A serpins dN , n = 10.

Z. Li et al. / BioSystems 82 (2005) 226–234

polynomial relationship between dN and the expression divergence of gene pairs in Clades A or B (Fig. 3b and d).

4. Discussion In this study, a weak statistical association is observed between synonymous mutation rates and temporal gene expression patterns. This holds true in human Clades A and B. In an earlier study, a negative association between spatial expression divergence and dS was detected in humans by a statistical analysis based on genome data (Makova and Li, 2003). Moreover, our study detects a temporal association at gene family level. During the last decade, much has been learned about the evolutionary relationships of the serpin superfamily. Serpins may have evolved by duplication within or between chromosomes. After duplication, each gene obtains a specific function through subfunctionalization, like zebrafish engrailed-1 and 1b genes. They express at different regions and partition the ancestor genes without neofunctionalization (Zhang, 2003). Therefore, the most obvious interpretation of our results is that the human serpin genes evolved under the DDC model. In this study, a temporal microarray data set of U937 cells with conditional AML1-ETO was used to analyse the association between synonymous mutation rate and gene temporal expression pattern. AML1-ETO may block myeloid differentiation, prevent apoptosis and cause the cellular transformation (McNeil et al., 1999). With the expression of AML-ETO, the expression patterns of serpins in U937 cells show some characteristic pattern in tumors. For example, Serpin A5 is one of the well-studied serpin genes (also named protein C inhibitor, shown in Table 1). The lack of Serpin A5 expression in a subpopulation of high-grade tumor cells in combination with maintained protease expression may facilitate invasive growth patterns of tumors (Cao et al., 2003). In our study, the expression of Serpin A5 was down regulated (Fig. 2). Moreover, U937 cells with conditional AML1-ETO were an ecdysone inducible system. This makes the effects of the introduced AML1-ETO gene easily be monitored in a rather controlled manner and simplifies the explanation of the result. Therefore, this system was selected for the analysis.

233

The association detected in this study may be used to predict the expression patterns of the genes that were omitted from the data set due to unreliability. For example, the expression data of Serpin B5 (also named Maspin, shown in Table 1) were omitted because of the incomplete data. Since Serpin B1 and Serpin B5 are in a sister group (Fig. 1), we predict that the expression pattern between them should be similar. That is to say Serpin B5 may also be down regulated. This supposition is supported by the study of Zhang and Zhang (2002). They used the tissue microaary to determine the expression of Serpin B5 in various normal tissues and cancers. The results indicated that Serpin B5 expression was consistently down regulated during tumor progression (Zhang and Zhang, 2002). Gene expression in humans is multi-dimensional (Hughes and Friedman, 2005). The analysis of normal tissues or tumor tissues other than those considered here may reveal different expression patterns than those shown in our data set. In addition, in some cases, the genes that only express in other tissues may not in our data set. Nonetheless, it is important that the substantial change in the gene coding region is accompanied by a change in expression pattern, as shown in this study. This suggests that functional divergences of duplicated genes evolve with sequence divergences, even just after duplication. Makova and Li (2003) found that the patterns of protein sequence divergence and divergence of spatial expression may initially couple. A weak association was also detected in this study. However, further work is needed since our analysis was based on a single cell line, and it is therefore not possible to deduce general guidelines for all cells. Moreover, a statistical model is needed to make the results more reliable (Gu, 2004). Nonetheless, our study provides a new strategy for comparative analysis of gene sequence and microarray data. In other words, more information may be used to assist in the analysis of the microarry data, such as phylogenetic information.

Acknowledgements We are grateful to Yupeng Geng, Li Wang, Yidong Lei, Qingbiao Wang, Jimei Liu for comments on the earlier versions of the manuscript; to Professor Bruce Anderson (Queen’s University, Canada) for help in

234

Z. Li et al. / BioSystems 82 (2005) 226–234

revising this manuscript and to two anonymous reviewers for valuable comments. This work was partially supported by the National Key Project for Basic Research (973)(2002CB5128 and 2003CB715904). References Altschul, S.F., Madden, T.L., Schaffer, A.A., Zhang, J.H., Zhang, Z., Miller, W., Lipman, D.J., 1997. Gapped BLAST and PSIBLAST: a new generation of protein database search programs. Nucleic Acids Res. 25, 3389–3402. Cao, Y., Becker, C., Lundwall, A., Christensson, A., Gadaleanu, V., Lilja, H., Bjartell, A., 2003. Expression of protein C inhibitor (PCI) in benign and malignant prostatic tissues. Prostate 57, 196–204. Eisenberg, E., Levanon, E.Y., 2003. Human housekeeping genes are compact. Trends Genet. 19, 362–365. Felsenstein, J., 1997. PHYLIP (Phylogeny Inference Package) Version 3. 75c. University of Washington, Seattle. Force, A., Lynch, M., Pickett, F.B., Amores, A., Yan, Y.L., Postlethwait, J., 1999. Preservation of duplicate genes by complementary, degenerative mutations. Genetics 151, 1531–1545. Gu, X., 2004. Statistical framework for phylogenomic analysis of gene family expression profiles. Genetics 167, 531– 542. Gu, Z.L., Nicolae, D., Lu, H.H.S., Li, W.H., 2002. Rapid divergence in expression between duplicate genes inferred from microarray data. Trends Genet. 18, 609–613. Hagglof, P., Bergstrom, F., Wilczynska, M., Johansson, L.B.A., Ny, T., 2004. The reactive-center loop of active PAI-1 is folded close to the protein core and can be partially inserted. J. Mol. Biol. 335, 823–832. Hughes, M.K., Hughes, A.L., 1993. Evolution of duplicate genes in a tetraploid animal Xenopus laevis. Mol. Biol. Evol. 10, 1360–1369. Hughes, A.L., Friedman, R., 2005. Expression patterns of duplicate genes in the developing root in Arabidopsis thaliana. J. Mol. Evol. 60, 247–256. Ito, Y., Yoshida, H., Tomoda, C., Uruno, T., Takamura, Y., Miya, A., Kobayashi, K., Matsuzuka, F., Matsuura, N., Kuma, K., Miyauchi, A., 2004. Maspin expression is directly associated

with biological aggressiveness of thyroid carcinoma. Thyroid 14, 13–18. Makova, K.D., Li, W.H., 2003. Divergence in the spatial pattern of gene expression between human duplicate genes. Genome Res. 13, 1638–1645. McNeil, S., Zeng, C.M., Harrington, K.S., Hiebert, S., Lian, J.B., Stein, J.L., van Wijnen, A.J., Stein, G.S., 1999. The t(8;21) chromosomal translocation in acute myelogenous leukemia modifies intranuclear targeting of the AML1/CBF alpha 2 transcription factor. Proc. Natl. Acad. Sci. USA 96, 14882–14887. Maxwell, M.K., Enrico, D.C., 2003. Conserved ser residues, the shutter region, and speciation in serpin evolution. J. Biol. Chem. 278, 37810–37814. Silverman, G.A., Bird, P.I., Carrell, R.W., Church, F.C., Coughlin, P.B., Gettins, P.G.W., Irving, J.A., Lomas, D.A., Luke, C.J., Moyer, R.W., Pemberton, P.A., Remold-O’Donnell, E., Salvesen, G.S., Travis, J., Whisstock, J.C., 2001. The serpins are an expanding superfamily of structurally similar but functionally diverse proteins—evolution, mechanism of inhibition, novel functions, and a revised nomenclature. J. Biol. Chem. 276, 33293–33296. Silverman, G.A., Whisstock, J.C., Askew, D.J., Pak, S.C., Luke, C.J., Cataltepe, S., Irving, J.A., Bird, P.I., 2004. Human Clade B serpins (ov-serpins) belong to a cohort of evolutionarily dispersed intracellular proteinase inhibitor clades that protect cells from promiscuous proteolysis. Cell Mol. Life Sci. 61, 301–325. Thompson, J.D., Gibson, T.J., Plewniak, F., Jeanmougin, F., Higgins, D.G., 1997. The Clustal-X windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools. Nucleic Acids Res. 24, 4876–4882. van Gent, D., Sharp, P., Morgan, K., Noor, K., 2003. Serpins: structure, function and molecular evolution. Intl. J. Biochem. Cell Biol. 35, 1536–1547. Wagner, A., 2000. Decoupled evolution of coding region and mRNA expression patterns after gene duplication: Implications for the neutralist-selectionist debate. Proc. Natl. Acad. Sci. USA 97, 6579–6584. Yang, Z., 1997. PAML: a program package for phylogenetic analysis by maximum likelihood. CABIOS 13, 555–556. Zhang, J.Z., 2003. Evolution by gene duplication: an update. Trends Ecol. Evol. 18, 292–298. Zhang, W.G., Zhang, M., 2002. Tissue microarry analysis of maspin expression and its reverse correlation with mutant p53 in various tumors. Int. J. Oncol. 20, 1145–1150.