The James Webb Space Telescope has ushered in an era of abundant high-redshift observations of young stellar populations characterized by strong emission lines, motivating us to integrate nebular emission into the new Maraston stellar population model which incorporates the latest Geneva stellar evolutionary tracks for massive stars with rotation. We use the photoionization code Cloudy to obtain the emergent nebular continuum and line emission for a range of modelling parameters, then compare our results to observations on various emission line diagnostic diagrams. We carry out a detailed comparison with several other models in the literature assuming different input physics, including modified prescriptions for stellar evolution and the inclusion of binary stars, and find close agreement in the H\(\rm \beta\), H\(\rm \alpha\), [N II]\(\lambda 6583\), and [S II]\(\lambda 6731\) luminosities between the models. However, we find significant differences in lines with high ionization energies, such as He II\(\lambda\)1640 and [O III]\(\lambda 5007\), due to large variations in the hard ionizing photon production rates. The models differ by a maximum of \(\hat{Q}_{\rm [O III]\lambda 5007} = \rm 6 \times 10^9 \; s^{-1} \, M_{\odot}^{-1}\), where these differences are mostly caused by the assumed stellar rotation and effective temperatures for the Wolf Rayet phase. Interestingly, rotation and uncorrected effective temperatures in our single star population models alone generate [O III] ionizing photon production rates higher than models including binary stars with ages between 1 to 8 Myr. These differences highlight the dependence of derived properties from SED fitting on the assumed model, as well as the sensitivity of predictions from cosmological simulations.