  We study a three-orbital Hubbard-Kanamori model relevant for iron-based
superconductors using variational wave functions, which explicitly include
spatial correlations and electron pairing. We span the nonmagnetic sector
ranging from a filling $n=4$, which is representative of undoped iron-based
superconductors, to $n=3$. In the latter case, a Mott insulating state is
found, with each orbital at half filling. In the strong-coupling regime, when
the electron density is increased, we find a spontaneous differentiation
between the occupation of $d_{xz}$ and $d_{yz}$ orbitals, which leads to an
orbital-selective state with a nematic character that becomes stronger at
increasing density. One of these orbitals stays half filled for all densities
while the other one hosts (together with the $d_{xy}$ orbital) the excess of
electron density. Most importantly, in this regime, long-range pairing
correlations appear in the orbital with the largest occupation. Our results
highlight a strong link between orbital-selective correlations, nematicity, and
superconductivity, which requires the presence of a significant Hund's
coupling.
