[proofplan]
The Besov space $B^s_{p,q}(\mathbb{R}^n)$ has norm $\|f\|_{B^s_{p,q}} = \|(2^{js}\|\Delta_j f\|_{L^p})_{j \ge 0}\|_{\ell^q}$ relative to a fixed Littlewood--Paley resolution. Up to norm equivalence, this realises $B^s_{p,q}$ as the weighted vector-valued sequence space $\ell^q_s(L^p) := \ell^q(\mathbb{N}_0; 2^{js}; L^p)$, i.e. $L^p$-valued sequences weighted by $2^{js}$ in $\ell^q$. The Lions--Peetre $K$-method computes the real-interpolation space of a couple of Banach spaces in terms of its $K$-functional, and the Holmstedt formula computes the $K$-functional for weighted $\ell^q$ pairs explicitly. Combining these, the real-interpolation space $(B^{s_0}_{p,q_0}, B^{s_1}_{p,q_1})_{\theta, q}$ coincides up to equivalent norms with the weighted $\ell^q$ pair $(\ell^{q_0}_{s_0}(L^p), \ell^{q_1}_{s_1}(L^p))_{\theta, q}$, which by the Holmstedt formula equals $\ell^q_{s_\theta}(L^p) \cong B^{s_\theta}_{p,q}$. This four-step chain yields the asserted identification.
[/proofplan]
[step:Fix a Littlewood--Paley resolution and recall the Besov norm]
Fix $\phi_0, \phi \in C_c^\infty(\mathbb{R}^n)$ with $\phi_0$ supported in $\{|\xi| \le 2\}$ and equal to $1$ on $\{|\xi| \le 1\}$, $\phi$ supported in the annulus $\{1/2 \le |\xi| \le 2\}$, and
\begin{align*}
\phi_0(\xi) + \sum_{j=1}^\infty \phi(2^{-j}\xi) = 1 \qquad \text{for all } \xi \in \mathbb{R}^n.
\end{align*}
Set $\phi_j(\xi) := \phi(2^{-j}\xi)$ for $j \ge 1$. Define the Littlewood--Paley projectors
\begin{align*}
\Delta_j : \mathcal{S}'(\mathbb{R}^n) &\to \mathcal{S}'(\mathbb{R}^n) \\
f &\mapsto \mathcal{F}^{-1}\bigl(\phi_j\, \hat f\,\bigr), \qquad j \ge 0.
\end{align*}
For $s \in \mathbb{R}$, $1 \le p \le \infty$, $1 \le q \le \infty$, the Besov norm is
\begin{align*}
\|f\|_{B^s_{p,q}(\mathbb{R}^n)} := \Bigl\| \bigl(2^{js}\, \|\Delta_j f\|_{L^p(\mathbb{R}^n)}\bigr)_{j \ge 0} \Bigr\|_{\ell^q(\mathbb{N}_0)}.
\end{align*}
Different choices of resolution $(\phi_0, \phi)$ yield equivalent norms; the norm equivalence constants depend on the chosen resolution.
[/step]
[step:Construct an isometric (up to equivalence) embedding into a weighted vector-valued sequence space]
Define the weighted vector-valued $\ell^q$ space
\begin{align*}
\ell^q_s(L^p) := \Bigl\{ a = (a_j)_{j \ge 0} : a_j \in L^p(\mathbb{R}^n), \; \|a\|_{\ell^q_s(L^p)} := \bigl\|(2^{js}\|a_j\|_{L^p})_{j \ge 0}\bigr\|_{\ell^q} < \infty \Bigr\},
\end{align*}
which is a Banach space (with the obvious modification for $q = \infty$). Define the projection-and-extension pair
\begin{align*}
P : \mathcal{S}'(\mathbb{R}^n) &\to \prod_{j \ge 0} \mathcal{S}'(\mathbb{R}^n), & f &\mapsto (\Delta_j f)_{j \ge 0}, \\
S : \prod_{j \ge 0} \mathcal{S}'(\mathbb{R}^n) &\to \mathcal{S}'(\mathbb{R}^n), & (a_j)_{j \ge 0} &\mapsto \sum_{j=0}^\infty \widetilde\Delta_j a_j,
\end{align*}
where $\widetilde\Delta_j$ is convolution with $\mathcal{F}^{-1}\widetilde\phi_j$ for a fixed fattened bump $\widetilde\phi_j$ satisfying $\widetilde\phi_j \equiv 1$ on $\operatorname{supp}\phi_j$ and supported in a slightly enlarged annulus (so that, by the construction, $\widetilde\Delta_j \Delta_j = \Delta_j$ and $\Delta_j \widetilde\Delta_k = 0$ when $|j - k|$ is large enough relative to the resolution).
By the construction of the resolution, $\sum_j \Delta_j f = f$ in $\mathcal{S}'(\mathbb{R}^n)$, and the operator $S$ is a left inverse of $P$ on the closed subspace $\{a : a_j = \Delta_j a_j\} \subseteq \ell^q_s(L^p)$. The operator norm bounds
\begin{align*}
\|P f\|_{\ell^q_s(L^p)} = \|f\|_{B^s_{p,q}}, \qquad \|S a\|_{B^s_{p,q}} \le M\, \|a\|_{\ell^q_s(L^p)}
\end{align*}
hold for every $f \in B^s_{p,q}$ and every $a \in \ell^q_s(L^p)$, where $M = M(p, q, s, n)$ is a constant arising from the finite overlap of the fattened bumps $\widetilde\phi_j$. The first equality is the definition of the Besov norm; the second follows from a standard finite-overlap argument together with the Mihlin multiplier theorem applied uniformly in $j$ to the symbol $\widetilde\phi_j$. Hence $B^s_{p,q}(\mathbb{R}^n)$ is, up to equivalent norms, a complemented subspace of $\ell^q_s(L^p)$ via the projection-extension pair $(P, S)$.
[/step]
[step:Apply the retract principle to reduce to interpolation of weighted $\ell^q$ pairs]
Consider the Banach pair $(\ell^{q_0}_{s_0}(L^p), \ell^{q_1}_{s_1}(L^p))$. By the retract principle for real interpolation (see the [Retract Theorem for Real Interpolation](/theorems/???)): if $(P, S)$ is a retract pair such that $P$ is bounded $X_i \to Y_i$ and $S$ is bounded $Y_i \to X_i$ with $S P = \mathrm{id}_{X_i}$ for $i = 0, 1$, then $(X_0, X_1)_{\theta, q}$ is a retract of $(Y_0, Y_1)_{\theta, q}$ via the same pair, and the norm equivalence is preserved.
Applying this with $X_i = B^{s_i}_{p, q_i}$ and $Y_i = \ell^{q_i}_{s_i}(L^p)$ and the pair $(P, S)$ from the previous step (which is bounded with the appropriate norms for both $i = 0$ and $i = 1$ by the same construction), we conclude
\begin{align*}
\bigl(B^{s_0}_{p,q_0}, B^{s_1}_{p,q_1}\bigr)_{\theta, q} = \bigl(\ell^{q_0}_{s_0}(L^p), \ell^{q_1}_{s_1}(L^p)\bigr)_{\theta, q} \cap \mathrm{Image}(P) \quad \text{up to equivalent norms},
\end{align*}
and the equivalence constants depend only on $\theta$, $q$, the operator norms of $P, S$, and the parameters of the resolution.
It thus suffices to compute the right-hand interpolation space and its restriction to $\mathrm{Image}(P) = \{a \in \ell^q_s(L^p) : a_j = \Delta_j a_j\}$, then transfer back via $S$.
[/step]
[step:Compute the $K$-functional of the weighted $\ell^q(L^p)$ pair via the Holmstedt formula]
Let $A = \ell^{q_0}_{s_0}(L^p)$ and $B = \ell^{q_1}_{s_1}(L^p)$. For a couple of Banach spaces $(A, B)$ embedded in a common Hausdorff topological vector space, the [Lions--Peetre $K$-functional](/theorems/???) is
\begin{align*}
K : (A + B) \times (0, \infty) &\to [0, \infty) \\
(a, t) &\mapsto \inf \{\|a_0\|_A + t\, \|a_1\|_B : a = a_0 + a_1,\, a_0 \in A,\, a_1 \in B\}.
\end{align*}
For the weighted $\ell^q$ pair, the [Holmstedt Formula](/theorems/???) gives the explicit computation: for $a = (a_j)_{j \ge 0}$ in $A + B$,
\begin{align*}
K(a, t)^{\min(q_0, q_1)} \asymp \sum_{j : 2^{j(s_0 - s_1)} \ge t} \bigl(2^{js_0}\|a_j\|_{L^p}\bigr)^{q_0}\, \cdot \mathbb{1}_{q_0 < \infty} + t^{q_1} \!\! \sum_{j : 2^{j(s_0 - s_1)} < t} \bigl(2^{js_1}\|a_j\|_{L^p}\bigr)^{q_1} \cdot \mathbb{1}_{q_1 < \infty},
\end{align*}
with the obvious supremum modifications when $q_0 = \infty$ or $q_1 = \infty$. The implicit constants depend on $s_0, s_1, p, q_0, q_1$. The threshold scale is $j_0(t) := \log_2(t)/(s_0 - s_1)$, splitting the dyadic axis into "low scales" (where $\ell^{q_0}_{s_0}$ dominates) and "high scales" (where $\ell^{q_1}_{s_1}$ dominates).
The $K$-functional simplifies (still by Holmstedt) to
\begin{align*}
K(a, t) \asymp \biggl( \sum_{j \ge j_0(t)} \bigl(2^{js_0}\|a_j\|_{L^p}\bigr)^{q_0} \biggr)^{1/q_0} + t \cdot \biggl( \sum_{j < j_0(t)} \bigl(2^{js_1}\|a_j\|_{L^p}\bigr)^{q_1} \biggr)^{1/q_1},
\end{align*}
under the mild proviso $s_0 > s_1$ (the case $s_0 < s_1$ is symmetric, swapping the roles of $A$ and $B$).
[/step]
[step:Compute the real-interpolation norm via the $K$-method and identify it with $\ell^q_{s_\theta}(L^p)$]
By the [Lions--Peetre $K$-Method](/theorems/???), for $0 < \theta < 1$ and $1 \le q \le \infty$, the real-interpolation space $(A, B)_{\theta, q}$ has norm
\begin{align*}
\|a\|_{(A, B)_{\theta, q}} := \biggl( \int_0^\infty \bigl(t^{-\theta}\, K(a, t)\bigr)^q \frac{d\mathcal{L}^1(t)}{t} \biggr)^{1/q}, \qquad 1 \le q < \infty,
\end{align*}
with the obvious essential-supremum modification for $q = \infty$. Substituting the Holmstedt expression for $K(a, t)$ and changing variables $t = 2^{j(s_0 - s_1)}$ (giving $dt/t = (s_0 - s_1)\ln 2 \cdot dj$ as a discrete measure on $j \in \mathbb{Z}$), the integral $\int_0^\infty (t^{-\theta} K(a, t))^q\, dt/t$ becomes a discrete sum over $j \in \mathbb{Z}$ of terms of the form
\begin{align*}
\bigl(2^{j(s_0 - s_1)}\bigr)^{-\theta q}\, \cdot \bigl[\text{Holmstedt expression at threshold } j\bigr]^q.
\end{align*}
After unwinding (this is the standard computation that the Holmstedt formula was designed for), the Lions--Peetre $K$-method on a weighted $\ell^q$ pair recovers the weighted $\ell^q$ space at the interpolated weight:
\begin{align*}
(A, B)_{\theta, q} = \ell^q_{s_\theta}(L^p) \quad \text{with equivalent norms}, \qquad s_\theta = (1 - \theta)\, s_0 + \theta\, s_1.
\end{align*}
The equivalence constants depend on $\theta, q, q_0, q_1, p, s_0, s_1$; they are finite because $1 \le q_0, q_1, q \le \infty$ and $s_0 \ne s_1$ (which makes the change of variable $t \leftrightarrow j$ non-degenerate, with Jacobian $|s_0 - s_1| > 0$).
[/step]
[step:Transfer the identification back to Besov spaces via the retract pair $(P, S)$]
By the retract principle invoked in the third step,
\begin{align*}
(B^{s_0}_{p,q_0}, B^{s_1}_{p,q_1})_{\theta, q} = S\bigl((\ell^{q_0}_{s_0}(L^p), \ell^{q_1}_{s_1}(L^p))_{\theta, q}\bigr) \cap \mathcal{S}'(\mathbb{R}^n)
\end{align*}
with equivalent norms (the retract maps on both sides are the same $(P, S)$ from the second step). By the previous step, the right-hand side equals $S(\ell^q_{s_\theta}(L^p)) \cap \mathcal{S}'(\mathbb{R}^n)$, which by the construction of $S$ is exactly $B^{s_\theta}_{p, q}(\mathbb{R}^n)$ (with norm equivalence constant depending only on the resolution).
Concretely: a tempered distribution $f \in (B^{s_0}_{p, q_0}, B^{s_1}_{p, q_1})_{\theta, q}$ corresponds via $P$ to its Littlewood--Paley sequence $a = (\Delta_j f)_{j \ge 0}$, which lies in $(\ell^{q_0}_{s_0}(L^p), \ell^{q_1}_{s_1}(L^p))_{\theta, q} = \ell^q_{s_\theta}(L^p)$. The $\ell^q_{s_\theta}(L^p)$-norm of this sequence is, by definition, $\|f\|_{B^{s_\theta}_{p, q}}$. Conversely, $f \in B^{s_\theta}_{p, q}$ has Littlewood--Paley sequence in $\ell^q_{s_\theta}(L^p) = (\ell^{q_0}_{s_0}(L^p), \ell^{q_1}_{s_1}(L^p))_{\theta, q}$, and applying $S$ recovers $f$ in $(B^{s_0}_{p, q_0}, B^{s_1}_{p, q_1})_{\theta, q}$.
Hence
\begin{align*}
(B^{s_0}_{p, q_0}(\mathbb{R}^n), B^{s_1}_{p, q_1}(\mathbb{R}^n))_{\theta, q} = B^{s_\theta}_{p, q}(\mathbb{R}^n)
\end{align*}
with norms equivalent up to constants depending only on $n, s_0, s_1, p, q_0, q_1, q, \theta$, and the resolution of unity. This is the asserted identification.
[/step]