Typo.
There is a typo in the 2006 version of Bishop's PRML, which has been corrected in later versions.
Instead, equation (9.39) on p443 should read "the expectation of the indicator variable $z_{nk}$ under this posterior distribution is then given by:
$$\begin{align}
\mathbb{E}[z_{nk}] &= \frac{\sum_{\mathbf{z}_n} z_{nk} \prod_{k'} \left[\pi_{k'} \mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_{k'}, \boldsymbol{\Sigma}_{k'}) \right]^{z_{nk'}}}{\sum_{\mathbf{z}_n} \prod_{j} \left[\pi_{j} \mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_{j}, \boldsymbol{\Sigma}_{j}) \right]^{z_{nj}}} \\
&= \frac{\pi_{k} \mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)}{\sum^K_{j=1} \pi_j \mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)} = \gamma(z_{nk})
\end{align}$$
If you want further confirmation that this correction is indeed the case, see p19 of the publisher errata on the 2006 version of Bishop here.
Derivation.
When I worked through this I found it tricky getting my head around the indicators and the selection of the appropriate dimensionality of the posterior for computing expectations. Here is one way of deriving the result.
The graphical model for the Gaussian mixture model, from Bishops's PRML is the following:

Notice that if we are interested in computing the posterior distribution $p(\mathbf{z}_n | \mathbf{x}_n, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})$, this amounts to reversing the direction of the arrow so that it goes from $\mathbf{x}_n$ to $\mathbf{z}_n$. The conditional independencies in the diagram also imply that we can factorise the posterior over all latent variables $\mathbf{Z}$ into a product of $n$ independent posterior distributions:
$$p(\mathbf{Z} | \mathbf{X}, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})= \prod^N_{n=1} p(\mathbf{z}_n | \mathbf{x}_n, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})$$
We are now interested in computing the expectation of the indicator variable $z_{nk}$, that is, $\mathbb{E}[z_{nk}]$, where the expectation is with respect to the posterior $p(\mathbf{Z} | \mathbf{X}, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})$.
Now for clarity, I am going to explicitly choose to work at the dimensionality of vectors that is, rewriting everything in terms of $\mathbf{z}_n \in \mathbb{R}^K$, rather than matrices $\mathbf{Z}$, or scalars $z_{nk}$.
The indicator variable $z_{nk}$ takes the value 1 when the $n$th latent mixture assignment $\mathbf{z}_n$ is component $k$, and $0$ otherwise. Explicitly writing the indicator variable in terms of an indicator function on $\mathbf{z}_n$, we have that $z_{nk} = \mathbb{I}(\mathbf{z}_n = k)$. Due to the factorisation of the posterior over the latent variable matrix $\mathbf{Z}$ into $n$ independent posterior distributions over vectors $\mathbf{z}_n$, the expectation simplifies so that we only need to consider it with respect to the posterior $p(\mathbf{z}_n | \mathbf{x}_n, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})$:
$$\begin{align}
\mathbb{E}_{p(\mathbf{Z} | \mathbf{X}, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})}[z_{nk}] &= \mathbb{E}_{p(\mathbf{Z} | \mathbf{X}, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})}[\mathbb{I}(\mathbf{z}_n = k)] \\
&= \mathbb{E}_{\prod^N_{n=1} p(\mathbf{z}_n | \mathbf{x}_n, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})}[\mathbb{I}(\mathbf{z}_n = k)] \\
&= \mathbb{E}_{p(\mathbf{z}_n | \mathbf{x}_n, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})}[\mathbb{I}(\mathbf{z}_n = k)] \\
&= p(\mathbf{z}_n = k | \mathbf{x}_n, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi}) \\
\tag{1}
\end{align}$$
Where in the last line we have used the fact that the expectation of an indicator function is a probability. As you have sought further clarity on the steps behind this derivation, I have put further exposition at the bottom of the post.
Using standard probability formulae, the posterior probability is:
$$p(\mathbf{z}_n | \mathbf{x}_n, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi}) = \frac{p(\mathbf{x}_n, \mathbf{z}_n)}{p(\mathbf{x}_n)} = \frac{p(\mathbf{x}_n, \mathbf{z}_n)}{\sum_{\mathbf{z}_n} p(\mathbf{x}_n, \mathbf{z}_n)}$$
Using the factorisation induced by the graphical rerpesentation of the Gaussian mixture model, and also the distributions specified by the model in equation $(9.10)$ and $(9.11)$, the joint distribution is:
$$\begin{align}
p(\mathbf{x}_n, \mathbf{z}_n) &= p(\mathbf{z}_n) \cdot p(\mathbf{x}_n | \mathbf{z}_n) \\
&= \prod^K_{k'=1} \pi_k^{z_{nk'}} \cdot \prod^K_{k'=1} \mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_{k'}, \boldsymbol{\Sigma}_{k'})^{z_{nk'}} \\
&= \prod^K_{k'=1} [\pi_{k'} \mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_{k'}, \boldsymbol{\Sigma}_{k'})]^{z_{nk'}}
\end{align}$$
The marginal distribution is:
$$p(\mathbf{x}_n) = \sum_{\mathbf{z}_n} p(\mathbf{x}_n, \mathbf{z}_n) = \sum_{\mathbf{z}_n} \prod^K_{j=1} [\pi_k \mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)]^{z_{nj}} = \sum^K_{j=1} \pi_j \mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j) $$
Hence the posterior probability distribution is:
$$p(\mathbf{z}_n | \mathbf{x}_n, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi}) = \frac{\prod^K_{k'=1} [\pi_{k'} \mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_{k'}, \boldsymbol{\Sigma}_{k'})]^{z_{nk'}}}{\sum^K_{j=1} \pi_j \mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)}$$
Putting this all together, and using the above distribution to query the probability that the $n$th latent mixture assignment is component $k$, we have:
$$\mathbb{E}[z_{nk}] = p(\mathbf{z}_n = k | \mathbf{x}_n, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi}) = \frac{\pi_k \mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k})}{\sum^K_{j=1} \pi_j \mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)} = \gamma(z_{nk})$$
Further derivation of $(1)$ in response to comments.
This derivation relies on two results that are quite useful, which I will not prove formally with any rigour, but state.
a. Computing expectations with respect to a factorisable joint distribution.
If a joint distribution over a random vector $\mathbf{x} \in \mathbb{R}^N$ can be factorised into a product of marginal distributions over each of its $N$ components, that is, $p(\mathbf{x}) = p(X_1, ..., X_N) = \prod^N_{i=1} p(X_i)$, and we wish to compute the expectation of a deterministic function of one of the components with respect to the joint distribution, then we only need to compute the expectation with respect to the relevant marginal distribution. To see this:
$$\begin{align}
\mathbb{E}_{p(\mathbf{x})}[f(X_i)] &= \int_{\mathbb{R}^N} f(x_i) p(\mathbf{x}) d\mathbf{x} \\
&= \int ... \int f(x_i) p(x_1, ..., x_N) dx_1 ... dx_N \\
&= \int ... \int f(x_i) \prod^N_{i=1} p(x_i) dx_1 ... dx_N \\
&= \int ... \int f(x_i) p(x_i) \prod_{-i} p(x_{-i}) dx_i dx_{-i} \\
&= \int ... \int \prod_{-i} p(x_{-i}) \int f(x_i) p(x_i) dx_i dx_{-i} \\
&= \mathbb{E}_{p(x_i)}[f(X_i)]\int ... \int \prod_{-i} p(x_{-i}) dx_{-i} \\
&= \mathbb{E}_{p(x_i)}[f(X_i)]
\end{align}$$
Where $\int ... \int d_{x_{-i}}$ refers to the multiple integrals with respect to all components other than component $i$, forming a total of $(N-1)$ integrals. And where I have assumed that all distributions are appropriately normalised - meaning that the rightmost term in the penultimate equality is equal to 1.
This argument is being generalised to the case where the posterior distribution over a matrix $\mathbf{Z}$, $p(\mathbf{Z} | \mathbf{X}, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})$ can be factorised into $N$ independent posterior distributions over vectors $\mathbf{z}_n$, $p(\mathbf{z}_n | \mathbf{x}_n, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})$. Hence
$$\mathbb{E}_{p(\mathbf{Z} | \mathbf{X}, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})}[f(\mathbf{z}_n)] = \mathbb{E}_{\prod^N_{n=1} p(\mathbf{z}_n | \mathbf{x}_n, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})}[f(\mathbf{z}_n)] = \mathbb{E}_{p(\mathbf{z}_n | \mathbf{x}_n, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})}[f(\mathbf{z}_n)]$$
Normally you would use the unconscious/lazy statistician rule to evaluate this, but in the case that $f(\cdot)$ is an indicator function, then we use the next result.
b. Expectation of an indicator function.
There is a fairly standard result that the expectation of an indicator function of a random variable is its probability. For a set $A$ and random variable $X$:
$$\mathbb{E}[\mathbb{I}(X \in A)] = 1 \cdot P(X \in A) + 0 \cdot P(X \notin A) = P(X \in A)$$
Applying this result, we have that:
$$\begin{align}\mathbb{E}_{p(\mathbf{z}_n | \mathbf{x}_n, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})}[z_{nk}] &= \mathbb{E}_{p(\mathbf{z}_n | \mathbf{x}_n, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})}[\mathbb{I}(\mathbf{z}_n = k)] \\
&= 1 \cdot p(\mathbf{z}_n = k | \mathbf{x}_n, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi}) + 0 \cdot p(\mathbf{z}_n \neq k | \mathbf{x}_n, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})\\
&= p(\mathbf{z}_n = k | \mathbf{x}_n, \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})
\end{align}$$
Lastly, if you find this confusing, you can think of $z_{nk}$ as a Bernoulli random variable. The mean of a Bernoulli random variable is the parameter $p$, which in this case is the posterior probability that $\mathbf{z}_n$ is component $k$.
Hopefully that addresses the difficulties you are experiencing.