We report extensive numerical simulations of growth models belonging to the nonlinear molecular beam epitaxy (nMBE) class, on flat (fixed-size) and expanding substrates (ES). In both $d=1+1$ and $2+1$, we find that growth regime height distributions (HDs), and spatial and temporal covariances are universal, but are dependent on the initial conditions, while the critical exponents are the same for flat and ES systems. Thus, the nMBE class does split into subclasses, as does the Kardar-Parisi-Zhang (KPZ) class. Applying the KPZ ansatz to nMBE models, we estimate the cumulants of the $1+1$ HDs. Spatial covariance for the flat subclass is hallmarked by a minimum, which is not present in the ES one. Temporal correlations are shown to decay following well-known conjectures.