We introduce the data-driven nested Operator Inference method for learning projection-based reduced-order models (ROMs) from snapshot data of high-dimensional dynamical systems. Projection-based ROMs exploit the intrinsic low-dimensionality of a full-order solution manifold. They typically 1) achieve significant computational savings, 2) guarantee approximation accuracy through established error theory, and 3) remain interpretable through the governing equations. However, constructing ROMs via projection requires access to the full-order operators -- a significant shortcoming for applications with legacy codes or commercial solvers. Operator Inference (OpInf) circumvents this requirement by learning the intrusive ROM from available full-order data and the structure of the governing equations. Under certain conditions, OpInf guarantees the exact reconstruction of the intrusive ROM, though meeting its data requirements in practice can be challenging, especially for highly non-linear operators because the degrees of freedom in the classic OpInf regression problem scale polynomially in the dimension of the reduced space. Consequently, classic OpInf requires precise regularization, balancing the numerical stability of its learning problem and the structural stability of its inferred ROM. In contrast, our nested OpInf approach partitions the learning problem into multiple regression problems that are provably better conditioned, thus alleviating the need for additional regularization. The partition is based upon a nested structure in the projection-based reduced-order matrices. It exploits a hierarchy in the reduced space's basis vectors to guarantee that the ROM's dominant dynamics are learned accurately. Since only few unknowns are learned at a time, nested OpInf is particularly applicable to higher-order polynomial systems. We demonstrate our method for the shallow ice equations with eighth-order polynomial operators.