I am talking about the angular momentum. Is there any deep reason? For some other group, it is not the case?
In addition, compactness of SO(3) is crucial. Compactness implies that the eigenvalues of $J_\pm$ under $J_z$ in the adjoint representation ($J_z \circ J_\pm = [J_z,J_\pm]=\pm J_\pm$) must be real (in order for all one-parameter subgroups to form closed loops, it is necessary for their action to be unitary). Since the eigenvalues of the raising and lowering operators (which generalize to a 'root diagram' for arbitrary compact simple Lie groups) are used to project representations out of tensor products exactly like the SO(3) Clebsch Gordan procedure, the elements of the generalized Clebsch-Gordan transition matrices should be real for arbitrary compact groups, with the appropriate initial choice of phase.
For non-compact Lie groups, things are more complicated. You are allowed to have non-unitary group actions, for example, in the action of the Lorentz group on normal 4 vectors, and the representations aren't constrained by the requirement for global consistency (or rather, the requirement for global consistency is less strict).
I have read through various sources and my own notes, based on Quantum Mechanics by Robinette, page 493 and it seems that the fact that they are taken as real is either assumed, or explained as outlined below.
The transformation coefficients $\langle j_1 ,m_1 ; j_2 ,m_2 |j,m; j_1 , j_2\rangle $ are known as the Clebsch-Gordon (CG) coefficients (or the vector coupling coefficients).
The Clebsch-Gordan matrix is unitary (since it just transforms a vector from one basis to another) and by convention its elements are chosen real because the phase of this ket $|j,m; j_1, j_2 \rangle $is arbitrary.
This follows because the matrix elements of the ladder operators $L_\pm$ are chosen to be real. It is possible to choose CGs to be real whenever this observation on the matrix elements of ladder operators holds.
It is possible to have real CGs even for non-compact algebras: an elementary example is $su(1,1)$. Decomposing the tensor product of two irreps in the positive discrete series gives an infinite sum of irreps in the positive discrete series, for which basis states can be expressed as product states using real Clebsch's. Again this is possible because the matrix elements of ladder operators of $su(1,1)$ can be chosen to be real.
Roughly speaking, the "deep" reason is as follows. Since the highest (or lowest) weight of an irrep must be killed by raising (or lowering) operators, this higest state will be a real linear combination of states whenever the raising operators have real matrix elements. Once the highest states is constructed, the other states in the irrep are reached by applying lowering operators, which will generate real combination of states states if their matrix elements are real.
Note that the phase of the raising and lowering operators of $su(2)$ or $su(1,1)$ need not be real in principle, but chosen to be real because if ensuing simplifications.