The determinants of binding specificities of peptide recognition domains and their evolution remain important problems in molecular systems biology. Here, we present a new methodology to analyze the coevolution between a domain and its ligands by combining high-throughput phage display with deep sequencing. First, from a library of PDZ domains with diversity introduced at ten positions in the binding site, we evolved domains for binding to 15 distinct peptide ligands. Interestingly, for a given peptide many different functional domains emerged, which exhibited only limited sequence homology, showing that many different binding sites can recognize a given peptide. Subsequently, we used peptide-phage libraries and deep sequencing to map the specificity profiles of these evolved domains at high resolution, and we found that the domains recognize their cognate peptides with high affinity but low specificity. Our analysis reveals two aspects of evolution of new binding specificities. First, we were able to identify some common features amongst domains raised against a common peptide. Second, our analysis suggests that cooperative interactions between multiple binding site residues lead to a diversity of binding profiles with considerable plasticity. The details of intramolecular cooperativity remain to be elucidated, but nonetheless, we have established a general methodology that can be used to explore protein evolution in a systematic yet rapid manner.