The origin of the nonlinear solvatochromic shift of para-nitroaniline was investigated using a mean-field sequential QM/MM method, with electron transitions computed at the CASPT2/cc-pVDZ level. Experimental data shows that the solvatochromic shift has a strong nonlinear behavior in certain solvent mixtures. We studied the case of cyclohexane-triethylamine mixtures. The results are in good agreement with the experiments and correctly reproduce the nonlinear variation of the solvent shift. Preferential solvation is clearly observed, where the local solvent composition in the neighborhood of the solute is significantly different from the bulk. It is found that even at low triethylamine concentrations a strong hydrogen bond is formed between para-nitroaniline and triethylamine, and cyclohexane is practically absent from the first solvation layer already at a molar fraction of 0.6 in triethylamine. The hydrogen bond formed is sufficiently long-lived to determine an asymmetric environment around the solute molecule. The resulting nonlinear solvent effect is mainly due to this hydrogen bond influence, although there is also a small contribution from dielectric enrichment.