  Measurements of proton and nuclear collisions at the Large Hadron Collider at
nucleon-nucleon c.m. energies up to $\sqrt{s_{NN}}=$ 13 TeV, have improved our
understanding of hadronic interactions at the highest energies reached in
collisions of cosmic rays with nuclei in the earth atmosphere, up to
$\sqrt{s_{NN}}\approx 450$ TeV. The Monte Carlo event generators (EPOS, QGSJET,
and SIBYLL) commonly used to describe the air showers generated by
ultrahigh-energy cosmic rays (UHECR, with $E_{CR}\approx 10^{17}$--$10^{20}$
eV) feature now, after parameter retuning based on LHC Run-I data, more
consistent predictions on the nature of the cosmic rays at the tail of the
measured spectrum. However, anomalies persist in the data that cannot be
accommodated by the models. Among others, the total number of muons (as well as
their maximum production depth) remains significantly underestimated
(overestimated) by all models. Comparisons of EPOS, QGSJET, and SIBYLL
predictions to the latest LHC data, and to collider MC generators such as
PYTHIA, indicate that improved description of hard multiple minijet production
and nuclear effects may help reduce part of the data--model discrepancies, shed
light on the UHECR composition approaching the observed $E_{CR}\approx 10^{20}$
eV cutoff, and uncover any potential new physics responsible of the observed
anomalies.
