I do not know who first formulated consciously the idea that there should exist a good formal definition of what is a “family of 
-functions”, or “family of cusp forms”. Such ideas were obviously present in the late 1990’s when it was discovered that there were links between special values of -functions and statistics of random matrices in families of compact classical groups, and when this was given a rigorous expression in the case of -functions of algebraic origin over finite fields by Katz and Sarnak. Specific families of various kinds and their fundamental properties were also very clearly present in the work of Iwaniec and Sarnak on non-vanishing of central values of Hecke -functions of modular forms around the same time. However, these families were one of these mathematical objects which people are happy to recognize when they see them, but which seemed difficult to define beforehand in the right generality.
For some reason, this is a question I’ve been thinking about every once in a while since that “old” time. A paper of Conrey-Farmer-Keating-Rubinstein-Snaith proposed a good heuristic proposal of what families of -functions should be, in order to be able to make robust conjectures concerning averages of products/ratios of shifted values of these -functions. However (to my mind) this proposal did not really give any tool to attempt a general theory, and it certainly was not a Bourbakist-style definition. Around 2006, when I was finishing my (almost award-winning) book on the large sieve, I noticed a clear analogy between certain problems about families of cusp forms and the formalism I was using for sieve. This convinced me that the starting point for any good notion of “family of -functions” should be the underlying cusp forms, or rather automorphic representations when dealing with the general case, and that the first principle must be: if some collection of automorphic representations (say over the rationals) deserves the name of “family”, it must be the case that, for any prime number 
, the family of -components must also be well-behaved. As it turns out, this local behavior can be given a clear meaning, and this provides a “test” for any attempted definition. In fact, for classical modular forms, this good local behavior amounts to the question of equidistribution of eigenvalues of Hecke operators, which had been treated (in various generality) by Bruggeman (a bit implicitly), Sarnak, Serre, Conrey-Duke-Farmer and Royer (and a few others). This seemed a good indication that this was the right track. (To be precise, these results only handle unramified primes; especially when the local components are or can be supercuspidals, it is not so easy to get a handle on the underlying asymptotic distribution properties…)
I happened to talk about this in 2009 with P. Sarnak (during the Verbania conference) and he told me that he had also formulated a definition of families, which was sufficient to predict a “symmetry type” for any associated family of -functions (built in the Langlands style using arbitrary representations of the -group). This was explained in an unpublished letter he had written a few months ago, which is now available on his new “blog”.
We discussed the issue further in late 2009 when I was in Princeton. Later, the philosophy of “local spectral equidistribution” (as I had decided to call the idea behind the analogy with sieve) was the background of the papers I wrote with A. Saha and J. Tsimerman last year, where some special families of Siegel modular forms were considered. But it is only recently (i.e., during the blessed summer time of mellow fruitful-writingness) that I found time to attempt to write a consistent and somewhat detailed account of this point of view on families of cusp forms. This account can be found here. It is very informal, but hopefully readable. The summary of those papers with Saha and Tsimerman (in Section 13) might interest those readers who know — and like! — Böcherer’s fascinating conjecture concerning the arithmetic nature of Fourier coefficients of Siegel cuspidal eigenforms of genus 2, and their expected relation with special -values of spinor -functions.