The history of the least squares mean, its appearance in SAS, and its interpretation is discussed in Searle, Milliken, and Speed (1979).
Some discussion of the concepts around least squares means (population marginal means) is found in Searle, Speed, and Milliken (1980). Earliest mention of the concept that they note is Damon et al (1959). They provide some other references, but I do not have access to the full article.
The initial implementations of the calculation seem to have been worked out explicitly in Harvey (1960) and some subsequent publications, including but probably not limited to Harvey (1977), Goodnight (1979), Harvey (1982), and Goodnight and Harvey (1997).
It looks like the computation routines were first developed as
LSML 76 and
LSML GP before the user-contributed
J.H. Goodnight (1979) A tutorial on the SWEEP operator. The American Statistician 33 (3): 149-159.
J.H. Goodnight and W.R. Harvey (1997) Least squares means in the fixed effects general model. SAS Technical Report R-103. SAS Institute Inc.
W.R. Harvey (1960) Least-squares analysis of data with unequal subclass numbers. USDA National Agricultural Library ARS-20-8.
Harvey, W.R (1977) User's guide for LSML 76. Mixed model least-squares and maximum likelihood computer program. Ohio State Univ., Colarubus (Mimeo).
W.R. Harvey (1982) Mixed model capabilities of LSML76. Journal of Animal Science 54:1279-1285.
S.R. Searle, F.M. Speed, and G.A. Milliken (1980) Population marginal means in the linear model: An alternative to least squares means. The American Statistician 34 (4):216-221.
Searle, S. R., Milliken, G. A., and Speed, F. M. (1979). Expected Marginal Means in the Linear Model. Cornell University Biometrics Unit Technical Reports: Number BU-672-M.