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](https://arxiv.org/abs/1110.6430) 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](https://escholarship.org/content/qt98k5g03d/qt98k5g03d.pdf) 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](https://www.maths.ed.ac.uk/~v1ranick/papers/conwaysens.pdf)), 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}$.