We study the sign problem in the Hubbard model on the hexagonal lattice away from half-filling using the Lefschetz thimbles method. We identify the saddle points, reduce their amount, and perform quantum Monte Carlo (QMC) simulations using the holomorphic gradient flow to deform the integration contour in complex space. Finally, the results are compared against exact diagonalization (ED). We show that the sign problem can be substantially weakened, even in the regime with low temperature and large chemical potential, where standard QMC techniques exhibit an exponential decay of the average sign.