We report quantum Monte Carlo calculations of single-$Lambda$ hypernuclei for $A<50$ based on phenomenological two- and three-body hyperon-nucleon forces. We present results for the $Lambda$ separation energy in different hyperon orbits, showing that the accuracy of theoretical predictions exceeds that of currently available experimental data, especially for medium-mass hypernuclei. We show the results of a sensitivity study that indicates the possibility to investigate the nucleon-isospin dependence of the three-body hyperon-nucleon-nucleon force in the medium-mass region of the hypernuclear chart, where new spectroscopy studies are currently planned. The importance of such a dependence for the description of the physics of hypernuclei, and the consequences for the prediction of neutron star properties are discussed.