Have the explicit Poisson-type formulas of Guinand and Meyer been observed before? In a recent paper of Meyer Measures with locally finite support and spectrum PNAS vol. 113 no. 12:3152–3158 (behind a paywall, but see also these seminar notes) some new explicit Poisson-type formulas are presented (the existence of such new formulas had been established by Lev and Olevskii). 
The text introducing the paper of Meyer says :

An important problem in harmonic analysis is solved in this article: Is the Poisson summation formula unique or does it belong to a wider class? The latter is true. The method that is used to prove this statement is surprising. Our new Poisson’s formulas were hidden inside an old and almost forgotten paper published in 1959 by A. P. Guinand. The role of number theory in this issue is fascinating. 

Indeed, the function $r_3$ that counts the number of decompositions of an integer into a sum of three squares makes an appearance. 
While that paper of Guinand was only cited 3 times as of 1999 according to Google Scholar, one of the three is in a paper of Berndt Character analogues of the Poisson and Euler-MacLaurin summation formulas with applications JNT Volume 7, Issue 4, November 1975, Pages 413–445 which was much more cited. Yet it appears that the type of formulas obtained by Berndt are different from those of Meyer. 
Hence the question in the title : 

1) do you know of any reference (published before 2016) where a Poisson-type formula involving $r_3$ appears in number theory ?

To which I'll add :

2) are these type of formulas surprising to a number theorist, or looking rather natural ?

 A: Here's an recent result that comes close [1]  The Fourier transform takes dirac combs supported on a lattices $\Gamma$ to a dirac comb supported on the dual lattices:
$$ \widehat{\delta_\Gamma} = |\Gamma| \cdot \delta_{\Gamma^*} $$
The dual to $\mathbb{Z}^2$ is itself.  Dual to the triangular lattice $\mathbb{Z}[\omega]$ with $\omega^3 = 1$ is another hexagonal lattice.
Then look at numbers which are sums of two squares as distances to the origin:
$$D_\square = \{ \sqrt{m^2 + n^2} : m, n \in \mathbb{Z} \} = \{ 0,1,\sqrt{2}, 2, \sqrt{5}, 2 \sqrt{2}, 3, \dots \} $$
How many lattice points on there of a circle of radius $r$?
Then  $\eta_\square(r) = 4 \times \sigma_{ \mathbb{Z}[i]}(r^2)$ the number of divisors of $r^2$ in $\mathbb{Z}[i]$.  In this measure we are going to weight each circle of radius $r$ by the number of lattice points.
$$ \sum_{r \in D_\square}\eta_\square(r) \mu_r = \left[ \sum_{r \in D_\square}\eta_\square(r) \mu_r \right]^\hat{} = \sum_{r \in D_\square}\eta_\square(r) J_0(2\pi |k| r) $$
This Fourier transform is in two dimensions (in $\mathbb{R}^2$).  I would have to read more carefully that this function is its own Fourier transform.
The rest of that paper deals with sets of points - arising from tilings - which are circularly symmetric on average and possible number theoretic consequences there.

Yves Meyer does not rotate the tilings, but deals with the lattice points themeselves.  He says cut-and-project tilings will not lead to crystalline measures and Poisson summation formulas where the spectrum is locally finite.
