Without the trick in the answer by fretty (essentially recognising the expression asked for can be rewritten as $f(\alpha+\beta+\gamma+\delta)$) you can obtain the result through the approach you had originally chosen. For that it suffices to express $(\alpha+\beta+\gamma)(\alpha+\beta+\delta)(\alpha+\gamma+\delta)(\beta+\gamma+\delta)$, which is clearly a symmetric polynomial of $\alpha,\beta,\gamma,\delta$, in terms of the elementary symmetric polynomials (the right hand sides of your four equations; call them $e_1,e_2,e_3,e_4$ respectively).
Expansion of the product gives $3^4=81$ terms, which can be grouped into minimal (also called monomial) symmetric polynomials, the sum of all distinct permutations of a given monomial. I will denote these by $m_\lambda$ where $\lambda$ is the pattern of exponents; for instance $m_{2,1,1}$ is the sum of monomials like $\alpha^2\beta\gamma$ or $\alpha\gamma^2\delta$. With this notation one gets
$$
(\alpha+\beta+\gamma)
(\alpha+\beta+\delta)
(\alpha+\gamma+\delta)
(\beta+\gamma+\delta)=
m_{3,1}+2m_{2,2}+4m_{2,1,1}+9m_{1,1,1,1}.
$$
Now each minimal symmetric polynomial contains the leading term of a unique product of elementary symmetric polynomials, as explained here, and this leads to relations
$$\begin{align}
m_{3,1}&=e_2e_1^2 - 2m_{2,2} - 5m_{2,1,1} - 12m_{1,1,1,1}, \\
m_{2,2}&=e_2^2 - 2m_{2,1,1} - 6m_{1,1,1,1}, \\
m_{2,1,1}&=e_3e_1 - 4m_{1,1,1,1},\\
m_{1,1,1,1}&=e_4.
\end{align}
$$
Thus our equation becomes, after some more computation,
$$
(\alpha+\beta+\gamma)
(\alpha+\beta+\delta)
(\alpha+\gamma+\delta)
(\beta+\gamma+\delta)=
e_2e_1^2 - e_3e_1 + e_4 = qp^2 - (pq)p + 1=1,
$$
where your equations $p=e_1$, $q=e_2$, $pq=e_3$ and $1=e_4$ were substituted.
This is of course relatively laborious, but has the advantage of being a method that works for arbitrary symmetric polynomials of $\alpha,\beta,\gamma,\delta$.
The coefficients for the expression of products of elementary symmetric polynomials in terms of minimal ones are given by small combinatorial counting problems: the number of matrices with entries $0$ or $1$ with column sums given by the indices $i$ of the factors $e_i$, and row sums given by the indices of the $m_\lambda$. For instance the coefficient $5$ of $m_{2,1,1}$ in $e_2e_1^2$ equals the number of such matrices with row and column sums both given by the sequence $2,1,1$.