We study deviations between MSSM and $Z_3$-invariant NMSSM, with respect to their predictions in $\Delta F=2 $ processes. We find that potentially significant effects arise either from the well known double-penguin diagrams, due to the extra scalar NMSSM states, or from neutralino-gluino box contributions, due to the extended neutralino sector. Both are discussed to be effective in the large $\tan\beta$ regime. Enhanced genuine-NMSSM contributions in double penguins are expected for a light singlet spectrum (CP-even,CP-odd), while the magnitude of box effects is primarily controlled through singlino mixing. The latter is found to be typically subleading (but non-negligible) for $\lambda \lesssim 0.5$, however it can become dominant for $\lambda\sim \mathcal{O}(1)$. We also study the low $\tan\beta$ regime, where a distinction between MSSM and NMSSM can come instead due to experimental constraints, acting differently on the allowed parameter space of each model. To this end, we incorporate the LHC Run-I limits from $H\rightarrow Z Z$, $A \rightarrow hZ$ and $H^\pm \rightarrow \tau \nu $ non-observation along with Higgs observables and set (different) upper bounds for new physics contributions in $\Delta F=2 $ processes. We find that a $\sim 25\%$ contribution in $\Delta M_{s(d)}$ is still possible for MFV models, however such a large effect is nowadays severely constrained for the case of MSSM, due to stronger bounds on the charged Higgs masses.