A framework for low-order predictive statistical modeling and uncertainty quantification in turbulent dynamical systems is developed here. These reduced-order, modified quasilinear Gaussian (ROMQG) algorithms apply to turbulent dynamical systems in which there is significant linear instability or linear nonnormal dynamics in the unperturbed system and energy-conserving nonlinear interactions that transfer energy from the unstable modes to the stable modes where dissipation occurs, resulting in a statistical steady state; such turbulent dynamical systems are ubiquitous in geophysical and engineering turbulence. The ROMQG method involves constructing a low-order, nonlinear, dynamical system for the mean and covariance statistics in the reduced subspace that has the unperturbed statistics as a stable fixed point and optimally incorporates the indirect effect of non-Gaussian third-order statistics for the unperturbed system in a systematic calibration stage. This calibration procedure is achieved through information involving only the mean and covariance statistics for the unperturbed equilibrium. The performance of the ROMQG algorithm is assessed on two stringent test cases: the 40-mode Lorenz 96 model mimicking midlatitude atmospheric turbulence and two-layer baroclinic models for high-latitude ocean turbulence with over 125,000 degrees of freedom. In the Lorenz 96 model, the ROMQG algorithm with just a single mode captures the transient response to random or deterministic forcing. For the baroclinic ocean turbulence models, the inexpensive ROMQG algorithm with 252 modes, less than 0.2% of the total, captures the nonlinear response of the energy, the heat flux, and even the one-dimensional energy and heat flux spectra.dynamical systems with many instabilities | nonlinear response and sensitivity | reduced-order modified quasilinear Gaussian closure T urbulent dynamical systems are characterized by both a large dimensional phase space and a large dimension of instabilities (i.e., a large number of positive Lyapunov exponents on the attractor). They are ubiquitous in many complex systems with fluid flow such as, for example, the atmosphere, ocean, and coupled climate system, confined plasmas, and engineering turbulence at high Reynolds numbers. In these, linear instabilities are mitigated by energy-conserving nonlinear interactions that transfer energy to the linearly stable modes where it is dissipated, resulting in a statistical steady state. Uncertainty quantification (UQ) in turbulent dynamical systems is a grand challenge where the goal is to obtain statistical estimates such as the change in mean and variance for key physical quantities in the nonlinear response to changes in external forcing parameters or uncertain initial data. These key physical quantities are often characterized by the degrees of freedom that carry the largest energy or variance, and an even more ambitious grand challenge is to develop truncated low-order models for UQ for a reduced set of important variables with the largest variance. Thi...