After revisiting my question, I think I have managed to find a proof that there are no other examples of lattices with the requested property. I'm posting it as a self-answer in case someone is interested. Here is the outline:
The case of even dimension $\mathrm{dim}\: \Lambda \ge 4$ was already described$^\dagger$ in the "Attempt" section of my question; the most important part (bounding the generalized Bernoulli numbers) is done in Section 3 of this 2011 paper by M. Johnson, which I cited above.
We can also adapt the strategy in that paper to the odd-dimensional case. As in the even case, we start by observing that the theta series must be a certain Eisenstein series of half-integer weight and prescribed level $4N$ and character $\chi$, normalized such that the zeroth coefficient is $1$. One needs an expression for the Fourier coefficients of these series, which I found for $\mathrm{dim}\: \Lambda = 2k \ge 5$ ($k\in \frac12 \mathbb{Z}$) in page 17 of this 2015 thesis by M. Owen. We are only interested in the coefficient of $q^1$, which simplifies to
$$c_1 = \left(\frac{2\pi}{Ni}\right)^k \Gamma(k)^{-1} A(1) X(1),$$
where
$$A(1) = \sum_{r=1}^{4N} \epsilon_r^{2k} \left(\frac{4N}{r}\right) \chi(r) e^{2\pi ir/N}, \quad X(1)= \prod_{l\nmid 4N}(1+\chi(l) \epsilon_l^{2k-1} l^{1/2-k})$$
(here $l$ is prime and $\epsilon_n$ denotes the principal branch of $\sqrt{\left(\frac{-1}{n}\right)}$. Note also that I use $4N$ where the author uses $N$). This $c_1$ is an analogue of the generalized Bernoulli numbers for half-integral weight.
The $A(1)$ factor is a sum of roots of unity whose modulus is obviously bounded by $4N$, and by the same argument as in page 9 of Johnson's paper, $|X(1)|$ is bounded by $\prod_{l}(1+l^{1/2-k})=\frac{\zeta(k-1/2)}{\zeta(2k-1)}$. Putting everything together, we see that the bound on $|c_1|$ will be lower than $2$ for high enough $N$ and $k$ (concretely we have either $N=1$ and $k \le 9/2$ or $1 < N \le 3$ and $k=5/2$), so recalling that $c_1$ must be a nonzero even integer, we can restrict our analysis to these cases only. All of them are then ruled out by manual computation of the coefficients.
The case $\mathrm{dim}\: \Lambda = 2$ is easily dealt with by using the relationship between 2D integer lattices and rings of integers of imaginary quadratic number fields, and applying the class number formula.
This only leaves the case $\mathrm{dim}\: \Lambda = 3$. It is known that there are no isospectral lattices in dimension lower than $4$ (see e.g. the end of Section 2 here), so in this case it suffices to check the single-class examples in dimension $3$ with the LMFDB online tool, which I already did (this could also be an alternative way to treat the case $\mathrm{dim}\: \Lambda = 2$). This completes the proof.
$\dagger$- The extra possible example of dimension $4$ and level $42$ that I mentioned does not correspond to any lattice; perhaps the easiest way to see this is to inspect the coefficients of the candidate theta series, which starts as
$$1+2+2q^2+2q^3+2q^4+12q^5+2q^6+2q^7+2q^8+2q^9+12q^{10}+\ldots$$
(the general coefficient is twice the sum of the divisors of $n$ that don't divide $42$).
Such a lattice would have a vector $\mathbf{a}$ of squared norm $2$ and another linearly independent vector $\mathbf{b}$ of squared norm $10$; we can also see that all vectors of squared norm lower than $10$, and those strictly between $10$ and $20$, must be integer multiples of $\mathbf{a}$, since there are only two vectors of each norm.
The squared norm of the sum $\mathbf{a}+\mathbf{b}$ is a positive even integer bounded by $(\sqrt{2}+\sqrt{10})^2 \approx 20.944$; by the above observations, it can only be $20$. But then by the parallelogram law, $|\mathbf{a}-\mathbf{b}|^2 = 2|\mathbf{a}|^2+2|\mathbf{b}|^2-|\mathbf{a}+\mathbf{b}|^2 = 4+20-20=4$, which is a contradiction since $\mathbf{a}-\mathbf{b} \not\propto \mathbf{a}$.