We address adaptive multivariate polynomial approximation by means of the discrete least-squares method with random evaluations, to approximate in the L2 probability sense a smooth function depending on a random variable distributed according to a given probability density. The polynomial least-squares approximation is computed using random noiseless pointwise evaluations of the target function. Here noiseless means that the pointwise evaluation of the function is not polluted by the presence of noise. Recent works Migliorati et al. (Found Comput Math 14:419–456, 2014), Cohen et al. (Found Comput Math 13:819–834, 2013), and Chkifa et al. (Discrete least squares polynomial approximation with random evaluations – application to parametric and stochastic elliptic PDEs, EPFL MATHICSE report 35/2013, submitted) have analyzed the univariate and multivariate cases, providing error estimates for (a priori) given sequences of polynomial spaces. In the present work, we apply the results developed in the aforementioned analyses to devise adaptive least-squares polynomial approximations. We build a sequence of quasi-optimal best n-term sets to approximate multivariate functions that feature strong anisotropy in moderately high dimensions. The adaptive approximation relies on a greedy selection of basis functions, which preserves the downward closedness property of the polynomial approximation space. Numerical results show that the adaptive approximation is able to catch effectively the anisotropy in the function.