We consider the least-squares variational kernel-based methods for numerical solution of partial differential equations. Indeed, we focus on least-squares principles to develop meshfree methods to find the numerical solution of a general second order ADN elliptic boundary value problem in domain Ω ⊂ R d under Dirichlet boundary conditions. Most notably, in these principles it is not assumed that differential operator is self-adjoint or positive definite as it would have to be in the Rayleigh-Ritz setting. However, the new scheme leads to a symmetric and positive definite algebraic system allowing us to circumvent the compatibility conditions arising in standard and mixed-Galerkin methods. In particular, the resulting method does not require certain subspaces satisfying any boundary condition.The trial space for discretization is provided via standard kernels that reproduce H τ (Ω), τ > d/2, as their native spaces. Therefore, the smoothness of the approximation functions can be arbitrary increased without any additional task. The solvability of the scheme is proved and the error estimates are derived for functions in appropriate Sobolev spaces. For the weighted discrete least-squares principles, we show that the optimal rate of convergence in L 2 (Ω) is accessible. Furthermore, for d ≤ 3, the proposed method has optimal rate of convergence in H k (Ω) whenever k ≤ τ . The condition number of the final linear system is approximated in terms of discterization quality. Finally, the results of some computational experiments support the theoretical error bounds.