We construct a nonrelativistic effective field theory description of heavy quarkonium hybrids from QCD. We identify the symmetries of the system made of a heavy quark, a heavy antiquark, and glue in the static limit. Corrections to this limit can be obtained order by order in an expansion in the inverse of the mass $m$ of the heavy quark. At order $1/m$ in the expansion, we obtain at the level of potential Non-Relativistic QCD a system of coupled Schr\"odinger equations that describes hybrid spin-symmetry multiplets, including the mixing of different static energies into the hybrid states, an effect known as $\Lambda$-doubling in molecular physics. In the short distance, the static potentials depend on two nonperturbative parameters, the gluelump mass and the quadratic slope, which can be determined from lattice calculations. We adopt a renormalon subtraction scheme for the calculation of the perturbative part of the potential. We numerically solve the coupled Schr\"odinger equations and obtain the masses for the lowest lying spin-symmetry multiplets for $c\bar{c}$, $b\bar{c}$, and $b\bar{b}$ hybrids. The $\Lambda$-doubling effect breaks the degeneracy between opposite parity spin-symmetry multiplets and lowers the mass of the multiplets that get mixed contributions of different static energies. We compare our findings to the experimental data, direct lattice computations, sum rules calculations, and discuss the relation to the Born-Oppenheimer approximation.