I think there is indeed a "right" proof of Riemann's quartic theta relation. Namely, as a corollary of Mumford's first and second "fundamental identities" (described in 6.4, 6.5 of Tata III). Mumford shows us that there is a simple machinery (thetas with quadratic forms) producing an infinite number of theta relations. This abundance of theta relations is not a totally obvious fact (at least not until one finds the fundamental identities).

Riemann's quartic theta relation is apparently one very remarkable identity, or at least as remarkable as the $4\times 4$ rational symmetric positive definite and orthogonal matrix $$A=\frac{1}{2}\begin{pmatrix} 1 & 1 & 1 & 1 \cr 1 &1 &-1&-1 \cr 1&-1&1&-1 \cr 1&-1&-1&1 \end{pmatrix}.$$

Let us displace the question of "proving" riemann's quartic relation with confirming the identity, i.e. as a physical process of displacing and shuffling thetas. Then we should look maybe for a reason why $A$ is so remarkable. Or, at least see where else $A$ occurs.

But let's agree on *what* a theta function is. This requires replacing the upper half plane $\mathbb{H}$ consisting of $z=x+iy$ with $y>0$ with Siegel's upper half space $\mathbb{H}_g$ consisting of symmetric $g\times g$ complex matrices $T$ whose imaginary part is positive definite. An element $T\in \mathbb{H}_g$ yields a hermitian form on $\mathbb{C}^g$, by defining $T[z]:={}^tzTz$. Likewise for any $g\times h$ matrix $N$ we set $T[N]={}^t N T N$. It's also useful for a ring $R$ to let $R(g,h)$ denote all $g\times h$ matrices over $R$.

Computing from Mumford's second identity requires (explicit) coset representatives for the quotient $$\mathbb{Z}(g,4){}^t A / \mathbb{Z}(g,4){}^tA \cap \mathbb{Z}(g,4).$$ But the rational matrix $A$ is peculiar enough to have $A={}^tA=A^{-1}$ leading the coset representatives to take the form $(\eta, \eta, \eta, \eta)$, where $\eta$ varies over $\frac{1}{2}\mathbb{Z}^g / \mathbb{Z}^g$.

Identifying the coset representatives with the diagonal in $$\frac{1}{2}\mathbb{Z}^g / \mathbb{Z}^g \times \frac{1}{2}\mathbb{Z}^g / \mathbb{Z}^g \times\frac{1}{2}\mathbb{Z}^g / \mathbb{Z}^g \times\frac{1}{2}\mathbb{Z}^g / \mathbb{Z}^g $$ is perhaps reason enough to affirm $A$'s singularity. The identification of this diagonal with these cosets may in fact have some other root system interpretation (in sense of Jeff Harvey's answer involving string things).

But in a sense, after Mumford, it is the matrix $$A=\begin{pmatrix} 1 &1 \cr 1&-1 \cr \end{pmatrix}$$ which is responsible for riemann's quartic relation. This is Theorem 7.4 in Tata III.

Now a comment regarding the OP's secondary question concerning, for lack of a better expression, *a priori* knowledge of holomorphic functions on the complex torus $\mathbb{C}^g / \mathbb{Z}^g + T \mathbb{Z}^g$, for $T \in \mathbb{H}_g$.

At the beginning of the theory, by whatever divine interference, we define a theta function $\theta:\mathbb{C}^g \times \mathbb{H}_g \to \mathbb{C}$ by the expression
$$\theta(z,T)=\sum_{n \in \mathbb{Z}^g} \exp \pi i(T[n]+2{}^tnz).$$ In maintaining one's honesty (and i can only speak for myself) it is worth thinking about the necessity of our condition $im(T)>0$ in guaranteeing that $\theta(z,T)$ is a holomorphic function in $z$. One should also derive the periodicity and quasi-periodicity relations of $\theta(z,T)$ with respect to translations by $\mathbb{Z}^g$ and $T\mathbb{Z}^g$ in $\mathbb{C}^g$.

From the definition we find ourselves with a holomorphic function on $\mathbb{C}^g$ periodic under translations by $\mathbb{Z}^g$ and satisfying the quasi-periodicity relation $$f(z+Tm) =\exp \pi i(-T[m] - 2{}^t mz) f(z).$$ The following observation is already contained on pp.121 of Tata I.

If $f:\mathbb{C}^g \to \mathbb{C}$ is any entire function satisfying the quasiperiodicity relations of $\theta(z,T)$, then $f= c\theta(z,T)$ for some constant $c\in \mathbb{C}$. Indeed $\mathbb{Z}^g$-periodicity permits a fourier expansion $$f(z)=\sum_{n\in \mathbb{Z}^g} c_n \exp 2\pi i {}^tn z.$$ Now quasiperiodicity permits us to derive recursive relations amongst the coefficients $c_n, \\ n\in \mathbb{Z}^g$. Explicitly, we can relate $c_{n}$ to $c_{n + \epsilon_k}$, where $\epsilon_k$ is the $k^{th}$ unit vector. It's amusing to think of how running around the lattice generated by $\mathbb{Z}^g T$ varies the cofficients in the fourier expanision.

So knowing that $f$ is a quasiperiodic holomorphic function already yields the formula for the theta function $\theta(z,T)$ -- whose usual formula we should recognize already as a fourier expansion. But I also know nothing about $\theta(z,T)$ as eigenfunction.

Now concerning "promises" relating the Heisenberg representations to facts on theta functions. From what I can tell, the only promise is Stone-von Neumann-Mackey's theorem on irreducible unitary representations (which are faithful on the centre) of the Heisenberg group: all such representations are equivalent. So given two representations (and there are many), there is an "intertwining" among them. But for some facts, one suffers. In 7.4 Tata III there is a description of riemann's quartic relation as arising from such an intertwining. However the Heisenberg group is defined over the finite adeles, i.e. $Heis(2g, \mathbb{A}_f)$. I would like to know where triality, spinors, or $so(8)$ interact in $\mathbb{A}_f$.