As current gravitational wave (GW) detectors increase in sensitivity, and particularly as new instruments are being planned, there is the possibility that ground-based GW detectors will observe GWs from highly eccentric neutron star binaries. We present the first detailed study of highly eccentric BNS systems with full (3+1)D numerical relativity simulations using consistent initial conditions, i.e., setups which are in agreement with the Einstein equations and with the equations of general relativistic hydrodynamics in equilibrium. Overall, our simulations cover two different equations of state (EOSs), two different spin configurations, and three to four different initial eccentricities for each pairing of EOS and spin. We extract from the simulated waveforms the frequency of the f-mode oscillations induced during close encounters before the merger of the two stars. The extracted frequency is in good agreement with f-mode oscillations of individual stars for the irrotational cases, which allows an independent measure of the supranuclear equation of state not accessible for binaries on quasicircular orbits. The energy stored in these f-mode oscillations can be as large as 10−3  M⊙∼1051  erg, even with a soft EOS. In order to estimate the stored energy, we also examine the effects of mode mixing due to the stars’ offset from the origin on the f-mode contribution to the GW signal. While in general (eccentric) neutron star mergers produce bright electromagnetic counterparts, we find that for the considered cases with fixed initial separation the luminosity decreases when the eccentricity becomes too large, due to a decrease of the ejecta mass. Finally, the use of consistent initial configurations also allows us to produce high-quality waveforms for different eccentricities which can be used as a test bed for waveform model development of highly eccentric binary neutron star systems.