You know the maximum possible digit sum (ignoring the prime condition for a moment) is $36$, due to the digitsum from $9999$ within your range.
The possible digitsums under that barrier that are divisible by $8$ (ignoring $0$) are $8, 16, 24, 32$.
How many ways can you write each of these as a sum of $4$ digits where each digit falls within $0 \leq d \leq 9$?
If you were working in a team and had to do this by hand, you could probably divide-and-conquer and have each person tackle a different digitsum in order to enumerate all the cases more quickly.
However, you already know the $24$ digitsum is a waste of time, since any number whose digitsum is divisible by $24$ will also be divisible by $3$, which means the number itself will be divisible by $3$ and therefore not prime.
Then you can eliminate the cases that can't possibly be rearranged to form a prime number (such as those missing a $1, 3, 7$, or $9$ to form the last digit of the prime).
This leaves you with the following cases:
$8 = [[0, 0, 1, 7], [0, 0, 3, 5], [0, 1, 1, 6], [0, 1, 2, 5], [0, 1, 3, 4], [0, 2, 3, 3], [1, 1, 1, 5], [1, 1, 2, 4], [1, 1, 3, 3], [1, 2, 2, 3]]$
$16 = [[0, 0, 7, 9], [0, 1, 6, 9], [0, 1, 7, 8], [0, 2, 5, 9], [0, 2, 7, 7], [0, 3, 4, 9], [0, 3, 5, 8], [0, 3, 6, 7], [0, 4, 5, 7], [1, 1, 5, 9], [1, 1, 6, 8], [1, 1, 7, 7], [1, 2, 4, 9], [1, 2, 5, 8], [1, 2, 6, 7], [1, 3, 3, 9], [1, 3, 4, 8], [1, 3, 5, 7], [1, 3, 6, 6], [1, 4, 4, 7], [1, 4, 5, 6], [1, 5, 5, 5], [2, 2, 3, 9], [2, 2, 5, 7], [2, 3, 3, 8], [2, 3, 4, 7], [2, 3, 5, 6], [3, 3, 3, 7], [3, 3, 4, 6], [3, 3, 5, 5], [3, 4, 4, 5]]$
$32 = [[5, 9, 9, 9], [6, 8, 9, 9], [7, 7, 9, 9], [7, 8, 8, 9]]$
Then for each unique permutation of each case that falls within your acceptable range (also ending in $1, 3, 7$, or $9$), test if the number is actually prime. There are only $276$ numbers to test.
One way you can do this is by listing the (relevant) primes under $\sqrt{10000} = 100$:
$7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97$
Then you test each one for divisibility into your candidate number. If they all fail the test, then the candidate number is prime.
I assume you can at least use a calculator for this step, because otherwise this is probably too tedious a problem to do by hand. The size of your team would need to be large enough to be able to perform all the primality-check calculations in the time allotted.
Another approach (assuming you had a large enough team) is to divide the $1000-10000$ range up into chunks and have each person enumerate their section of the range and then use the segmented Sieve of Eratosthenes approach using the primes under $100$ to cross out the composite numbers. Then you can check the digitsums of the resulting primes directly.
Best Answer
I have a proposal for the solution of the more general problem for the sum of like powers of arbitrary natural exponents. This is based on the Faulhaber-/Bernoulli-polynomials.
Consider that polynomials organized in a matrix (which I called "Gp" when I found them). Each row r contain the coefficients for the r 'th Faulhaber polynomial. I introduce also the notation for a "Vandermondevector" -type: $ V(x) = \text{columnvector}(1,x,x^2,x^3,...)$; if I use it as diagonalmatrix I prefix it with a "d" like $\ ^dV(x)$.
The matrix $Gp$ begins like $$ \Tiny \begin{bmatrix} 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1/2 & 1/2 & 0 & 0 & 0 & 0 \\ 0 & 1/6 & 1/2 & 1/3 & 0 & 0 & 0 \\ 0 & 0 & 1/4 & 1/2 & 1/4 & 0 & 0 \\ 0 & -1/30 & 0 & 1/3 & 1/2 & 1/5 & 0 \\ 0 & 0 & -1/12 & 0 & 5/12 & 1/2 & 1/6 \end{bmatrix}$$ and along the rows we see the coefficients for the integrals of the Bernoulli-polynomials. Using that in a matrix-formula we get $$ Gp \cdot V(m) = V(1) + V(2) + V(3) + ... + V(m) $$ To have the sum $V(3)+V(6)+...+V(3m)$ instead it suffices to write $$ ( ^dV(3) \cdot Gp \cdot \ ^dV(1/3))\cdot V(3m) = V(3) + V(6) + V(9) + ... + V(3m) $$ Thus, introducing a matrix-valued-function $$H(x) \overset{\text{def}}= \ ^dV(x) \cdot Gp \cdot \ ^dV(1/x) $$ allows to write $$ \begin{eqnarray} H(3) \cdot V(3m) &=& V(3)+V(6)+...+V(15)+ ... +V(3m) \\ H(5) \cdot V(5m) &=& V(5)+V(10)+...V(15) + ... +V(5m) \\ H(15) \cdot V(15m) &=& V(15)+V(30)+... +V(15m) \\ \end{eqnarray} $$ and finally, with a notation for the highest integer equal or below n divisible by m for convenience $$ n:m \overset{\text{def}}{=} m \cdot \Big\lfloor { n \over m} \Big\rfloor $$ we can write $$ S_{3,5}(m) = H(3)\cdot V( m:3) + H(5) \cdot V(m:5) - 2H(15) \cdot V(m:15 ) $$ For $m=1000$ we get the vector of solutions for the sums with exponents 0,1,2,3,4,5,6,7 $$ \begin{matrix} S_{3,5}(1000) &=& \small \begin{bmatrix}\begin{array} {rll} 401 &\Tiny = 3^0+6^0+...+999^0 + 5^0+10^0+...+1000^0 -2(15^0+30^0+...+990^0)\\ 201003 &\Tiny = 3+6+...+999 + 5+10+15+...+1000 -2(15+30+...+990)\\ 134335661 &\Tiny = 3^2+6^2+...+999^2 + 5^2+10^2+...+1000^2 -2(15^2+30^20+...+990^2)\\ 101003482917 \\ 81004632555017 \\ 67672443055260693 \\ 58149771588796814081 \\ 51008046700741091931597 \end{array} \end{bmatrix} \end{matrix}$$ The formulae are $$ \begin{matrix} S_{3,5}&=& \small \begin{bmatrix} \begin{array} {l} \frac 13 m_3+\frac 15 m_5- \frac 2{15} m_{15} \\ \frac 1 6 m_3^2 +\frac 12 m_3 + \frac 1{10} m_5^2 +\frac 12 m_5 - \frac 1{15} m_{15}^2- m_{15} \\ \frac 1 9 m_3^3+ \frac 12 m_3^2+\frac 12 m_3 +\frac 1{15}m_5^3+\frac 12 m_5^2+ \frac 56m_5 -\frac 2{45} m_{15}^3-m_{15}^2-5 m_{15} \\ \vdots \end{array} \end{bmatrix} \end{matrix}$$ where the $m_k$ denote the $(m:k) = k \cdot \lfloor \frac mk \rfloor$.
You're interested in the second row $S_{3,5;1}$.
With $m=1000 $ inserted you get $$\small \frac 16 m_3^2+\frac 12 m_3+ \frac 1{10}m_5^2+\frac 12 m_5-\frac 1{15} m_{15}^2-m_{15} $$ and with $m_3=999,m_5=1000,m_{15}=990$ inserted you get $$ S_{3,5;1}(1000) = \frac 12(\frac 1 3 999^2+ 999 + \frac 1 5 1000^2+ 1000 -\frac 2{15} 990^2-2 \cdot 990) $$ which is finally $$ S_{3,5;1}(1000) = 201003 $$
One can arrive at the last formula for the second row by far simpler means (as shown in other answers); however for the general problem: to get the rows for the sums-of-like-powers with higher exponents, I think this (known) array of Faulhaber-polynomials (handled by the matrix-formulae) should be the general and most convenient one.
(as someone remarked already: this (without the removal of the multiples of 15) is also the "project-euler-problem-1", but much generalized)