Dark matter subhalos are the remnants of (incomplete) halo mergers. Identifying them and establishing their evolutionary links in the form of merger trees is one of the most important applications of cosmological simulations. The Hierachical Bound-Tracing (hbt) code identifies halos as they form and tracks their evolution as they merge, simultaneously detecting subhalos and building their merger trees. Here we present a new implementation of this approach, hbt+, that is much faster, more user friendly, and more physically complete than the original code. Applying hbt+ to cosmological simulations we show that both the subhalo mass function and the peak-mass function are well fit by similar double-Schechter functions.The ratio between the two is highest at the high mass end, reflecting the resilience of massive subhalos that experience substantial dynamical friction but limited tidal stripping. The radial distribution of the most massive subhalos is more concentrated than the universal radial distribution of lower mass subhalos. Subhalo finders that work in configuration space tend to underestimate the masses of massive subhalos, an effect that is stronger in the host centre. This may explain, at least in part, the excess of massive subhalos in galaxy cluster centres inferred from recent lensing observations. We demonstrate that the peak-mass function is a powerful diagnostic of merger tree defects, and the merger trees constructed using hbt+ do not suffer from the missing or switched links that tend to afflict merger trees constructed from more conventional halo finders. We make the hbt+ code publicly available.