This paper is concerned with two types of two-grid stabilized finite element methods (FEMs) based on Newton iteration for the steady-state nature convection problem. The first method needs to solve one small nonlinear natural convection system on the coarse mesh with mesh width $H$, and then to solve one large linearized natural convection system on the fine mesh with mesh width $h=\mathcal {O}(H^2)$ based on Newton iteration. The other method needs to solve one small nonlinear natural convection system on the same coarse mesh, and then to solve two large linearized systems on the fine mesh with mesh width $h=\mathcal {O}(H^{\frac{7-\varepsilon}{2}})$ based on Newton iteration which have the same stiffness matrix with only different right-hand side. In both methods, the stabilization terms are defined via two local Gauss integrations at element level which has no need to introduce additional variables comparing with the standard variational multiscale stabilized FEMs. The stability estimates and the convergence analysis for both methods are derived strictly. Ample numerical results are presented to confirm the theoretical predictions and demonstrate the efficiency of the new methods.