New observational data of neutron stars since GW170817 have helped improve our knowledge about nuclear symmetry energy especially at high densities. We have learned particularly: (1) The slope parameter $L$ of nuclear symmetry energy at saturation density $rho_0$ of nuclear matter from 24 new analyses is about $Lapprox 57.7pm 19$ MeV at 68% confidence level consistent with its fiducial value, (2) The curvature $K_{rm{sym}}$ from 16 new analyses is about $K_{rm{sym}}approx -107pm 88$ MeV, (3) The magnitude of nuclear symmetry energy at $2rho_0$, i.e. $E_{rm{sym}}(2rho_0)approx 51pm 13$ MeV at 68% confidence level, has been extracted from 9 new analyses of neutron star observables consistent with results from earlier analyses of heavy-ion reactions and the latest predictions of the state-of-the-art nuclear many-body theories, (4) while the available data from canonical neutron stars do not provide tight constraints on nuclear symmetry energy at densities above about $2rho_0$, the lower radius boundary $R_{2.01}=12.2$ km from NICERs very recent observation of PSR J0740+6620 of mass $2.08pm 0.07$ $M_{odot}$ and radius $R=12.2-16.3$ km at 68% confidence level sets a tight lower limit for nuclear symmetry energy at densities above $2rho_0$, (5) Bayesian inferences of nuclear symmetry energy using models encapsulating a first-order hadron-quark phase transition from observables of canonical neutron stars indicate that the phase transition shift appreciably both the $L$ and $K_{rm{sym}}$ to higher values but with larger uncertaintie , (6) The high-density behavior of nuclear symmetry energy affects significantly the minimum frequency necessary to rotationally support GW190814s secondary component of mass (2.50-2.67) $M_{odot}$ as the fastest and most massive pulsar discovered so far.