No one single person came up with the GEE formulation, but if we had to pick the most well known reference, it would be Liang Zeger (1986). They penned the most widely used terminology and notation for these models.
You're on the right track to connect the GEE to the quasilikelihood. So Wedderburn 1974 deserves some credit. About the same time, several highly generalized theories were being codified. McCullogh and Nelder (1983) released the first edition of Generalized Linear Models which, in part, emphasized the GLM as an abstract, method of moments type estimator which was derived using a link function and a variance function and, like quasi-likelihood, may not correspond to a proper set of likelihood equations, but has well defined limiting behaviors for inference and estimation.
There were several completely unrelated efforts to derive robust variance estimation which was unified under Liang, Zeger's GEE framework. Eicker (1963, 1967), Huber (1967), and White (1980) explored the "Heteroscedasiticity-consistent" variance estimator (HC). As pointed out by Boos and Stefanski (2013), Freedman (2006), and others, sandwich variance estimation very closely resembles model-based inference when assumptions are met. This is because in GLMs as it is for GEE, the maximum likelihood technique requires a very specific relationship between the link and the variance for exponential families.
As a motivating example, consider modeling a binomial response with logistic regression. The logistic link function is: $ g(x) = \log x / (1-x)$ and its inverse is the sigmoid which expresses the probability of response as a function of the linear predictor. $g^{-1}(x) = \exp(x) / (1-\exp(x)) = \mu$. Taking the derivative of the inverse link with respect to $x$ gives $\frac{\partial}{\partial x} g^{-1}(x) = g^{-1}(x)(1-g^{-1}(x)) = \mu(1-\mu)$ which is readily recognized as the mean-variance relationship in a binomial random variable. For this reason, the general GEE framework applied to this logistic regression scenario:
$$ \mathcal{U}(\beta) = D^TV^{-1}(y - g^{-1}(X^T\beta))$$
would yield the usual likelihood equations for a logistic regression.
References
Eicker, F. (1963), Asymptotic normality and consistency of the least squares estimator for families of linear regressions. Annals of Mathematical Statistics, 34(2), 447–456.
Eicker, F. (1967), Limit theorems for regressions with unequal and dependent errors. Pp. 59–82 in Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability. Berkeley: University of California Press.
Huber, P.J. (1967), The behavior of maximum likelihood estimates under non-standard conditions. Pp. 221–233 in Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability. Berkeley: University of California Press.
Liang, K.-L., & Zeger, S. L. (1986). Longitudinal data analysis using generalized linear models. Biometrika, 73(1), 13–22.
Wedderburn, R. W. M. (1974). Quasi-Likelihood Functions, Generalized Linear Models, and the Gauss-Newton Method. Biometrika, 61(3), 439–447.
White, H. (1980), A heteroskedastic-consistent covariance matrix estimator and a direct test of heteroskedasticity. Econometrica, 48(4), 817–838.
