In Quantum field theory by M. Srednicki, we find on p. 337 the statement regarding the Maxwell Lagrangian
The action we seek should be Lorentz invariant, gauge invariant, parity and time-reversal invariant, and no more than second order in derivatives. The only candidate is $S=\int d^4x\mathcal{L}$, where $$\mathcal{L}=-\frac{1}{4}F^{\mu\nu}F_{\mu\nu}+J^\mu A_\mu.$$
Let us ignore the source term for now and expanding the field-strength tensor $F^{\mu\nu}$, we find $$ \mathcal{L} = - \frac{1}{2} \left(\partial^\mu A^\nu\right) \left(\partial_\mu A_\nu\right) + \frac{1}{2} \left(\partial^\mu A^\nu\right) \left(\partial_\nu A_\mu\right) .\tag{1} $$
How do I arrive at eq. (1) from first principles?
We know that the action of a relativistic field theory must be Lorentz-invariant scalar. We can construct Lorentz-invariant scalars by contracting Lorentz-invariant tensors.
From polarization experiments, we should be able to conclude that photons are spin-1 particles. Hence, we are looking at vector fields $A^\mu$.
Photons should be massless otherwise the QED Lagrangian is not invariant under local gauge transformations which would then spoil electric charge conservation. Hence, our Lagrangian cannot contain a mass term and the Lagrangian should be invariant under local gauge transformations.
Photons don't self-interact directly (though they can interact through vacuum fluctuations). Hence, there are no power terms, e.g., $(A_\mu A^\mu)^2$, in the Lagrangian.
With the assumptions so far, the action contains only derivatives of $A^\mu$. From dimensional analysis, see Constraints on scalar field theories in $n$ dimensional spacetime, we reduce the Lagrangian to $$ \mathcal{L} = c_1 \left( \partial_\mu A^\nu \right) \left( \partial^\mu A_\nu \right) + c_2 \left( \partial_\mu A^\mu \right) \left( \partial_\nu A^\nu \right) + c_3 \left( \partial_\mu A^\nu \right) \left( \partial_\nu A^\mu \right) \tag{2}. $$
For the action to exist, we need the vector field $A^\mu$ to be integrable, i.e., to falloff sufficiently rapid at infinity. Hence, boundary terms are zero and we can perform partial integration to move around the derivatives. We thereby see that we can rewrite the third term as the first term. The argument was used in C. de Rham. Massive Gravity. 2014. and we reduced the Lagrangian density to $$ \mathcal{L} = c_1 \left( \partial_\mu A^\nu \right) \left( \partial^\mu A_\nu \right) + c_2 \left( \partial_\mu A^\mu \right) \left( \partial_\nu A^\nu \right) \tag{3}. $$
Comparing eq. (3) with eq. (1), we are left to show that $c_1=-1/2$ and $c_2=-c_1$.
We still have not used the time-reversal and parity invariance.
How do we see that time-reversal and parity invariance require $c_1=-c_2=-1/2$?
According to this post on physicsforums.com, parity transformations only change the field argument, i.e.,
$$ \hat{A}(x,t) \xrightarrow{x\to -x} \hat{A}(-x,t) \tag{4} $$
but inserting eq. (4) into the action and substituting the integration variable to "undo" the parity transformation does not change anything (or I missed something).
The same would be true for the time-reversal.