We present two observations of the high-mass X-ray binary GX 301−2 with NuSTAR, taken at different orbital phases and different luminosities. We find that the continuum is well described by typical phenomenological models, like a very strongly absorbed NPEX model. However, for a statistically acceptable description of the hard X-ray spectrum we require two cyclotron resonant scattering features (CRSF), one at ∼35 keV and the other at ∼50 keV. Even though both features strongly overlap, the good resolution and sensitivity of NuSTAR allows us to disentangle them at ≥99.9% significance. This is the first time that two CRSFs have been seen in GX 301−2. We find that the CRSFs are very likely independently formed, as their energies are not harmonically related and, if the observed feature were due to a single line, the deviation from a Gaussian shape would be very large. We compare our results to archival Suzaku data and find that our model also provides a good fit to those data. We study the behavior of the continuum as well as the CRSF parameters as function of pulse phase in seven phase bins. We find that the energy of the 35 keV CRSF varies smoothly as a function of phase, between 30 and 38 keV. To explain this variation, we apply a simple model of the accretion column, taking into account the altitude of the line-forming region, the velocity of the in-falling material, and the resulting relativistic effects. We find that in this model the observed energy variation can be explained as being simply due to a variation of the projected velocity and beaming factor of the line-forming region towards us.