According to Hund's first rule, a set of degenerate orbitals are singly occupied first, before the second slot in any of the orbitals are populated. This is quite intuitive because electron-electron repulsions would make an atom more unstable if the electrons start filling two at a time in a single orbital.
However, the rule also states that in the ground state, electrons that singly fill the orbitals will adopt the same spin. Why is it necessary for the electrons to have the same spin?
Disclaimer
My following answer is the "traditional" explanation of Hund's first rule, which is based on a smaller value of $V_\mathrm{ee}$ (electron-electron repulsions) in the triplet state arising from Fermi holes. According to Levine's Quantum Chemistry 7th ed.:
This traditional explanation turns out to be wrong in most cases. It is true that the probability that the two electrons are very close together is smaller for the helium $^3S$ $1s2s$ term than for the $^1S$ $1s2s$ term. However, calculations with accurate wave functions show that the probability that the two electrons are very far apart is also less for the $^3S$ term. The net result is that the average distance between the two electrons is slightly less for the $^3S$ term than for the $^1S$ term, and the interelectronic repulsion is slightly greater for the $^3S$ term.
The calculations show that the $^3S$ term lies below the $^1S$ term because of a substantially greater electron nucleus-attraction in the $^3S$ term as compared with the $^1S$ term. Similar results are found for terms of the atoms beryllium and carbon.
Wildcat has an answer which briefly mentions this. Aside from that Levine gives several references:
With that said, here is the old-fashioned argument for why parallel spins are favoured over paired spins (all else being equal).
The origin of Hund's first rule lies in "exchange energy", which is a way of saying that electrons with like spin repel each other less than electrons with unlike spin.
To see why, let's consider an excited state of the helium atom: $\mathrm{1s^1 2s^1}$. (The ground state configuration $\mathrm{1s^2}$ is of no use to us because those electrons have to be paired with unlike spin.) In the excited state, the electrons can either have parallel spin (the triplet case) or paired spin (the singlet case).
The total electronic wavefunction comprises a spatial part, which describes which orbitals the electrons are in (i.e. it specifies $n,l,m_l$), and a spin part, which describes the, well, spin ($m_s$). According to the Pauli exclusion principle, the total wavefunction must be antisymmetric with respect to interchange of the electron labels: $\Psi(1,2) = -\Psi(2,1)$.
The spin wavefunctions are simultaneous eigenfunctions of the $\hat{S}^2$ and $\hat{S}_z$ operators. The total spin quantum number, $S$, can be found from the eigenvalue of the operator $\hat{S}^2$:
$$\hat{S}^2\psi = \hbar^2 S(S+1)\psi$$
For more detailed discussion, refer to any QM text on the coupling of angular momenta.
There are three possible triplet spin wavefunctions, which are degenerate.
$$\psi_{\mathrm{spin,triplet}}(1,2) = \begin{cases} \alpha(1)\alpha(2) \\ \beta(1)\beta(2) \\ \frac{1}{\sqrt{2}}[\alpha(1)\beta(2) + \alpha(2)\beta(1)] \end{cases}$$
It turns out that $\hat{S}^2\psi_\mathrm{triplet} = 2\hbar^2\psi_{\mathrm{triplet}}$, which means that these have the quantum number $S = 1$. This corresponds to a multiplicity of $M = 2S + 1 = 3$, hence "triplet". These wavefunctions represent the case where the electrons have parallel spin. (Loosely speaking, each electron has a spin of $1/2$, and if these are aligned parallel the total spin is $1/2 + 1/2 = 1$.)
The singlet state is
$$\psi_{\mathrm{spin,singlet}}(1,2) = \frac{1}{\sqrt{2}}[\alpha(1)\beta(2) - \alpha(2)\beta(1)]$$
and as you might guess, $\hat{S}^2\psi_\mathrm{singlet} = 0$, meaning that $S = 0$ and $M = 2S+1 = 1$. This represents the case where the electrons are paired - one can think of it as the two $1/2$ electrons spins “cancelling each other out”.
The important thing to notice is that the triplet spin wavefunctions are symmetric with respect to permutation of the labels, whereas the singlet spin wavefunction is antisymmetric. For example:
$$\begin{align} \psi_\mathrm{singlet}(2,1) &= \frac{1}{\sqrt{2}}[\alpha(2)\beta(1) - \alpha(1)\beta(2)] \\ &= -\frac{1}{\sqrt{2}}[\alpha(1)\beta(2) - \alpha(2)\beta(1)] \\ &= -\psi_\mathrm{singlet}(1,2) \end{align}$$
The only way of constructing these are the linear combinations below.
$$\begin{align} \psi_\mathrm{space,symm} &= \frac{1}{\sqrt{2}}[\mathrm{1s}(1)\mathrm{2s}(2) + \mathrm{1s}(2)\mathrm{2s}(1)] \\ \psi_\mathrm{space,antisymm} &= \frac{1}{\sqrt{2}}[\mathrm{1s}(1)\mathrm{2s}(2) - \mathrm{1s}(2)\mathrm{2s}(1)] \end{align}$$
Remember that the atomic orbitals above are functions of the coordinates of the electron: $\mathrm{1s}(1) = \mathrm{1s}(\vec{r}_1) = \mathrm{1s}(r_1,\theta_1, \phi_1)$. Now, the antisymmetric spatial wavefunction has an interesting property. Let's explore what happens when the electrons get very close together, i.e. $\vec{r}_1 = \vec{r}_2 = \vec{r}$.
$$\begin{align} \psi_\mathrm{space,antisymm} &= \frac{1}{\sqrt{2}}[\mathrm{1s}(\vec{r})\mathrm{2s}(\vec{r}) - \mathrm{1s}(\vec{r})\mathrm{2s}(\vec{r})] \\ &= 0 \end{align}$$
If the wavefunction is 0 when the two electrons are on top of each other, that means that there is no probability of this happening. The average distance between the two electrons $\left\langle\left|\vec{r}_2 - \vec{r}_1\right|\right\rangle$ will therefore be larger, and on average the electrons will repel each other less, since they are on average further away from each other. This is known as a Fermi hole.
In fact, one can also show that the symmetric wavefunction has a larger probability amplitude when $\vec{r}_1 = \vec{r}_2$. This is known as a Fermi heap.
As I mentioned earlier, the Pauli exclusion principle dictates that the total wavefunction be antisymmetric. This means that the triplet spin wavefunctions, which are symmetric, must be paired with the antisymmetric spatial wavefunction. Likewise, the singlet spin wavefunction, which is antisymmetric, must be paired with the symmetric spatial wavefunction.
$$\begin{array}{c|c|c} \text{Spin} & \text{Spatial} & \text{Total}\\ \hline \text{Triplet (S)} & \text{AS} & \text{S } \times \text{ AS} = \text{AS} \\ \text{Singlet (AS)} & \text{S} & \text{AS } \times \text{ S} = \text{AS} \end{array}$$
So, there is a decrease in electron-electron repulsion when the electrons are parallel - because this means they have a triplet spin wavefunction, and hence an antisymmetric spatial wavefunction, which has a Fermi hole.
There is another, perhaps more mathematical, way of looking at it, and that is to calculate the expectation value of the electronic energy for both wavefunctions. For the helium atom, under the assumption of an infinitely heavy nucleus (essentially an atomic version of the Born-Oppenheimer approximation), the total Hamiltonian (in atomic units) is
$$\hat{H} = -\frac{1}{2}(\nabla_1^2+\nabla_2^2) - \frac{2}{r_1} - \frac{2}{r_2} + \frac{1}{r_{12}}$$
where $r_i$ is the distance of electron $i$ from the nucleus, and $r_{12}$ is the distance between the two electrons. The expectation value of the energy is simply $\left\langle\Psi\middle|\hat{H}\middle|\Psi\right\rangle$ (assuming our wavefunction is already normalised). It turns out that the only important term in our present discussion is the electron-electron repulsion term $1/r_{12}$.
Let's look at the triplet spin wavefunction (of course paired with the antisymmetric spatial wavefunction) first. The operator doesn't depend on spin degrees of freedom, so we can separate those, and the inner product of the spin wavefunction with itself is just unity.
$$\begin{align} \left\langle\Psi_\mathrm{triplet}\middle|\frac{1}{r_{12}}\middle|\Psi_\mathrm{triplet}\right\rangle &= \left\langle\psi_\mathrm{spin,triplet}\middle|\psi_\mathrm{spin,triplet}\right\rangle\left\langle\psi_\mathrm{space,antisymm}\middle|\frac{1}{r_{12}}\middle|\psi_\mathrm{space,antisymm}\right\rangle \\ &= \frac{1}{2} \left\langle \mathrm{1s}(1)\mathrm{2s}(2) - \mathrm{1s}(2)\mathrm{2s}(1) \middle| \frac{1}{r_{12}} \middle| \mathrm{1s}(1)\mathrm{2s}(2) - \mathrm{1s}(2)\mathrm{2s}(1) \right\rangle \\ &= \frac{1}{2} \Biggr\{ \left\langle \mathrm{1s}(1)\mathrm{2s}(2) \middle|\frac{1}{r_{12}}\middle| \mathrm{1s}(1)\mathrm{2s}(2) \right\rangle \\ &\quad\quad\quad - \left\langle \mathrm{1s}(2)\mathrm{2s}(1) \middle|\frac{1}{r_{12}}\middle| \mathrm{1s}(1)\mathrm{2s}(2) \right\rangle \\ &\quad\quad\quad - \left\langle \mathrm{1s}(1)\mathrm{2s}(2) \middle|\frac{1}{r_{12}}\middle| \mathrm{1s}(2)\mathrm{2s}(1) \right\rangle \\ &\quad\quad\quad + \left\langle \mathrm{1s}(2)\mathrm{2s}(1) \middle|\frac{1}{r_{12}}\middle| \mathrm{1s}(2)\mathrm{2s}(1) \right\rangle\Biggr\} \end{align}$$
Ugly, but that's what it is.
We can simplify it however by noting that the first and the fourth brackets are the same (since the variables in the integral are just dummy variables, interchanging the electron labels does not affect the value of the integral). This integral is traditionally called a Coulomb integral, since it simply reflects the average Coulombic repulsion between two independent (uncorrelated) electrons, one in the 1s orbital and the second in the 2s orbital. It is given the symbol $J_\mathrm{1s2s}$ (the subscripts indicating the orbitals that the electrons are in). Because this is an expectation value for a potential energy arising from repulsion, it is necessarily a positive quantity.
The second and third terms are actually also the same, again because the variables are just dummy variables and you can interchange the labels. This integral, however, does not have a physical interpretation; it is a purely quantum mechanical effect. It is called the exchange integral, and is given the symbol $K_\mathrm{1s2s}$. It's a bit harder to show, but the exchange integral is also a positive quantity (there is a question on it here). However, its magnitude is not as large as $J_\mathrm{1s2s}$.
Putting it all together (leaving out subscripts for clarity) we have
$$\left\langle\Psi_\mathrm{triplet}\middle|\frac{1}{r_{12}}\middle|\Psi_\mathrm{triplet}\right\rangle = \frac{1}{2}(J-K-K+J) = J-K$$
and if you were to perform the same calculation on the singlet spin wavefunction (not very difficult, just replace the minus signs with a plus!) you would get
$$\left\langle\Psi_\mathrm{singlet}\middle|\frac{1}{r_{12}}\middle|\Psi_\mathrm{singlet}\right\rangle = \frac{1}{2}(J+K+K+J) = J+K$$
Since $J > K > 0$, the expectation value of the energy for the singlet state is higher.
Well, you just have to do a ton more ugly maths and perhaps a few hundred more integrals. But the idea is the same. For every pair of electrons with parallel spins, there is a corresponding exchange integral that reduces the expectation value of the energy.
Therefore, in a $\mathrm{3d^5}$ manganese atom (we ignore the closed shells) there are
$$\frac{5!}{2!(5-2)!} = 10$$
different pairs of d-d electrons with parallel spins. Therefore if you evaluate the expectation value of the energy, there will be a $-10K_\mathrm{3d3d}$ term, which is stabilising. This is precisely the exchange energy effect that Geoff describes in his answer.
The lowest energy state has parallel spins to maximize the exchange energy.
As you say, there's a Coulomb repulsion between two electrons to put them in the same orbital.
There's also a quantum mechanical effect. The exchange energy (which is favorable) increases with the number of possible exchanges between electrons with the same spin and energy.
Going from the top state to the middle state, we remove the Coulomb repulsion between electrons in the same orbital.
Going between the middle to the bottom (most stable state and that predicted by Hund's rule), we gain the exchange energy, because these two electrons are indistinguishable.
Some understanding can be gained by looking at the symmetry of the orbital parts of the wavefunctions involved.
The total wavefunction for electrons must be anti-symmetric with respect to exchanging electrons, the Pauli Principle. The total wavefunction has spin and spatial parts, $\psi =\psi _{space}\psi _{spin}$.
When electrons are paired with total spin = 0, (a singlet state) the spatial part of the wavefunction is symmetric and the spin part therefore antisymmetric. Its the other way round for the triplet state, spins parallel and total spin = 1, $\psi_{triplet} =\psi _{space}(antisym)\psi _{spin}(sym)$.
The effect is that the electrons occupy different regions of space in triplet states; they have a zero probability of being in the same place.
In singlets the electrons crowd together; formally there is a non-zero probability of being in the same region of space. This is why the electrostatic repulsion, which gives energy $e^2/r_{1,2}$ and is large when the electrons are near to one another, raises the energy of the singlet. Thus the degeneracy between singlet (spin paired) and triplet (spins 'parallel') is split.
Therefore the energy difference between singlet and triplet ( and more generally paired and unpaired electrons) is a consequence of the wavefunction symmetry and so the value of the total spin.
The picture shows the effect on the spatial part of wavefunctions illustrated with 2D 'particle in a box' wavefunctions.
With respect to the coulomb and exchange integrals and energies, K and J respectively, both are positive, the total energy works out as $K+J$ and $K-J$ which shows the energy splitting.
(For completeness both Coulomb and Exchange integrals are forms of electrostatic repulsion. If W is the potential energy, which may contain several terms positive and negative of the form $e^2/r$ , then the Coulomb integral with electrons 1 and 2 and orbitals a and b, is
$\int \psi(1,a)\psi(2,b)~W~ \psi(1,a)\psi(2,b)$
and exchange integral is
$\int \psi(2,a)\psi(1,b)~W~ \psi(1,a)\psi(2,b)$
where the electrons are 'exchanged')
