Systems of coupled optical parametric oscillators (OPOs) forming an Ising machine are emerging as large-scale simulators of the Ising model. The advances in computer science and nonlinear optics have triggered not only the physical realization of hybrid (electro-optical) or all-optical Ising machines, but also the demonstration of quantum-inspired algorithms boosting their performances. To date, the use of the quantum nature of parametrically generated light as a further resource for computation represents a major open issue. A key quantum feature is the non-Gaussian character of the system state across the oscillation threshold. In this paper, we perform an extensive analysis of the emergence of non-Gaussianity in the single quantum OPO with an applied external field. We model the OPO by a Lindblad master equation, which is numerically solved by an ab initio method based on exact diagonalization. Non-Gaussianity is quantified by means of three different metrics: Hilbert-Schmidt distance, quantum relative entropy, and photon distribution. Our findings reveal a nontrivial interplay between parametric drive and applied field: (i) Increasing pump monotonously enhances non-Gaussianity, and (ii) Increasing field first sharpens non-Gaussianity, and then restores the Gaussian character of the state when above a threshold value.