Discrete eigenvalues can always be labeled by the integers, so we can always pick up $E_n$ with $n$ starting from $0$ or $1$, it is just a matter of convenction. I suppose this is where your question originates. For the harmonic oscillator the energies are given by $\hbar \omega (m +1/2)$ with $m=0,1,\ldots$ so it is natural to pick $n=m$ starting from $0$.
For the particle in the box, note that the problem is equivalent to a "free" particle with certain boundary conditions (wavefunction zero at the border). In this sense it is no surprise that "momentum" is conserved (people speak of quasi-momentum in this case). In this picture the Hamiltonian is simply
$$H = \frac{p^2}{2m}.$$
If you set the box as in $[0,L]$ you can write the eigenfunctions as
$$\psi_k(x) \propto \sin(kx)$$
where $k$ is this quasi-momentum to be determined. In this way one of the boundary conditions is automatically satisfied. To satisfy the other one we must impose $\sin(Lk)=0$ which in principle leads to $k= m\pi/L$, with $m\in \mathbb{Z}$. However, $k$ and $-k$ single out the same quantum state (a quantum state is an equivalence class of vectors which differ by a phase). Moreover $k=0$ is not allowed because the zero wave-function is not normalizable.
Hence we are led to
$$k= \frac{n\pi}{L}, \ \ n = 1,2,\ldots$$
In a sense, physically, this tells us that the state with zero momentum (a state that does not move left or right) is not allowed. This is clearly a quantum mechanical effect (a fixed particle is perfectly allowed classically). You can think of its origin in the Heisenberg's uncertainty principle. Since $H\propto p^2$ we can obviously label the energies with the same label that we use for the momentum (another good notation would be simply to label them with the quasi-momentum itself as in $E_k$ and remember what is the spectrum of $p$, $\sigma(p)$.
I suppose you knew all of this but I don't think you can say much more.
Edit
Let me briefly comment on the number of nodes of the wavefunctions.
For a Hamiltonian of the form
$$
H= \frac{p^2}{2m} + V(x) = - \frac{\hbar^2}{2m} \frac{\partial^2}{\partial^2 x} +V(x)
$$
a theorem by Courant and Hilbert (generally called nodal theorem) states that the ground state wave function has no nodes (can be chosen to be everywhere non-negative) and the $n$-th level has precisely $n-1$ nodes. Moreover if the potential goes to infinity as $|x| \to \infty$, then the eigenvalues form a discrete unbounded sequence. I imagine variations of this theorem exist for $\mathbb{R}^d$. So the behavior you observed for the number of nodes is in fact quite general.