Why are Lagrangians linear in $\dot{q}$ so ubiquitous? Gauge theory, Berry phase, Dirac Equation, and more

Let us start with a general remark. Why there are typically at most only first-order derivatives in the Lagrangian (density) is discussed in e.g. this Phys.SE post. This implies that the Euler-Lagrange EL equations are at most of second order, cf. e.g. this Phys.SE post.

Now let us return to OP's question. OP is interested in the case where the Lagrangian (density) is affine in the time derivatives. This is quite common. It has some interesting consequences:

• The EL equations are at most of first order.

• The main example is the Hamiltonian formulation: $L_H(q,\dot{q},p,t) ~=~\sum_{i=1}^n p_i \dot{q}^i - H(q,p,t).$ (This formula can be generalizes to field theory.)

• Given a Lagrangian (density) affine in time derivatives, if we try to construct the corresponding Hamiltonian formulation via a Legendre transformation following the Dirac-Bergmann analysis, we encounter primary constraints.

• Faddeev & Jackiw devised another method to construct a Hamiltonian formulation, see e.g. arXiv:hep-th/9306075. This is related to presymplectic geometry, cf. e.g. this Phys.SE post.

• For concrete examples of such systems, see e.g. this, this, this & this Phys.SE posts.

Let me discuss just one aspect of your question. I don't understand the statement about "first-order nature of the Dirac equation". Note that the Dirac equation is a system of four first-order partial differential equations (PDEs) for four components of the Dirac spinor. However, it is well-known that any system of PDEs can be rewritten as a system of first order PDEs. Furthermore, the Dirac equation in electromagnetic field is generally equivalent to one fourth-order equation for just one component (see my article http://akhmeteli.org/wp-content/uploads/2011/08/JMAPAQ528082303_1.pdf (J. Math. Phys. 52, 082303 (2011))).