Apologies: The following answer is not complete and may contain inaccuracies (although I have tried to identify potential problem areas). I do not have time due to work commitments to see this through to the end but thought given the amount I have already done it would be foolish not to post it in case it is of some help. Best of luck solving this, hopefully @whuber has something up his sleeve.
Here is how I would attempt to solve the problem:
As pointed out, the $X_{(i)}$ are not independent, since they must be ordered, so a simple convolution does not work. We will work with your ordering and hopefully @whuber will point out how much simpler it is with the alternative ordering and how to apply the results from one ordering to the other.
We must therefore calculate the PDF of $Y = \frac{1}{M}\sum_{i=1}^{M}{X_{(i)}}$ by first considering the CDF of $Y$:
$$
F_Y(y)=\Pr(Y\le y)=\int_0^\infty{\cdots\int_0^\infty{f_{X_{(1)},\ldots,X_{(M)}}(x_{(1)},\ldots,x_{(M)})I(x_{(1)}+x_{(2)}+\ldots+x_{(M)}\le yM)dx_{(1)}}\ldots dx_{(M)}}
$$
Where $I(\cdot)$ is an indicator function taking the value $1$ if the inequality holds and $0$ otherwise.
For now I think it is acceptable to integrate only over $x_{(1)}$ to $x_{(M)}$ but it may ultimately be necessary to integrate over all $N$.
Now $f_{X_{(1)},\ldots,X_{(M)}}(x_{(1)},\ldots,x_{(M)})$ is the formula you give above or zero if the $x_{(i)}$ are not ordered. To avoid making separate case statements we add another indicator function:
$$
F_Y(y)=\int_0^\infty{\cdots\int_0^\infty{f_{X_{(1)},\ldots,X_{(M)}}(x_{(1)},\ldots,x_{(M)})I\left(\sum_{i=1}^M{x_{(i)}}\le yM\right)I(x_{(1)}\ge x_{(2)}\ge \ldots \ge x_{(M)})dx_{(1)}}\ldots dx_{(M)}}
$$
Now we work on adjusting the integral limits in such a way to exclude the values of $x_{(i)}$ for which the indicator functions give $0$. Firstly we take the indicator function for the ordering. It will take the value $0$ if $x_{(1)}<x_{(2)}$ so we adjust the integral limits for $x_{(1)}$ to be $(x_{(2)},\infty)$. Likewise the limits for $x_{(2)}$ become $(x_{(3)},\infty)$ and so on. The limits for $x_{(M)}$ are a tricky one because if $M<N$ then a lower limit of $0$ is clearly problematic. I'll park this issue for now if you don't mind!
The other indicator function will set our upper limits. For $x_{(i)}, i=1,\ldots,M-1$ we have the inequality $x_{(i)} \le yM - \sum_{j=i+1}^M{x_{(j)}}$ and for $x_{(M)}$ the upper limit is simply $yM$. In fact, I think you can make it even simpler and use the inequality of $x_{(i)} \le yM - x_{(i+1)}$ and devolve the issue of maintaining the larger inequality "up the food chain" as it were.
Our final integral looks like this:
$$
F_Y(y)=\int_0^{yM}{\int_{x_{(M)}}^{yM-x_{(M)}}{\cdots\int_{x_{(2)}}^{yM-x_{(2)}}{f_{X_{(1)},\ldots,X_{(M)}}(x_{(1)},\ldots,x_{(M)})dx_{(1)}}\ldots dx_{(M-1)}}dx_{(M)}}
$$
Now consider that many of the terms in $f_{X_{(1)},\ldots,X_{(M)}}(x_{(1)},\ldots,x_{(M)})$ can actually be pulled out of the integral or at least distributed in such a way to make it less scary:
$$
F_Y(y)=\frac{\lambda^M N!}{(N-M)!}\int_0^{yM}{(1-e^{-\lambda x_{(M)}})^{N-M}e^{-\lambda x_{(M)}}\int_{x_{(M)}}^{yM-x_{(M)}}{e^{-\lambda x_{(M-1)}}\cdots\int_{x_{(2)}}^{yM-x_{(2)}}{e^{-\lambda x_{(1)}} dx_{(1)}}\ldots dx_{(M-1)}}dx_{(M)}}
$$
(You may want to check whether the limits on the product might actually be $i=1,\ldots,M-1$ and I have assumed you mean $x_{(i)}$ instead of $x_i$ throughout).
The first integral is solved as follows:
$$
\int_{x_{(2)}}^{yM-x_{(2)}}{e^{-\lambda x_{(1)}} dx_{(1)}} = \left[-\frac{e^{-\lambda x_{(1)}}}{\lambda}\right]_{x_{(2)}}^{yM-x_{(2)}} = \frac{e^{-\lambda x_{(2)}}-e^{-\lambda (yM-x_{(2)})}}{\lambda}
$$
Leaving us with:
$$
F_Y(y)=\frac{\lambda^M N!}{(N-M)!}\int_0^{yM}{(1-e^{-\lambda x_{(M)}})^{N-M}e^{-\lambda x_{(M)}}\int_{x_{(M)}}^{yM-x_{(M)}}{e^{-\lambda x_{(M-1)}}\cdots \int_{x_{(3)}}^{yM-x_{(3)}}{e^{-\lambda x_{(2)}}\frac{e^{-\lambda x_{(2)}}-e^{-\lambda (yM-x_{(2)})}}{\lambda}dx_{(2)}}\ldots dx_{(M-1)}}dx_{(M)}}
$$
We can then simplify the innermost integrand:
$$
e^{-\lambda x_{(2)}}\frac{e^{-\lambda x_{(2)}}-e^{-\lambda (yM-x_{(2)})}}{\lambda} = \frac{e^{-2\lambda x_{(2)}}-e^{-\lambda yM}}{\lambda}
$$
This integral can also be solved explicitly and so on.
As I said above, sorry the answer is not complete. Hopefully after doing a couple of further integrations you could use induction to finish the job.
