We present numerical simulations of binary neutron star mergers, comparing irrotational binaries to binaries of NSs rotating aligned to the orbital angular momentum. For the first time, we study spinning BNSs employing nuclear physics equations of state, namely the ones of Lattimer and Swesty as well as Shen, Horowitz, and Teige. We study mainly equal mass systems leading to a hypermassive neutron star (HMNS), and analyze in detail its structure and dynamics. In order to exclude gauge artifacts, we introduce a novel coordinate system used for post-processing. The results for our equal mass models show that the strong radial oscillations of the HMNS modulate the instantaneous frequency of the gravitational wave (GW) signal to an extend that leads to separate peaks in the corresponding Fourier spectrum. In particular, the high frequency peaks which are often attributed to combination frequencies can also be caused by the modulation of the m=2 mode frequency in the merger phase. As a consequence for GW data analysis, the offset of the high frequency peak does not necessarily carry information about the radial oscillation frequency. Further, the low frequency peak in our simulations is dominated by the contribution of the plunge and the first 1-2 bounces. The amplitude of the radial oscillations depends on the initial NS spin, which therefore has a complicated influence on the spectrum. Another important result is that HMNSs can consist of a slowly rotating core with an extended, massive envelope rotating close to Keplerian velocity, contrary to the common notion that a rapidly rotating core is necessary to prevent a prompt collapse. Finally, our estimates on the amount of unbound matter show a dependency on the initial NS spin, explained by the influence of the latter on the amplitude of radial oscillations, which in turn cause shock waves.