Another method is to estimate the time-invariant coefficients in a second stage equation, using the mean error as the dependent variable.
First, estimate the model with FE. From here you get an estimation of $\beta$ and $\gamma_{t}$. For simplicity, let's forget about the year-effects. Define the estimation error $\hat{u}_{it}$ as before:
$$ \hat{u}_{it} \equiv y_{it} - X_{it}\hat{\beta} $$
The linear predictor $\bar{u}_{i}$ is:
$$ \bar{u}_{i} \equiv \frac{\sum_{t=1}^{T}\hat{u}_{i}}{T} = \bar{y_{it}} - \bar{x}_{i}\hat{\beta} $$
Now, consider the following second stage equation:
\begin{equation}
\bar{u}_{i} = \delta male_{i} + c_{i}
\end{equation}
Assuming that gender is uncorrelated with unobserved factors $c_{i}$. Then, the OLS estimator of $\delta$ is unbiased and time-consistent (this is, it is consistent when $T \rightarrow \infty$).
To prove the above, replace the original model into the estimator $\bar{u}_{i}$:
$$ \bar{u}_{i} = \bar{x}_{i}\beta - \bar{x}_{i}\hat{\beta} + \delta male_{i} + c_{i} + \frac{\sum_{t=1}^{T}\epsilon_{it}}{T} $$
The expectation of this estimator is:
$$ E(\bar{u}_{i}) = \bar{x}_{i}\beta - \bar{x}_{i}E(\hat{\beta}) + \delta male_{i} + E(c_{i}) + \frac{\sum_{t=1}^{T}E(\epsilon_{it})}{T} $$
If assumptions for FE consistency hold, $\hat{\beta}$ is an unbiased estimator of $\beta$, and $E(\epsilon_{it}) = 0$. Thus:
$$ E(\bar{u}_{i}) = \delta male_{i} + E(c_{i}) $$
This is, our predictor is an unbiased estimator of the time-invariant components of the model.
Regarding consistency, the probability limit of this predictor is:
$$ p \lim\limits_{T \rightarrow \infty} \bar{u}_{i} = p \lim\limits_{T \rightarrow \infty} \left( \bar{x}_{i}\beta\right) - p \lim\limits_{T \rightarrow \infty} \left(\bar{x}_{i}\hat{\beta}\right) + p \lim\limits_{T \rightarrow \infty} \delta male_{i} + p \lim\limits_{T \rightarrow \infty} c_{i} + p \lim\limits_{T \rightarrow \infty} \left( \frac{\sum_{t=1}^{T}\epsilon_{it}}{T}\right) $$
Again, given FE assumptions, $\hat{\beta}$ is a consistent estimator of $\beta$, and the error term converges to its mean, which is zero. Therefore:
$$ p \lim\limits_{T \rightarrow \infty} \bar{u}_{i} = \delta male_{i} + c_{i} $$
Again, our predictor is a consistent estimator of the time-invariant components of the model.