We study a three-orbital Hubbard-Kanamori model relevant for iron-based superconductors using variational wave functions explicitly including spatial correlations and electron pairing. We span the nonmagnetic sector from filling n = 4, which is representative of undoped iron-based superconductors, to n = 3, where a Mott insulating state with each orbital at half filling is found. In the strong-coupling regime, when the electron density is increased, we find a spontaneous differentiation between the occupation of dxz and dyz orbitals, leading 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 dxy 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.
