Frequency domain filtering: must filter coefficients always be conjugate symetric?
I noticed that the coefficients, complex numbers, are not conjugate symetric.

Check what is the frequency associated with each coefficient, usually the coefficients are defined such that $a_i$ corresponds to the frequency $f_s i / N$, where $f_s$ is the sampling frequency, and $N$ is the number of coefficients. So what you have is $a_i = conj(a_{N-i})$, in some cases they may be circularly shifted, and so on.

the signal (audio_filtered) should again be real?

Yes, if you think in the time domain, the convolution involves only real values, so the result is a real signal. If you think in the frequency domain, both the filter and the signal have Hermitian symmetry, so the product of them will have hermitian symmetry as well.
Proof, let the signal with FFT $X$, and the filter with FFT $H$, the FFT of the filtered signal $Y$ is
$Y_{i} = X_i H_i = conj(X_{N-i})conj(H_{N-i}) = conj(X_{N-i}H_{N-i}) = Y_{N-i}

If the input filter is IIR, are things different with respect to conjugate symetry?

No, if you look at the fourier transform of an infinite signal
$$X(\omega) = \sum_{k=0}^{\infty} u_k exp(-2j\pi k \omega)$$
Assume $u_k$ is real, and will will notice that since $e^{-j x}$ is the complex conjugate of $e^{+jx}$, each term will be conjugated and so, $X(\omega)$ is the complex conjugate of $X(-\omega)$