Despite the long-lasting interest in the synthesis control of nanoparticles (NPs) in both fundamental and applied nanosciences, the driving mechanisms responsible for their size and shape selectivity in an environment (solution) are not completely understood, and a clear assessment of the respective roles of equilibrium thermodynamics and growth kinetics is still missing. In this study, relying on an efficient atomistic computational approach, we decipher the dependence of energetics, shapes and morphologies of gold NPs on the strength and nature of the metal-environment interaction. We highlight the conditions under which the energy difference between isotropic and elongated gold NPs is reduced, thus prompting their thermodynamic coexistence. The study encompasses both monocrystalline and multi-twinned particles and extends over size ranges particularly representative of the nucleation and early growth stages. Computational results are further rationalized with arguments involving the dependence of facet and edge energies on the metal-environment interactions. We argue that by determining the abundance and diversity of particles nucleated in solution, thermodynamics may constitute an important bias influencing their final shape. The present results provide firm grounds for kinetic simulations of particle growth.